Spatial / GIS support in DAL

[DenesL]

Spatial / GIS support (in latest trunk)
=====================
Sponsored by AidIQ for use by Sahana Eden

available for (so far)
++++++++++++++++++++++

1) Postgres + PostGIS
http://postgis.refractions.net/docs/
2) MS SQL
http://msdn.microsoft.com/en-us/library/ff848797.aspx

SQLite + Spatialite next.

Both geometry and geography fields are supported.

supported functions (so far)
++++++++++++++++++++++++++++

st_asgeojson (PostGIS only)
st_astext
st_contained
st_contains
st_distance
st_equals
st_intersects
st_overlaps
st_simplify (PostGIS only)
st_touches
st_within


Some MS SQL examples
++++++++++++++++++++

dbm = DAL(r"mssql://user:pass@host:db")

sp = dbm.define_table('spatial',
Field('geo1','geometry()')
)

# a point
sp.insert(geo1="POINT (1 2)")
1
# a line
sp.insert(geo1="LINESTRING (100 100, 20 180, 180 180)")
2
# a polygon (a square in this case)
sp.insert(geo1="POLYGON ((0 0, 150 0, 150 150, 0 150, 0 0))")
3

rr=dbm(sp.id>0).select()
print rr
spatial.id,spatial.geo1
1,
2,
3,

# why doesn't it display geo1?
# field geo1 has an undisplayable internal representation in MS SQL
# e.g. the first record's POINT(1 2) is stored as
# 0x00000000010C000000000000F03F0000000000000040
# to visualize as WKT use function st_astext

rr=dbm(sp.id>0).select(sp.id, sp.geo1.st_astext())
print rr
spatial.id,spatial.geo1.STAsText()
1,POINT (1 2)
2,"LINESTRING (100 100, 20 180, 180 180)"
3,"POLYGON ((0 0, 150 0, 150 150, 0 150, 0 0))"

# or using an alias
rr=dbm(sp.id>0).select(sp.id, sp.geo1.st_astext().with_alias('g1'))
for r in rr: print r.spatial.id, r.g1
...
1 POINT (1 2)
2 LINESTRING (100 100, 20 180, 180 180)
3 POLYGON ((0 0, 150 0, 150 150, 0 150, 0 0))

# other functions

# STContains
rr=dbm(sp.geo1.st_contains("POINT(1
1)")).select(sp.id,sp.geo1.st_astext())
print rr
spatial.id,spatial.geo1.STAsText()
3,"POLYGON ((0 0, 150 0, 150 150, 0 150, 0 0))"

# STDistance
rr=dbm(sp.id>0).select(sp.id, sp.geo1.st_distance("POINT(-1
2)").with_alias('d'))
for r in rr: print r.spatial.id, r.d
...
1 2.0
2 140.714249456
3 1.0

# STIntersects
rr=dbm(sp.geo1.st_intersects("LINESTRING(20 120,60
160)")).select(sp.id,sp.geo1.st_astext())
print rr
spatial.id,spatial.geo1.STAsText()
2,"LINESTRING (100 100, 20 180, 180 180)"
3,"POLYGON ((0 0, 150 0, 150 150, 0 150, 0 0))"


Some PostGIS examples
+++++++++++++++++++++

#pg is the postgres db connection

pg.define_table('franchise',
Field('code', 'string', length=1),
Field('name'),
)

pg.define_table('fastfood',
Field('franchise', 'reference franchise'),
Field('lat', 'double'),
Field('lon', 'double'),
Field('geom', 'geometry()'),
)

# parameters for geom fields are: schema, srid and dimension
# e.g. Field('geom', 'geometry('',4326,2)')
# hard-coded defaults: srid=4326, dimension=2
# migration works, fields are added/dropped as required.

pg.franchise.insert(code='b', name='Burger Joint')

pg.fastfood.insert(franchise=1,lat=25.8092,lon=-80.24,geom='POLYGON((0
0, 10 0, 10 10, 0 10, 0 0))')
pg.commit()

ff=pg.fastfood
q=ff.geom.st_overlaps('POLYGON((1 1,5 1,5 5,1 5,1 1))')
print pg(q).select()
fastfood.id,fastfood.franchise,fastfood.lat,fastfood.lon,fastfood.geom
# no records, there are no overlaps

q=ff.geom.st_overlaps('POLYGON((1 1,11 1,11 11,11 1,1 1))')
print pg(q).select()
fastfood.id,fastfood.franchise,fastfood.lat,fastfood.lon,fastfood.geom
1,1,25.8092,-80.24,0103000020E61000000100000005000000000000000000000000000000000
00000000000000000244000000000000000000000000000002440000000000000244000000000000
00000000000000000244000000000000000000000000000000000

print
pg(ff.id>0).select(ff.id,ff.franchise,ff.lat,ff.lon,ff.geom.st_simplify(1).st_astext())
fastfood2.id,fastfood2.franchise,fastfood2.lat,fastfood2.lon,"ST_AsText(ST_Simplify(fastfood2.geom,
1.0))"
1,1,25.8092,-80.24,"POLYGON((0 0,10 0,10 10,0 10,0 0))"

