[GIS] PostGIS to GeoJSON with GDAL: how to get lat/longs

gdalgeojsonlatitude longitudeogr2ogrpostgis

I'm trying to convert line features from PostGIS to GeoJSON files, as follows:

from psycopg2 import *
from subprocess import call

conn = connect(dbname='gis')
cur=conn.cursor()

cur.execute("SELECT distinct osm_id,route_name from planet_osm_line where route_name like '%Rail Trail' and network='rcn';")

for record in cur:
  print "%s" % record[1]
  call (["ogr2ogr", "-f", "GeoJSON", record[1] + ".json",  'PG:dbname=\'gis\'', "-sql", 'SELECT route_name,osm_id,tags::hstore->\'state\' as state,way from planet_osm_line where osm_id=%d' % record[0]])

That is, ultimately ogr2ogr gets called like this:

ogr2ogr -f GeoJSON blah.json "PG:dbname='gis'" -sql 'SELECT route_name,osm_id,tags::hstore->\'state\' as state,way from planet_osm_line where osm_id=-12345'

The GeoJSON files that get created have coordinates that look like this [ 16153119.78, -4561568.23 ]:

{ "type": "Feature", "properties": { "route_name": "Waverley Rail Trail", "osm_id": -64660, "state": null }, "geometry": { "type": "LineString", "coordinates": [ [ 16153119.78, -4561568.23 ], [ 16153117.1, -4561567.48 ], [ 16153114.81, -4561565.86 ], [ 16153113.24, -4561563.55 ], [ 16153112.64, -4561561.63 ], [ 16153112.64, -4561558.82 ], [ 16153113.6, -4561556.2 ], [ 16153115.39, -4561554.05 ], [ 16153117.82, -4561552.66 ], [ 16153119.98, -4561552.2 ], [ 16153122.76, -4561552.48 ], [ 16153125.28, -4561553.69 ], [ 16153127.35, -4561555.88 ], [ 16153128.43, -4561558.46 ], [ 16153128.62, -4561560.38 ], [ 16153128.08, -4561563.13 ], [ 16153126.64, -4561565.51 ], [ 16153124.47, -4561567.26 ], [ 16153121.82, -4561568.17 ], [ 16153119.78, -4561568.23 ] ] } }

Please forgive my ignorance, but what am I doing wrong here? I assume what is being written is some kind of projected value rather than the raw lat/long – so how do I get those?

I've tried using t_srs=ESPG:3857 but that didn't fix it.

Best Answer

OpenLayers uses the EPSG:3857 coordinate system, in meters, and not the WGS84 system, in degrees, look at OpenStreetMap Wiki: EPSG:3857

But why use subprocess and ogr2ogr?

1) you can use directly the PostGIS ST_AsGeoJSON function:

import psycopg2
conn = psycopg2.connect("dbname='osm' host='localhost' user='me'")
cur = conn.cursor()
# srid of the layer
sql = """SELECT Find_SRID('public', 'planet_osm_line', 'way');"""
cur.execute(sql)
print cur.fetchone()[0]
900913 # The srid is 900913 (now ESPG:3857)

What is the srid of your layer ?
Use of ST_AsGeoJSON:

cur = conn.cursor()
sql = """SELECT ST_AsGeoJSON(way) from planet_osm_line WHERE osm_id=34587377;"""
cur.execute(sql)
print cur.fetchone()[0]
{"type":"LineString","coordinates":[[286345.320000000006985,6623598.759999999776483],[286086.950000000011642,6623556.240000000223517],[286256.710000000020955,6622864.799999999813735],[286244.340000000025611,6622861.809999999590218],[286248.78000000002794,6622841.459999999962747],[286195.179999999993015,6622784.349999999627471],[286193.900000000023283,6622756.799999999813735],[286201.429999999993015,6622760.429999999701977],[286214.479999999981374,6622746.240000000223517],[286280.869999999995343,6622702.580000000074506],[286315.070000000006985,6622679.21999999973923],[286298.820000000006985,6622638.849999999627471],[286381.39000000001397,6622609.290000000037253],[286408.419999999983702,6622599.429999999701977],[286499.520000000018626,6622396],[286540.380000000004657,6622248.320000000298023],[286491.299999999988358,6622026.139999999664724],[286333.520000000018626,6621847.75],[285781.320000000006985,6621559.429999999701977],[285733.900000000023283,6621542.389999999664724],[285704.659999999974389,6621539.450000000186265],[285580.510000000009313,6621479.200000000186265],[285461.570000000006985,6621473.490000000223517],[285387.659999999974389,6621462.879999999888241],[285204.580000000016298,6621478.320000000298023],[285177.35999999998603,6621476.370000000111759],[285151.14000000001397,6621468.360000000335276],[285092.150000000023283,6621434.330000000074506],[285063.840000000025611,6621405.099999999627471],[285042.570000000006985,6621364.759999999776483],[284963.869999999995343,6621333.660000000149012],[284913.429999999993015,6621318.639999999664724],[284824.5,6621241.110000000335276],[284719,6621002.940000000409782],[284415.020000000018626,6620781.519999999552965],[283894.53000000002794,6620613.259999999776483],[283610.39000000001397,6620628.200000000186265],[283526.669999999983702,6620603.169999999925494]]}

2) or the osgeo.ogrmodule:

from osgeo import ogr
connString = """PG: host='localhost' dbname='osm' user'=me'"""
conn = ogr.Open(connString)
for layer in conn:
     print layer.GetName()
planet_osm_point
planet_osm_roads
planet_osm_line
planet_osm_polygon
# first layer = planet_osm_point
conn = ogr.Open(connString)
layer= conn[0]
print layer.GetName()
planet_osm_point
print layer.GetSpatialRef().ExportToWkt() 
'PROJCS["Popular Visualisation CRS / Mercator (deprecated)",GEOGCS["Popular Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular Visualisation Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]]'
for feature in layer:
   geom = feature.GetGeometryRef()
   print geom.ExportToJson()
{ "type": "Point", "coordinates": [ 283528.026815425022505, 6604019.682090770453215 ] }
{ "type": "Point", "coordinates": [ 283528.026815425022505, 6604019.664439089596272 ] }
 ....
Related Question