1

I am trying to convert from EPSG:3857 to ESRI:54009 and facing some issues (pyproj version -2.6.1).

ERROR 1: PROJ: proj_create_from_database: crs not found

ERROR 1: PROJ: proj_create: unrecognized format / unknown name

ERROR 6: Cannot find coordinate operations from EPSG:3857' to `'``

My code as following:

def project_geometry(geom , in_crs = 3857 ,out_crs = 54009):
    inSpatialRef = osr.SpatialReference()   
    inSpatialRef.ImportFromEPSG(in_crs)
    outSpatialRef = osr.SpatialReference()
    outSpatialRef.ImportFromEPSG(out_crs)
    coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
    geom =ogr.CreateGeometryFromWkb(geom.to_wkb())
    geom.Transform(coordTrans)
    return loads(geom.ExportToWkb())

Is there a direct way to do this transformation?

or do I need to use transformation from EPSG:3857 to other CRS and than from that to ESRI:54009?

Also I didn't find transformation files in this pyproj release.

4

Your output CRS is not defined by EPSG but by ESRI. EPSG:54009 does not exist. For ESRI:54009 use ImportFromESRI https://gdal.org/python/osgeo.osr.SpatialReference-class.html#ImportFromESRI.

EDIT

User @snowman showed in the comments the way to get SpatialReference by the ESRI code:

>>> from osgeo import osr
>>> outSpatialRef = osr.SpatialReference()
>>> outSpatialRef.SetFromUserInput('ESRI:54009')
0
>>> print outSpatialRef
PROJCS["World_Mollweide",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Mollweide"],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["ESRI","54009"]]

Proj does know ESRI:54009. You can check with the projinfo utility: https://proj.org/apps/projinfo.html.

projinfo esri:54009
PROJ.4 string:
+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs

Proj also knows what operation to use for transforming data between EPSG:3857 and ESRI:54009. That can be verified also with projinfo

projinfo -s epsg:3857 -t esri:54009
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of Popular Visualisation Pseudo-Mercator + axis order change (2D) + World_Mollweide, 0 m, World

PROJ string:
+proj=pipeline +step +inv +proj=webmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +step +proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84

WKT2:2019 string:
CONCATENATEDOPERATION["Inverse of Popular Visualisation Pseudo-Mercator + axis order change (2D) + World_Mollweide",
    SOURCECRS[
        PROJCRS["WGS 84 / Pseudo-Mercator",
            BASEGEOGCRS["WGS 84",
                DATUM["World Geodetic System 1984",
                    ELLIPSOID["WGS 84",6378137,298.257223563,
                        LENGTHUNIT["metre",1]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433]],
                ID["EPSG",4326]],
            CONVERSION["Popular Visualisation Pseudo-Mercator",
                METHOD["Popular Visualisation Pseudo Mercator",
                    ID["EPSG",1024]],
                PARAMETER["Latitude of natural origin",0,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8801]],
                PARAMETER["Longitude of natural origin",0,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8802]],
                PARAMETER["False easting",0,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8806]],
                PARAMETER["False northing",0,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8807]]],
            CS[Cartesian,2],
                AXIS["easting (X)",east,
                    ORDER[1],
                    LENGTHUNIT["metre",1]],
                AXIS["northing (Y)",north,
                    ORDER[2],
                    LENGTHUNIT["metre",1]],
            ID["EPSG",3857]]],
    TARGETCRS[
        PROJCRS["World_Mollweide",
            BASEGEOGCRS["WGS 84",
                DATUM["World Geodetic System 1984",
                    ELLIPSOID["WGS 84",6378137,298.257223563,
                        LENGTHUNIT["metre",1]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["Degree",0.0174532925199433]]],
            CONVERSION["World_Mollweide",
                METHOD["Mollweide"],
                PARAMETER["Longitude of natural origin",0,
                    ANGLEUNIT["Degree",0.0174532925199433],
                    ID["EPSG",8802]],
                PARAMETER["False easting",0,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8806]],
                PARAMETER["False northing",0,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8807]]],
            CS[Cartesian,2],
                AXIS["(E)",east,
                    ORDER[1],
                    LENGTHUNIT["metre",1]],
                AXIS["(N)",north,
                    ORDER[2],
                    LENGTHUNIT["metre",1]],
            ID["ESRI",54009]]],
    STEP[
        CONVERSION["Inverse of Popular Visualisation Pseudo-Mercator",
            METHOD["Inverse of Popular Visualisation Pseudo Mercator",
                ID["INVERSE(EPSG)",1024]],
            PARAMETER["Latitude of natural origin",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8802]],
            PARAMETER["False easting",0,
                LENGTHUNIT["metre",1],
                ID["EPSG",8806]],
            PARAMETER["False northing",0,
                LENGTHUNIT["metre",1],
                ID["EPSG",8807]],
            ID["INVERSE(EPSG)",3856]]],
    STEP[
        CONVERSION["axis order change (2D)",
            METHOD["Axis Order Reversal (2D)",
                ID["EPSG",9843]],
            ID["EPSG",15498]]],
    STEP[
        CONVERSION["World_Mollweide",
            METHOD["Mollweide"],
            PARAMETER["Longitude of natural origin",0,
                ANGLEUNIT["Degree",0.0174532925199433],
                ID["EPSG",8802]],
            PARAMETER["False easting",0,
                LENGTHUNIT["metre",1],
                ID["EPSG",8806]],
            PARAMETER["False northing",0,
                LENGTHUNIT["metre",1],
                ID["EPSG",8807]]]],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]]]
| improve this answer | |
  • Thanks for replay @user30184 ! when I try to use outSpatialRef.ImportFromESRI(54009) I get an error ` return _osr.SpatialReference_ImportFromESRI(self, *args) TypeError: not a sequence` – Toren Jun 22 at 16:57
  • 1
    By looking at the examples at github.com/OSGeo/gdal/blob/master/autotest/osr/osr_esri.py ImportFromESRI seems to require some sort of WKT as input. It may be easiest to use the reported Proj string +proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs for defining the out_crs. – user30184 Jun 22 at 19:33
  • 2
    I would recommend SetFromUserInput('ESRI:54009') – snowman2 Jun 23 at 1:11
  • Thanks @snowman, I hope I did it right in my edited answer. – user30184 Jun 23 at 7:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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