# geography example
pg.define_table('airport',
Field('code','string',3),
Field('geog','geography()'),
)
a=pg.airport

a.insert(code='LAX', geog='POINT(-118.4079 33.9434)')
1
a.insert(code='CDG', geog='POINT(2.5559 49.0083)')
2
a.insert(code='REK', geog='POINT(-21.8628 64.1286)')
3

lax="POINT(-118.4079 33.9434)"
rr=pg(a.id>0).select(a.code,a.geog.st_distance(lax).with_alias('from_lax'))
for r in rr: print r.airport.code, r.from_lax
...
LAX,0.0
CDG,9124665.26917
REK,6969660.54628
# results are in meters

rr=pg(a.id>0).select(a.code,a.geog.st_asgeojson().with_alias('geojson'))
for r in rr: print r.airport.code, r.geojson
...
LAX {"type":"Point","coordinates":[-118.4079,33.943399999999997]}
CDG {"type":"Point","coordinates":[2.5559,49.008299999999998]}
REK {"type":"Point","coordinates":[-21.8628,64.128600000000006]}


More to come...

Denes Lengyel.

[Massimo Di Pierro]

This is a amazing! Thanks Denes.

[Manuele Pesenti]

Hi DanesL,
many compliments for the very interesting work... I would like to ask if
is it possible to insert some spatial data in the exadecimal format? I
think it could be more OO in the code to manage object like the one I
can get from the ppygis library from whom it seams you cannot simply get
the postgres like text format you used for insertion. The exadecimal
format can be easily get by the write_ewkb() method.

Thank you very mutch

Manuele

[DenesL]

Hi Manuele,

you are right about PPyGIs, from their caveats:
"Only EWKB representations are supported — the WKB, EWKT and WKT
alternatives are not."

At the moment, the DAL can only insert WKT data, as required by the
sponsor.
But in the future we should support the other formats as well.

Regards,
Denes

[Manuele Pesenti]

Il 21/03/2012 01:24, DenesL ha scritto:
> Hi Manuele,
>
> you are right about PPyGIs, from their caveats:

> "Only EWKB representations are supported � the WKB, EWKT and WKT

> alternatives are not."
>
> At the moment, the DAL can only insert WKT data, as required by the
> sponsor.
> But in the future we should support the other formats as well.
>
> Regards,
> Denes

if you consider what I said in a recent recipe[1] I post on web2py[:]
defining the geometry field (in a raw way indeed) as 'text' the EWKB
format for insertion is naturaly accepted and work fine so maybe it
would be quite easy to support EWKB format for insertion. Or can you
suggest some code or library to convert into WKT?

thank you very mutch

cheers
Manuele

http://www.web2pyslices.com/slice/show/1487/a-gis-recipe

[DenesL]

Could you provide more info on what are you trying to accomplish?
where is your data and how it is stored, etc.

The current DAL GIS implementation takes advantage of the built-in
support provided in the DB for GIS data, be it PostGIS or MSSQL
spatial support, without any additional dependencies in web2py. The
input is provided in human-readable WKT text form, the output is more
flexible.

WKT is not a technical limitation, it is just the sponsor's
preference. Other formats are equally possible with some additional
modification to the DAL's code.

On Mar 21, 4:18 am, Manuele Pesenti <manuele.pese...@gmail.com> wrote:
> Il 21/03/2012 01:24, DenesL ha scritto:> Hi Manuele,
>
> > you are right about PPyGIs, from their caveats:

> > "Only EWKB representations are supported the WKB, EWKT and WKT

[Manuele Pesenti]

Il 21/03/2012 15:37, DenesL ha scritto:
> Could you provide more info on what are you trying to accomplish?
> where is your data and how it is stored, etc.
>
> The current DAL GIS implementation takes advantage of the built-in
> support provided in the DB for GIS data, be it PostGIS or MSSQL
> spatial support, without any additional dependencies in web2py. The
> input is provided in human-readable WKT text form, the output is more
> flexible.
>
> WKT is not a technical limitation, it is just the sponsor's
> preference. Other formats are equally possible with some additional
> modification to the DAL's code.

