Hi,
I wanted to see if anyone has opinions on adding ST_Transform in to the
list of available GIS commands. I've only looked at this from the point of
view of PostgreSQL and PostGIS, but there is precedent for having Postgre
only GIS functions. The basic argument is that coordinate transformation is
a fundamental part of GIS functionality that would be helpful to expose via
the DAL. In my case, I have users providing coordinates as latitude and
longitude and want to make distance calculations - we have st_distance but
it is essentially useless on geographic coordinates.
We can always use `db.executesql` to run more complex GIS commands
directly, but I don't think it is too hard to make it transformation
through the DAL. I've got some code I can push but I wanted to see if I'm
reinventing the wheel or if there is an implementation issue I haven't
considered.
First, some demo data:
from pydal import geoPoint, geoLine, geoPolygon
# table with two geometry columns, one as WGS84, one as World Equidistant
Cylindrical
trans_test = db.define_table('transform_test',
Field('point', 'geometry()'),
Field('point_wec', 'geometry(public, 4087, 2)'
))
trans_test.bulk_insert([{'point': geoPoint(0,0)}, {'point': geoPoint(3,0)}])
There are some caveats about the current implementation and the GIS data
formats that get passed. Record updates work expect a geometry to be
provided as Well Known Text (WKT) and it seems like the underlying Well
Known Binary (WKB) representation is always translated to WKT when a
geometry field is selected. You might expect the output of these to differ
in the data, but they don't:
db(db.transform_test.id).select(db.transform_test.point).as_list()
#[{'point': 'POINT(0 0)'}, {'point': 'POINT(3 0)'}]
db(db.transform_test).select(db.transform_test.point.st_astext()).as_list()
#[{'_extra': {'ST_AsText("transform_test"."point")': 'POINT(0 0)'}},
# {'_extra': {'ST_AsText("transform_test"."point")': 'POINT(3 0)'}}]
Other functions do return WKB. Compare:
db(db.transform_test).select(db.transform_test.point.st_simplify(0)).as_list
()
#[{'_extra': {'ST_Simplify("transform_test"."point",0.0)':
'0101000020E610000000000000000000000000000000000000'}},
# {'_extra': {'ST_Simplify("transform_test"."point",0.0)':
'0101000020E610000000000000000008400000000000000000'}}]
db(db.transform_test).select(db.transform_test.point.st_simplify(0).
st_astext()).as_list()
#[{'_extra': {'ST_AsText(ST_Simplify("transform_test"."point",0.0))':
'POINT(0 0)'}},
# {'_extra': {'ST_AsText(ST_Simplify("transform_test"."point",0.0))':
'POINT(3 0)'}}]
If we want to update a record with a simplified geometry, then we have to
pass in WKT, so have to do this:
simplified = db(db.transform_test.id == 2).select(db.transform_test.point.
st_simplify(0).st_astext().with_alias('simp')).first()
rec = db.transform_test[2]
rec.update_record(point = simplified['simp'])
# <Row {'point_wec': None, 'id': 2L, 'point': 'POINT(3 0)'}>
But you can do it much more simply using a direct update, because the
underlying WKB data never leaves the backend.
db(db.transform_test.id == 2).update(point = db.transform_test.point.
st_simplify(0))
I'd suggest that an st_transform function follows the model of st_simplify.
Given the existing code as model, I have an implementation that seems to
work:
db(db.transform_test.id).select(db.transform_test.point.st_transform(4087)).
as_list()
#[{'_extra': {'ST_Transform("transform_test"."point",4087)':
'0101000020F70F000000000000000000000000000000000000'}},
# {'_extra': {'ST_Transform("transform_test"."point",4087)':
'0101000020F70F00002589B7E3196214410000000000000000'}}]
db(db.transform_test.id).select(db.transform_test.point.st_transform(4087).
st_astext()).as_list()
#[{'_extra': {'ST_AsText(ST_Transform("transform_test"."point",4087))':
'POINT(0 0)'}},
# {'_extra': {'ST_AsText(ST_Transform("transform_test"."point",4087))':
'POINT(333958.472379821 0)'}}]
I've chosen the WEC projection because that is easy to verify that a point
3 degrees east along the equator from the origin has a X cooordinate in WEC
of 1/120th of the diameter of the Earth given the underlying radius at the
equator in the WGS84 datum:
import math
(6378137 * 2 * math.pi) / 120
333958.4723798207
We can then do either:
transformed = db(db.transform_test.id == 2).select(db.transform_test.point.
st_transform(4087).st_astext().with_alias('wec')).first()
rec = db.transform_test[2]
rec.update_record(point_wec = transformed['wec'])
#<Row {'point_wec': 'POINT(333958.472379821 0)', 'id': 2L, 'point':
'POINT(3 0)'}>
or
db(db.transform_test.id == 2).update(point_wec = db.transform_test.point.
st_transform(4087))
Thoughts? I haven't made a pull request to pydal but I'm happy to if this
seems sensible.
Cheers,
David
--
Resources:
- http://web2py.com
- http://web2py.com/book (Documentation)
- http://github.com/web2py/web2py (Source code)
- https://code.google.com/p/web2py/issues/list (Report Issues)
---
You received this message because you are subscribed to the Google Groups
"web2py-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
For more options, visit https://groups.google.com/d/optout.
