OGR Projection transformation error

Question

I am trying to convert shapefile that was created by FME in EPSG:900913 to EPSG:3857 (I know its the same definition) with OGR in a script.

datasource = ogr.Open(os.path.join(dirname,shapefileName))
layer = datasource.GetLayer(0)
srcSpatialRef = layer.GetSpatialRef()
dstSpatialRef = osr.SpatialReference()
dstSpatialRef.ImportFromEPSG(3857)
coordTransform = osr.CoordinateTransformation(srcSpatialRef, dstSpatialRef)

for i in range(layer.GetFeatureCount()):
    srcFeature = layer.GetFeature(i)
    srcGeometry = srcFeature.GetGeometryRef()
    srcGeometry.Transform(coordTransform)
    geometry = GEOSGeometry(srcGeometry.ExportToWkt())

I got this error:

GDAL_ERROR 6: No translation for Mercator_Auxiliary_Sphere to PROJ.4 format is known.

This is weird because the scripts works if the shapefile is not created by FME but by QGIS both in 900913.

this is the .prj of the not working shapefile

PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0],AUTHORITY["EPSG",3857]]

Does it means that OGR cannot transform file from FME?

Answer 1

For a complete history of the evolution from EPSG:900913 to EPSG 3785 and finally EPSG:3857, look at Análisis de Google Maps (in Spanish) and the conclusions are:

• The projection parameters defined for the EPSG: 900913 or 3785 describe a Mercator projection with geographic coordinates defined on a spherical model of the Earth, 6378.137m ratio, as specified in the following line (see EPSG:3785) :

  SPHEROID["Popular Visualisation Sphere", 6378137, 0, AUTHORITY["EPSG","7059"]]

• The official EPSG: 3857, in contrast, suggests a Mercator projection with coordinates defined on the WGS84 ellipsoid with a flattening (f) of 1/298.257223563:

  SPHEROID["WGS 84", 6378137, 298.257223563, AUTHORITY["EPSG","7030"]]

And the prj of the shapefile is EPSG 3857 and not EPSG 900913/3785 witch is (3785 ESRI .prj):

PROJCS["Popular Visualisation CRS / Mercator",GEOGCS["Popular Visualisation CRS",DATUM["D_Popular_Visualisation_Datum",SPHEROID["Popular_Visualisation_Sphere",6378137,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]

If I try with the definitions of the module osr, the transformation works:

from osgeo import osr
x,y = (6000,5000)
p3857 = osr.SpatialReference()
p3857.ImportFromEPSG(3857)
p900913 = osr.SpatialReference()
p900913.ImportFromEPSG(900913)
transformation = osr.CoordinateTransformation(p900913,p3857)
x2,y2 = transformation.TransformPoint(x, y)
print "%.3f, %.3f" % (x2,y2)
6000.000, 5000.000

and

p3857.ExportToWkt()   
'PROJCS["WGS_84_Pseudo_Mercator",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",
SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",   
0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["central_meridian",0],
PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],
PARAMETER["standard_parallel_1",0.0]]'

p900913.ExportToWkt()
'PROJCS["Google Maps Global Mercator",GEOGCS["WGS84",DATUM["WGS_1984",
SPHEROID["WG84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"],
UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_2SP"],
PARAMETER["standard_parallel_1",0],PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",0],,PARAMETER["false_easting",0],
PARAMETER["false_northing",0],UNIT["Meter",1],EXTENSION["PROJ4",
"+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 
+units=m  +nadgrids=@null +wktext  +no_defs"],AUTHORITY["EPSG","900913"]]'

But if you use the WKT provided by your .prj file, the result is:

p = osr.SpatialReference()
p.ImportFromWkt(yourwkt)
p.ExportToProj4()
ERROR 6: No translation for Mercator_Auxiliary_Sphere to PROJ.4 format is known.

So, it is the prj from FME that is the problem and not that of QGIS.

Answer 2

I think the problem here is simply one of competing definitions or at least names.

I'm told that WKT is fairly loose in some aspects (precision, order, names) and so you get vendor implementations that are technically the same coordinate system but with a slightly different definition.

In this case Esri's WKT for this coordinate system is slightly different to the EPSG:3857 standard. Not wrong, just different.

So the WKT string FME writes in the prj is valid (as valid as anything else) and is made to be recognized by ArcGIS. I guess it makes sense for FME to go with what Esri use in this scenario, since Shape is an Esri format!!!!

But apparently this prj is not recognized by OGR, which presumably sticks to the OpenGIS standard only. See here for another user with the same problem - who created their data as an Esri export rather than FME.

If you did want FME to write different content in your PRJ - that OGR would recognize - then you could simply add what we call an 'exception'. Look in the file <fme>/Reproject/Exceptions/esriwkt.db and you'll see how all these definitions are stored. You can add an exception to change the PRJ you get when writing Shape, or you can ask us to help you out.

I hope this helps. I find this article pretty helpful. But there is certainly a mountain of confusion over this one coordinate system (47 comments on that article) and if I've gotten anything wrong here it's because I'm not 100% certain myself!

I think on my deathbed, when they ask if I want a priest, I shall tell them to bring me a geodetics expert, as one final attempt to understand what it all means :-)