Hi DenesL,
actually I've not yet tested the DAL GIS implementation, I just read
your examples and I imagined to have import a shape file. Some libraries
like osgeo can read this kind of data but as long as I have looked for I
have not found a way to export geometries in WKT format. How did you
approach this problem?

[DenesL]

What is in the shape file?.

[DenesL]

If you are referring to an ESRI shapefile then this might help:

http://gis.stackexchange.com/questions/15920/how-to-convert-shapefile-geometries-to-wkb-using-ogr

[Manuele Pesenti]

Il 21/03/2012 19:39, DenesL ha scritto:
> If you are referring to an ESRI shapefile then this might help:
>
> http://gis.stackexchange.com/questions/15920/how-to-convert-shapefile-geometries-to-wkb-using-ogr

Hi,

thanks for the very usefull link... but the object seams not to answer
the quest. For inserting geometries in database using new DAL features
you need WKT format and not WKB... every library I found support
exporting WKB (osgeo.ogr, Shapely, ppygis) but none seams to export WKT
format. On the other hand a good point of your approach is that
OpenLayers support WKT if you need to insert geometries by drawing them
on a map.

cheers

M.

[DenesL]

You were too quick to dismiss it.
There is also an exportToWkt in OGR:
http://www.gdal.org/ogr/classOGRGeometry.html

On Mar 21, 5:35 pm, Manuele Pesenti <manuele.pese...@gmail.com> wrote:
> Il 21/03/2012 19:39, DenesL ha scritto:> If you are referring to an ESRI shapefile then this might help:
>

> >http://gis.stackexchange.com/questions/15920/how-to-convert-shapefile...

[Fran]

On Wednesday, 21 March 2012 21:35:32 UTC, Manuele wrote:

For inserting geometries in database using new DAL features
you need WKT format and not WKB... every library I found support
exporting WKB (osgeo.ogr, Shapely, ppygis) but none seams to export WKT
format.

Shapely can export WKT:

http://toblerity.github.com/shapely/manual.html#well-known-formats

F 

[Tim Michelsen]

Am 20.03.2012 01:24, schrieb DenesL:
> Spatial / GIS support (in latest trunk)
> =====================
> Sponsored by AidIQ for use by Sahana Eden

Very cool. Great news.

Has the Sahana Eden team also implemented a solution to the proxy issue
with the JS mapping frameworks such as openlayers?

beginnings of a proxy for web2py. See
http://trac.sahanapy.org/wiki/BluePrintGISProxy for more
https://gist.github.com/284772

proxy issue in (python) web frameworks
http://lists.osgeo.org/pipermail/openlayers-users/2009-April/011443.html

Kind regards,
Timmie

[Fran]

On Monday, 26 March 2012 22:05:14 UTC+1, Timmie wrote:

Am 20.03.2012 01:24, schrieb DenesL:
> Spatial / GIS support (in latest trunk)
> =====================
> Sponsored by AidIQ for use by Sahana Eden
Very cool. Great news.

Has the Sahana Eden team also implemented a solution to the proxy issue
with the JS mapping frameworks such as openlayers?

beginnings of a proxy for web2py. See
http://trac.sahanapy.org/wiki/BluePrintGISProxy for more

Wow, old URL!

Yes, we have a proxy:

https://github.com/flavour/eden/blob/master/controllers/gis.py#L2897

F

[Tim Michelsen]

> Wow, old URL!
Where's your current roadmap?

Maybe I go and pickup something.

[Nate]

I cannot get Postgis to work with

web2py™     Version 2.6.3-stable+timestamp.2013.09.15.17.01.20
Python     Python 2.7.3: C:\Python27\python.exe (prefix: C:\Python27)

pg.define_table('franchise',
  Field('code', 'string', length=1),
  Field('name'),
)

pg.define_table('fastfood',
  Field('franchise', 'reference franchise'),
  Field('lat', 'double'),
  Field('lon', 'double'),
  Field('geom', 'geometry()'),
)

results in

<class 'psycopg2.ProgrammingError'> function addgeometrycolumn(unknown, unknown, unknown, integer, unknown, integer) does not exist LINE 1: SELECT AddGeometryColumn ('', 'fastfood', 'geom', 4326, 'GEO... ^ HINT: No function matches the given name and argument types. You might need to add explicit type casts.

Assistance appreciated

[Niphlod]

this is rather a "how can I enable GIS support on postgres?" question...

http://trac.osgeo.org/postgis/wiki/UsersWikiNewbieAddgeometrycolumn

[Nate]

Absolutely! Thank you very much.