Fig 1 – DIA Looks a Lot Better!
With a helpful comment from SharpGIS I was able to finally pin down my baffling problem with EPSG:54004.
The problem is in the datums.
ESRI:54004
PROJCS["World_Mercator",
GEOGCS["GCS_WGS_1984",
DATUM["WGS_1984",
SPHEROID["WGS_1984",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["Degree",0.017453292519943295]],
PROJECTION["Mercator_1SP"],
PARAMETER["False_Easting",0],
PARAMETER["False_Northing",0],
PARAMETER["Central_Meridian",0],
PARAMETER["Standard_Parallel_1",0],
UNIT["Meter",1],
AUTHORITY["EPSG","54004"]]
As Morten pointed out the 54004 datum includes a flattening, 298.257223563 :
SPHEROID["WGS_1984",6378137,298.257223563]],
So 54004 should be treated as an Ellipsoid rather than a Sphere.
There is a subtle difference in 900913. If you notice 900913 also includes a flattening:
SPHEROID["WGS_1984",6378137,298.257223563]],
EPSG:900913
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"]]
However, you might not notice in addition it includes an explicit minor axis parameter.
PARAMETER["semi_minor",6378137.0],
And this minor axis is identical to the major axis. The axis definition overrides the flattening in the Datum and is probably technically incorrect, but the idea was just to get a spherical mercator into a definition that people could use to match Google Maps. I’ve seen this definition in PostGIS, GeoServer, and OpenLayers.
I had already noticed this and played with a MercatorEllipsoid function to see if that would fix my problem. However, sadly, I made an error and miscalculated eccentricity. The real equation for e goes something like this:
double f = 1 / 298.257223563;
e = Math.Sqrt(2*f – Math.Pow(f,2));
resulting in e = 0.081819190842621486;
Once I made the correction for the proper eccentricity in MercatorEllipsoid, ESRI:54004 lines up with the EPSG:3857. DIA is back in its rightful place.
My MercatorEllipsoid function now calculates correct BBOX parameters for GetMap requests, but only for com.esri.wms.Esrimap services. Looks like ESRI is the expert and correctly produces ESRI:54004 with Datum as defined. However, not so with GeoServer.
GeoServer seems to ignore the flattening 298.257223563 or else assume it is like the 900913 where flattening is overridden by a minor axis parameter:
semi_minor axis = major.PARAMETER[semi_minor,6378137.0]
This leads to some problems. My WMS client now has to decide which service correctly interprets the DATUM on 54004. For now I just check for “com.esri.wms.Esrimap” in the WMS url and change datums accordingly. This will undoubtedly lead to problems with other services since I don’t yet know how MapServer or others treat 54004.
Summary
- ESRI is right again!
- Always check your math one more time
- Community responses produce answers
Thanks everyone!