GIS in XML

Fig 1 – Something is not right with DIA?

Life is full of problems. With eight kids we have lots of problems, especially as school gets started up again. Life in the chaos lane has its moments, but the problem that has been baffling me this week has to do with EPSG coordinate systems and in particular the venerable but unofficial EPSG:54004.

EPSG:54004 is really not official. You won’t find it in the official EPSG registry. It is commonly used, however, ever since it was added as a “user extension” by Google, originally for use in Google maps.

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

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

And both are similar to the official, official EPSG:3857

They should be equivalent. All are using earth radius 6378137.0m with a spherical mercator approximation that happens to be very efficient for quad tree tile pyramids. It has been so successful that it’s used in all the online map services, Google, Bing, Yahoo, and Open Street Map. It is also pretty simple to convert latitude, longitude into spherical Mercator coordinates. Using my dog eared John Snyder’s Manual I get something like this:

private Point Mercator(double lon, double lat)
{
        /* spherical mercator for Google, Bing, Yahoo
        * epsg:900913, epsg:54004, epsg:3857 R= 6378137
        * x = longitude
        * y= R*ln(tan(pi/4 +latitude/2)
        */············
        double R = 6378137.0;
        double x = R * Math.PI / 180 * lon;
        double y = R * Math.Log(Math.Tan(Math.PI / 180 *
         (45 + lat / 2.0)));
        return new Point(x, y);
}


The Problem

Now on to the problem, apparently not all EPSG:54004s are equal?

Looking at the screen shot above you can see that the MRLC NLCD 2001 impervious Surface layer using EPSG:54004 shows DIA a few miles south of the Bing Maps DIA? Bing Maps, as noted above, uses the EPSG:3857 spherical mercator which should be the same, but apparently is not? A quick check against OSM and Yahoo base maps verifies the same discrepancy.

I’ve been known to get tangled up in coordinate systems, so just as a sanity check here is a screenshot with DRCOG’s Network layer also rendered as EPSG:54004:

Fig 2 – DRCOG EPSG:54004 looks right with Bing Maps but wrong with MRLC?

The DRCOG EPSG:54004 and Bing Maps align correctly so it is now 2 to 1 against the MRLC EPSG:54004. But what is going on here?

Taking a peek at DRCOG’s GetCapabilities I can see it is served out of GeoServer. “http://gis.drcog.org:80/geoserver/wms” Here is the user epsg section of GeoServer’s epsg.properties file from my local copy of GeoServer:

54004=PROJCS["WGS84 / Simple Mercator", GEOGCS["WGS 84", DATUM["WGS_1984",
SPHEROID["WGS_1984", 6378137.0, 298.257223563]],
PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295]],
PROJECTION["Mercator_1SP_Google"],
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["x", EAST], AXIS["y", NORTH], AUTHORITY["EPSG","54004"]]

It looks familiar.

The MRLC GetCapabilites shows an ESRI server,
“http://ims.cr.usgs.gov:80/wmsconnector/com.esri.wms.Esrimap.”
Perhaps this is a clue. I know of several other ESRI WMS services so going to the National Atlas, I can see, yes, it publishes an EPSG:54004 SRS: <SRS>EPSG:54004</SRS>

I plug the National Atlas into my handy WMSClient and take a look. I selected the major highways and interstates which are visible, with opacity turned up, against the Impervious Surface backdrop. You can see a high degree of correlation between the two layers, and both are no where near the Bing base map?

Fig 3 – National Atlas EPSG:54004 agrees with impervious surface?

OK, now it’s 2 to 2 and the end of the ninth. I am clear that the EPSG:54004 discrepancy has something to do with a difference in interpretation between WMS servers. It appears that ESRI is consistently different from Bing Maps and GeoServer. ESRI is the expert in this area, so my thought is there may be some unusual EPSG consideration that ESRI has caught and others possibly have not.

A bit of Googling brings up an older ESRI technical article that may have some bearing. However, nothing very clear about this large of a discrepancy. Going with a hunch I tried pretending ESRI EPSG:54004 was ellipsoidal rather than spherically based. Doing a quick switch to an Ellipsoidal Mercator calculation, however, makes no discernable difference.

private Point MercatorEllipsoid(double lon, double lat)
{
        double a = 6378137.0;// meters equitorial radius            
        double e = 1 / 298.257223563;
        double ctrmerid = 0.0;

        double px = lon * Math.PI / 180;
        double py = lat * Math.PI / 180;
        double pc = ctrmerid * Math.PI / 180;

        double x = a * (px - pc);
        double y = a * (Math.Log(Math.Tan((45.0 * Math.PI / 180) + py / 2.0)
        * Math.Pow((1.0 - e * Math.Sin(py))
        / (1.0 + e * Math.Sin(py)), e / 2.0)));
        return new Point(x, y);

}

Switching to the universal EPSG:4326 lined things up well enough. Unfortunately it does nothing to solve the mystery.

Fig 4 – Switching to EPSG:4326 brings layers into alignment at this zoom level

Summary

Mystery is not resolved. If others have some insight, I’d like to hear it, especially from experts with ESRI com.esri.wms.Esrimap experience. If anyone wants to play with these EPSG:54004 layers I have a simple WMS client here: Silverlight WMSClient.

This entry was posted on Friday, August 28th, 2009 at 6:21 pm and is filed under ArcMap, GIS, Silverlight. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.