Revisiting Web Mercator

[warmerdam]

ESRI now uses the following definition for web mercator:

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["ESRI","102100"]]

NOTE: The Auxiliary Sphere Type parameter accepts 0 (use semimajor axis or radius of the geographic coordinate system), 1 (use semiminor axis or radius), 2 (calculate and use authalic radius), or 3 (use authalic radius and convert geodetic latitudes to authalic latitudes).

Erdas is now using the EPRJ_PSEUDO_MERCATOR:

"Pseudo Mercator" {
	Internal 68
	"Spheroid" <spheroid>
	"Longitude of central meridian" <4:angle ew:dd>
	"Latitude of true scale" <5:angle ns:dd>
	"False easting at central meridian" <6:distance ew:meters>
	"False northing at origin" <7:distance ns:meters>
}

EPSG now has a custom projection method, 1024, for google mercator (assuming major axis sphere I believe).

GDAL currently translates EPSG:3857 as:

PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    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],
    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","3857"],
    AXIS["X",EAST],
    AXIS["Y",NORTH]]

The EXTENSION element with a special proj.4 string that forces use of the semi-major sphere for projection calculations, but treated as WGS84.

I am thinking that we should move to treating the pseudo-mercator projection method better in WKT, and tying the various systems together properly. Either use a custom projection like ESRI ("Mercator_Auxiliary_Sphere") or introduce the "Auxilary_Sphere_Type" parameter to standard Mercator.

[mrosen]

I think the name of the projection method we want is "Popular Visualisation Pseudo Mercator" (no hyphen) ... that is what the EPSG Registry for 3857 ​http://www.epsg-registry.org/report.htm?type=selection&entity=urn:ogc:def:crs:EPSG::3857&reportDetail=short&style=urn:uuid:report-style:default-with-code&style_name=OGP%20Default%20With%20Code&title=epsg%203857 says.

I believe this issue has affected customers using GDAL to reproject files to 3857 ... the resulting WKT is misinterpretted by ArcGIS causing displacement errors.

[mrosen]

Apparently Spatialite users are also affected: ​http://lists.osgeo.org/pipermail/gdal-dev/2012-December/034978.html

[mrosen]

enhanced support for Web Mercator

[mrosen]

I am attaching a patch which adds the custom projection, Mercator_Auxiliary_Sphere as suggested above. The overall affect is to change that way that GDAL exposes the Web Mercator SRS. The old conversion is shown above. The patched conversion is quite similar to the ESRI definition (also above) and looks like this:

PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_Auxiliary_Sphere"],
    PARAMETER["central_meridian",0],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","3857"]]

Applications which depend on the WKT (e.g. MrSID) to expose the SRS will see improved recognition. I've verified the results in QGIS, ERDAS and ArcGIS.

The new projection, SRS_PT_MERCATOR_AUXSPHERE, has been integrated with the Proj4 converter and the OGR SRS validation mechanisms. Accordingly, I've removed the PROJ4 Extension 'override hack for google mercator.'

I have not attempted to update any format drivers to include support here. In particular, since it's unclear to me that there is any proper way to expose web mercator as a collection of geotiff tags, the GeoTiff conversion here is unchanged.

[warmerdam]

As dug up by Even, GeoToolkit apparently uses:

PROJCS["Google Mercator",
  GEOGCS["WGS 84",
    DATUM["World Geodetic System 1984",
      SPHEROID["WGS 84", 6378137, 298.257223563]],
    PRIMEM["Greenwich", 0],
    UNIT["degree", 0.017453292519943295]],
  PROJECTION["Mercator (1SP)"],
  PARAMETER["semi_major", 6378137],
  PARAMETER["semi_minor", 6378137],
  UNIT["m", 1]]

(from ​http://www.geotoolkit.org/modules/referencing/faq.html#Google) and geotools uses:

PROJCS["Google Mercator",
  GEOGCS["WGS 84",
    DATUM["World Geodetic System 1984",
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic latitude", NORTH],
    AXIS["Geodetic longitude", EAST],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Mercator_1SP"],
  PARAMETER["semi_minor", 6378137.0],
  PARAMETER["latitude_of_origin", 0.0],
  PARAMETER["central_meridian", 0.0],
  PARAMETER["scale_factor", 1.0],
  PARAMETER["false_easting", 0.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m", 1.0],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH],
  AUTHORITY["EPSG","900913"]]

from ​http://docs.geotools.org/latest/userguide/library/referencing/crs.html#google-maps

This leaves me still a uncertain what the best solution is.

[dfogel]

Melita Kennedy at ESRI offered some comments in an email thread where sharing b/w ArcGIS and GDAL was hindered. For what it's worth, this was with respect to both GeoTIFFs and Shapefiles.

Melita @ ESRI writes,

"GDAL is using the EPSG name, “WGS 84 / Pseudo-Mercator”, but not the operation name of “Popular Visualisation Pseudo Mercator”, which is what we would key on. In fact, at ArcGIS 10.2.1, we will match up “Popular Visualisation Pseudo Mercator” as being the same as our “Mercator Auxiliary Sphere” implementation.

This information is almost sufficient to guide a stable GDAL implementation. PROJ4 too. Superfluous parameters may be eliminated, or not.

[dfogel]

The EPSG Registry will show "Popular Visualisation Pseudo-Mercator" which will be parsed and identified as Web Mercator too ("[ESRI ignores] case, whitespace, and the characters “()/-_”."

[greenlaw]

ESRI apparently expects a semiminor axis of 6356752.314245179 as shown in the ArcGIS Desktop GUI:

WGS_1984_Web_Mercator_Auxiliary_Sphere
WKID: 3857 Authority: EPSG

Projection: Mercator_Auxiliary_Sphere
False_Easting: 0.0
False_Northing: 0.0
Central_Meridian: 0.0
Standard_Parallel_1: 0.0
Auxiliary_Sphere_Type: 0.0
Linear Unit: Meter (1.0)

Geographic Coordinate System: GCS_WGS_1984
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_WGS_1984
  Spheroid: WGS_1984
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314245179
    Inverse Flattening: 298.257223563

Thus the following proj4 string can be used to reproject to what ESRI expects:

+proj=merc +a=6378137 +b=6356752.314245179 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs

or using this WKT:

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["latitude_of_origin",0.0],
	PARAMETER["Standard_Parallel_1",0.0],
	PARAMETER["Auxiliary_Sphere_Type",0.0],
	UNIT["Meter",1.0],
	EXTENSION["PROJ4",
		"+proj=merc +a=6378137 +b=6356752.314245179 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
	],
	AUTHORITY["EPSG","3857"]
]

[mikrit]

I think the previous comment is a misunderstanding.

greenlaw wrote:

ESRI apparently expects a semiminor axis of 6356752.314245179 as shown in the ArcGIS Desktop GUI:

Well, but that semiminor axis is used only to specify the ellipsoid shape for the underlying geodetic datum (which is WGS84, so not spherical). However, the projection method Mercator_Auxiliary_Sphere will ignore the flattening and just use spherical formulas.

On the other hand, if you write +proj=merc in Proj.4 and use different semi-major axes specified by +a and +b, then Proj.4 will use ellipsoid formulas, which is not what you want in Web Mercator. So the EXTENSION suggested by greenlaw does not describe Web Mercator.

The following two examples, using first the original EXTENSION and then greenlaw's, shows that the results are different:

C:\>cs2cs +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 +to +proj=longlat +datum=WGS84
7000000 7000000
62d52'55.452"E  53d5'30.548"N 0.000

C:\>cs2cs +proj=merc +a=6378137 +b=6356752.314245179 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +to +proj=longlat +datum=WGS84
7000000 7000000
62d52'55.452"E  53d16'34.723"N 0.000

[wolf]

I'm working on a fix for this. Maybe anterim solution to recognize the Mercator_Auxilary_Sphere as EPSG:3857

[Even Rouault]

@wolf Yes in morphFromESRI(), if the projection method is Mercator_Auxilary_Sphere, resolving it by importFromEPSG(3857) is a reasonable approach (we want to make sure that there's a 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"] node in the WKT, what importFromEPSG(3857) will do)

[wolf]

Fixed in r27594

[greenlaw]

Screenshot showing misalignment of GDAL/OGR-produced shapefile with ESRI basemap

[greenlaw]

r27594 does not seem to fix the display offset problem when viewing a GDAL/ORG-produced EPSG:3857 GeoTIFF or Shapefile in an ESRI application, such as ArcMap.

Steps to reproduce:

1. Download U.S. States shapefile (EPSG:4269) from ​http://dds.cr.usgs.gov/pub/data/nationalatlas/statesp010g.shp_nt00938.tar.gz

2. Extract and reproject shapefile to EPSG:3857 using ogr2ogr:

   ogr2ogr -t_srs EPSG:3857 statesp010g_epsg3857.shp statesp010g.shp
   

3. Open ArcMap, set Data Frame projection to Web Mercator

4. Add the original (EPSG:4269) statesp010g.shp to map - it should align correctly using on-the-fly reprojection

5. Add statesp010g_epsg3857.shp to map - note the incorrect vertical alignment offset

Example of misalignment viewed over ESRI Web Mercator Topographic Basemap:

[wolf]

The problem is with the export part. I only fixed the importing of shapefiles with Mercator_Auxilary_Sphere to be recognized as Spherical mercator. Now the output of your sample command produces something that is WGS_84_Pseudo_Mercator with projection Mercator, which AFAIK ArcGIS doesn't recognize as spherical mercator.

Now, I could also implement the support for the Auxilary_Sphere_Type, and then would need to export a proper EPSG:3857 .prj file when producing that.

[greenlaw]

@wolf It would be great if exporting to ESRI-compatible EPSG:3857 format was supported in addition to import support. Has a separate ticket been opened for that?

Happy to help however I can.

[Even Rouault]

wolf, I believe there's a typo. Isn't there a missing I in auxilIary (SRS_PT_MERCATOR_AUXILARY_SPHERE) ?

[Jukka Rahkonen]

In the file /trunk/gdal/ogr/ogr_srs_esri.cpp, beginning from line 838 there is still this typo:

 else if( EQUAL(osProj, SRS_PT_MERCATOR_AUXILARY_SPHERE) )
    {
       // This is EPSG:3875 Pseudo Mercator. We might as well import it from
       // the EPSG spec.
       CPLString osAuxilarySphereType;
       importFromEPSG(3857);

[wolf]

Typo fixed in r28374

This still only fixes the input. Generating EPSG:3857 .prj is not yet done.

[bishop]

It seem to me that ogr2ogr produce wrong prj file:

  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]]

but should be

  PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    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","3857"]]

[Even Rouault]

Similar discussion dedicated to EPSG:3857 encoded in GeoTIFF in #5924

[bishop]

Shapes saved from 4326 to 3857 in different soft. Also images with shift demo

[Even Rouault]

This ticket has been automatically closed because Trac is no longer used for GDAL bug tracking, since the project has migrated to GitHub. If you believe this ticket is still valid, you may file it to ​https://github.com/OSGeo/gdal/issues if it is not already reported there.