Opened 9 years ago

Closed 7 months ago

#3962 closed defect (wontfix)

Revisiting Web Mercator

Reported by: warmerdam Owned by: wolf
Priority: normal Milestone: closed_because_of_github_migration
Component: OGR_SRS Version: unspecified
Severity: normal Keywords: web mercator.
Cc: mrosen

Description

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.

Attachments (3)

gdal_19_ticket3962.patch (12.2 KB) - added by mrosen 6 years ago.
enhanced support for Web Mercator
3857_esri_misalignment.jpg (143.4 KB) - added by greenlaw 5 years ago.
Screenshot showing misalignment of GDAL/OGR-produced shapefile with ESRI basemap
3857.zip (236.4 KB) - added by bishop 5 years ago.
Shapes saved from 4326 to 3857 in different soft. Also images with shift demo

Download all attachments as: .zip

Change History (23)

comment:1 Changed 7 years ago by 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.

comment:2 Changed 7 years ago by mrosen

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

Changed 6 years ago by mrosen

Attachment: gdal_19_ticket3962.patch added

enhanced support for Web Mercator

comment:3 Changed 6 years ago by 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.

comment:4 Changed 6 years ago by warmerdam

Cc: mrosen added

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.

comment:5 Changed 6 years ago by 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.

comment:6 Changed 6 years ago by 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 “()/-_”."

comment:7 Changed 6 years ago by 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"]
]

comment:8 Changed 5 years ago by 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

comment:9 Changed 5 years ago by wolf

Owner: changed from warmerdam to wolf

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

comment:10 Changed 5 years ago by 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)

comment:11 Changed 5 years ago by wolf

Resolution: fixed
Status: newclosed

Fixed in r27594

Changed 5 years ago by greenlaw

Attachment: 3857_esri_misalignment.jpg added

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

comment:12 Changed 5 years ago by 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
  1. Extract and reproject shapefile to EPSG:3857 using ogr2ogr:
    ogr2ogr -t_srs EPSG:3857 statesp010g_epsg3857.shp statesp010g.shp
    
  2. Open ArcMap?, set Data Frame projection to Web Mercator
  1. Add the original (EPSG:4269) statesp010g.shp to map - it should align correctly using on-the-fly reprojection
  1. Add statesp010g_epsg3857.shp to map - note the incorrect vertical alignment offset

Example of misalignment viewed over ESRI Web Mercator Topographic Basemap:

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

comment:13 Changed 5 years ago by 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.

comment:14 Changed 5 years ago by 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.

comment:15 Changed 5 years ago by Even Rouault

Resolution: fixed
Status: closedreopened

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

comment:16 Changed 5 years ago by 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);

comment:17 Changed 5 years ago by wolf

Typo fixed in r28374

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

comment:18 Changed 5 years ago by 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"]]

comment:19 Changed 5 years ago by Even Rouault

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

Changed 5 years ago by bishop

Attachment: 3857.zip added

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

comment:20 Changed 7 months ago by Even Rouault

Milestone: closed_because_of_github_migration
Resolution: wontfix
Status: reopenedclosed

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.

Note: See TracTickets for help on using tickets.