I had a look at the OGR source, and I think option 2 applies.
#ifdef notdef
if (poSRS)
{
const char* pszAuthority = poSRS->GetAuthorityName(NULL);
const char* pszAuthorityCode = poSRS->GetAuthorityCode(NULL);
if (pszAuthority != NULL && pszAuthorityCode != NULL && strcmp(pszAuthority, "EPSG") == 0)
{
json_object* poObjCRS = json_object_new_object();
json_object_object_add(poObjCRS, "type", json_object_new_string("name"));
json_object* poObjProperties = json_object_new_object();
json_object_object_add(poObjCRS, "properties", poObjProperties);
if (strcmp(pszAuthorityCode, "4326") == 0)
{
json_object_object_add(poObjProperties, "name", json_object_new_string("urn:ogc:def:crs:OGC:1.3:CRS84"));
}
else
{
/* FIXME?: the issue is that for geographic SRS, OGR will expose a latitude/longitude axis order */
/* which is probably not what is written in the file! */
json_object_object_add(poObjProperties, "name", json_object_new_string(CPLSPrintf("urn:ogc:def:crs:EPSG::%s", pszAuthorityCode)));
}
const char* pszCRS = json_object_to_json_string( poObjCRS );
VSIFPrintfL( fpOut_, "\"crs\": %s,\n", pszCRS );
json_object_put(poObjCRS);
}
}
#endif
That #ifdef notdef was added by this commit but I couldn't say why.
Perhaps you can try rebuilding gdal without the #ifdef and #endif?
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 ] }
....
Best Answer
Looks like ogr2ogr has support for exactly this problem.
The below is copied directly from that page:
How do I flip coordinates when they are not in the expected order
The EPSG has a recommanded order for geographic SRS where the coordinates tuples of a geometry must appear in the (latitude, longitude) order, whereas most GIS will properly display such geometries if they appear in the (longitude, latitude) order. This issue can be often encountered in situation with GML3 files, WFS 1.1 data, etc... that adhere to the (latitude, longitude) order.
When, for some reason, the coordinate order isn't the one that is wished, the following trick can be used to flip them:
... An alternative way of achieving the same result, providing using GDAL 2.0 compiled with Spatialite support to benefit from SQLite SQL dialect, is: