Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question
    
Looks like the basic cause of this is a bug in ogr2ogr - there is no CRS (coordinate reference system) element included in the GeoJSON output. Apparently fixed in 1.10.0, but I only have 1.9.2. trac.osgeo.org/gdal/ticket/4995 –  Steve Bennett Apr 14 at 1:34

2 Answers 2

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 ] }
 ....
share|improve this answer
    
Thank you! I didn't know about st_asgeojson - that will probably work, although one hiccup is there are several line segments per osm_id, which ogr2ogr automatically merges (or at least combines into one file). I still don't quite understand why ogr2ogr is producing projected coordinates which, say, Github.io is incapable of rendering. –  Steve Bennett Apr 13 at 22:15

Gene's answer pointed me in the right direction. The simple answer seems to be that the data in PostGIS is in Web Mercator (EPSG 3786) and so if I want lat/long it needs to be re-projected (or de-projected...), with an SRID like EPSG 4326.

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], "-t_srs", "EPSG: 4326"])

enter image description here

EDIT Turns out you need to explicitly state the source SRS as well:

call (["ogr2ogr", "-f", "GeoJSON", record[1] + ".json",  'PG:dbname=\'gis\'', "-sql", 'SELECT ...', "-t_srs", "EPSG: 4326", "-s_srs", "EPSG:3857"])
share|improve this answer
    
Actually this is still wrong - if I'd even looked at the screenshot I was taking I'd notice. They're lat-longs, but way off. (The left end of the orange line in the above pic should be exactly at Skipton, top left.) –  Steve Bennett Apr 14 at 5:32
    
Don't use the subprocess module to run ogr2ogr, use the psycopg2 module (How to change the SRID of existing data in PostGIS?) or the osgeo/osr module in Python. –  gene Apr 14 at 16:08

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.