ST_Transform in GIS support

[David Orme]

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 coordinate transformation available 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 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 it 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

[Leonel Câmara]

That looks great to me!

[David Orme]

I guess more generally there is an issue with having to do a round trip through WKT to update a record - it loses precision and is there is an explicit warning about using it for passing data:

http://postgis.org/docs/ST_AsText.html

I don't see any point in altering how the select represents the record (and WKT is more interpretable) but it would be better if you could pass a WKB straight back to the record. The update method avoids the problem but I think it might be fixable by updating the representer here:

https://github.com/web2py/pydal/blob/master/pydal/representers/postgre.py#L22

If the first character of the value is 0 (zero) then the code could attempt to treat the value as WKB and just return the value (possibly with SRID) rather than having to put it through ST_GeomFromText(). The backend validates WKB inputs, so you can't just pass in any old nonsense starting with a zero.

[David Orme]

So, in addition to st_transform, I would add st_aswkb. This isn't an official Postgis command - it provides a method to bypass the automatic conversion to WKT found in here:

https://github.com/web2py/pydal/blob/master/pydal/adapters/base.py#L550

With a suitably modified PostgreRepresenter, that then allows this:

from pydal import geoPoint, geoLine,

geoPolygon

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)}])

# [1L, 2L]

db(db.transform_test).select().as_list()
# [{'id': 1L, 'point': 'POINT(0 0)', 'point_wec': None},
#  {'id': 2L, 'point': 'POINT(3 0)', 'point_wec': None}]

copy_wkb = db(db.transform_test.id == 1).select(db.transform_test.point.st_aswkb().with_alias('wkb')).first()

rec = db.transform_test[2]
rec.update_record(point=copy_wkb['wkb'])
# <Row {'point_wec': None, 'id': 2L, 'point': '0101000020E610000000000000000000000000000000000000'}>

db(db.transform_test).select().as_list()
# [{'id': 1L, 'point': 'POINT(0 0)', 'point_wec': None},
#  {'id': 2L, 'point': 'POINT(0 0)', 'point_wec': None}]

This avoids the round trip through WKT. 

Note also that the extended WKB provided by PostGIS contains a reference to the geometry SRID, so trying to update a field with a different SRID gets caught, as you might hope.

In [10]: rec.update_record(point_wec=copy_wkb['wkb'])
---------------------------------------------------------------------------
ProgrammingError                          Traceback (most recent call last)
# <snip>
ProgrammingError: ('ERROR', '22023', 'Geometry SRID (4326) does not match column SRID (4087)')

I don't think this breaks any existing functions. 

Cheers,

David