EPSG Problems


EPSG:54004 problem fig 1
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.

Here are the basics of EPSG: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"]]
It is similar to 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"]]

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:

EPSG:54004 problem fig 2
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?

EPSG:54004 problem fig 3
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.

EPSG:4326  fig 3
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.

Neo vs Paleo Geography

Jurassic Dinosaur skeleton
Fig 1 – Paleo, Neo – what about Jurassic Geography?

I gather that there is some twittering about neo versus paleo geography. See Peter Batty’s blog entry or James Fee’s blog. I myself don’t Twitter, but in general I’m happy for Peter’s paleo accomodation of the non twitterers, repeating the conversation in a blog entry. Peter has also updated comments with a new post questioning, “Are we now in a post neogeography era?” The dreaded paradigm shifts are coming fast and furiously.

I am not really able to comment on neo vs paleo as I myself fall even further back into “Jurassic Geography.” Looking at connotations we have this accounting:

····neo - 1990 – present, new, recent, different, Obama, Keynesian, Apple, Google Earth, Cloud, Java C# RubyRails Twitter

····paleo - as in paleolithic 2.8m – 10k old, prehistoric, ancient, early, primitive, Nixon, Supply Side, Microsoft, Windows Desktop, ESRI Arc???, C C++ Javascript telephone

Obviously the “paleo” label is not carried with quite the honor of “neo.” It’s reminiscent of the Galen / Myers-Brigg personality typology characterized as Lion, Otter, Beaver, and Golden Retriever. What do you want to be? Obviously not the beaver, but there has to be a significant part of the world in that category, like it or not. After all what would lions eat for dinner without a few of us beavers? Likewise there is a significant branch of paleo in the GIS kingdom.

However, in the pre-paleolithic era there are still a few of us left, falling into the “long tail” of the Jurassic. So carrying on down the connotation stream here is the Jurassic geography equivalent:

····jurassic – 206m-144m dinosaurs, fossils, pre paleolithic, Hoover, laissez faire, IBM Big Iron, Assembly Cobol, open source

Wait “Open Source” – Jurassic Geography? How did that get in there? The notoriously frugal days of Hoover never made it into the paleolithic era’s “Supply Side” economy. It’s Keynesian economics all over the neo world, so Jurassic geography is the frugal end of the spectrum and how can you get more frugal than free! Obviously Open Source is as Jurassic as they come in Geography circles.

As I’ve recently been in a gig hunting mode, I’ve been having quite a few in depth conversations about GIS stacks. As a small businessman outside the corporate halls of paleo geography, I’ve had few occasions to get an in depth education on the corporate pricing world. So I spent the past couple of days looking into it.

Let’s start at the bottom of the stack. Here is some retail pricing on a few popular GIS databases:

  • Oracle Standard Spatial $17,500 + $3850 annual
  • Oracle Enterprise Locator $47,500 + $10,450 annual
  • SQL Server 2008 Web edition ~ $3500
  • PostgreSQL/PostGIS $0.00

If you’re a Jurassic geographer which do you choose? Probably not Oracle Enterprise Locator. If your Paleo you look at that and think, “Man, I am the negotiator! We don’t pay any of that retail stuff for the masses.” Neo? – well how would I know how a neo thinks?

Next take a look at the middle tier:

  • ESRI ArcGIS Server standard workgroup license
    ····Minimum $5000 2cores + $1250 2core annual
    ····Additional cores $2500/core + $625/core annual
  • ESRI ArcGIS hosted application server license
    ····Minimum $40,000 4 cores + $10,000 4 core annual
    ····Additional cores $10,000/core + $2500/core annual
  • OWS GeoServer or MapServer minimum $0 + annual $0
    But, this is GIS isn’t it? We want some real analytic tools not just a few hundred spatial functions in JTS Topology suite. OK, better throw in a few QGIS or GRASS installs and add a few $0s to the desktop production. Oh, and cores, we need some, “make that a 16core 64 bit please” – plus $0.

I think you catch the Jurassic drift here. How about client side.

  • ESRI Silverlight free, well sort of , if you’re a developer, NGO, educational, or non-profit otherwise take a look at that ArcGIS license back a few lines.
  • Google API it’s Neo isn’t it? $10k per annum for a commercial use, maybe its Paleo after all.
  • Virtual / Bing Maps api $8k per annum transaction based and in typical license obfuscation fashion impossible to predict what the final cost will be. Paleo, “Just send me the invoice.”
  • OpenLayers is a javascript api client layer too, just solidly Jurassic at $0
  • Silverlight well it can be Jurassic, try DeepEarth over at codeplex or MapControl from Microsoft with the Bing imageservice turned off, OSM on.

It’s been an interesting education. Here is the ideal Jurassic GIS stack:
Amazon EC2 Windows instance + PostGIS database + GeoServer OWS + IIS Silverlight MapControl client
The cost: starts at $100/mo(1 processor 1.7Gb 32bit) up to $800/mo(4 processor 15Gb 64bit)

So what does a Jurassic geographer get in this stack?

Hardware:
Amazon Cloud based virtual server, S3 Backup, AMI image replication, Elastic IP, AWS console, choice of OS, cores, memory, and drive space. Ability to scale in a big way with Elastic load balancing, auto scaling, and CloudWatch monitoring. Performance options like CloudFront edge service or something experimental like Elastic MapReduce Hadoop clusters.

Database:
PostgreSQL/PostGIS – Standards compliant SQL server with GIST spatial indexing on OGC “Simple Features for SQL” specification compliant geometry with extended support for 3DZ, 3DM and 4D coordinates. A full set of roughly 175 geometry, management, and spatial functions. It supports almost all projections. All this and performance? maybe a little vague but not shabby:

“PostGIS users have compared performance with proprietary databases on massive spatial data sets and PostGIS comes out on top.”

Middle Tier:
Geoserver – standards compliant OWS service for WMS, WFS, WCS.
Data sources: Shapefile, Shapefile Directory, PostGIS, external WFS, ArcSDE, GML, MySQL, Oracle, Oracle NG, SQL Server, VPF
Export formats: WFS GML, KML, SVG, PDF, GeoRSS, Png, Jpeg, Geotiff, OGR Output – MapInfo Tab and MID/MIF, Shp, CSV, GeoJSON …
OGC standard SLD styling, built in gwc tile caching – seeded or as needed, managed connection pools, RESTful configuration api, and ACEGI integrated security.

WCS adds :

  1. ArcGrid – Arc Grid Coverage Format
  2. ImageMosaic – Image mosaicking plugin
  3. WorldImage – A raster file accompanied by a spatial data file
  4. Gtopo30 – Gtopo30 Coverage Format
  5. GeoTIFF – Tagged Image File Format with Geographic information
“GeoServer is the reference implementation of the Open Geospatial Consortium (OGC) Web Feature Service (WFS) and Web Coverage Service (WCS) standards, as well as a high performance certified compliant Web Map Service (WMS). “

Browser client viewer:
Take your pick here’s a few:

Summary:
Well in these economic times Jurassic may in fact meet Neo. The GIS world isn’t flat and Jurassic going right eventually meets Neo going left, sorry Paleos. Will Obama economics turn out to be Hooverian in the end? Who knows, but here’s a proposition for the Paleos:

Let me do a GIS distribution audit. If I can make a Jurassic GIS Stack do what the Paleo stack is currently providing, you get to keep those annual Paleo fees from here to recovery. How about it?

rkgeorge @cadmaps.com

mxd to SLD – ArcMap2SLD


ArcMap2SLD screenshot
Fig 1 – ArcMap2SLD tool running with ArcMap 9.3

I had some more time to work out the ArcMap2SLD script tool.

Using an mxd in ArcMap 9.3, I first made sure all of my shape layers were connected. Then I ran the ArcMap2SLD tool which requires a session of ArcMap to be running.

This is a simple ArcMap project with a set of shp layers. It takes about 5-10 minutes to run through all the layers. The analysis seems to read through each symbol in a layer and then to loop through every record. After a few minutes it asks again for a file to store the sld and now the test.sld is ready.

The resulting SLD is a single file with all of the Feature types and rules listed together:

<?xml version="1.0" encoding="ISO-8859-1" standalone="yes"?>
<sld:StyledLayerDescriptor version="1.0.0"
xmlns:sld="http://www.opengis.net/sld"
xmlns:ogc="http://www.opengis.net/ogc"
xmlns:xlink="http://www.w3.org/1999/xlink">
  <sld:NamedLayer>
    <sld:Name>ovroads</sld:Name>
    <sld:UserStyle>
      <sld:Name>Style1</sld:Name>
      <sld:FeatureTypeStyle>
        <sld:FeatureTypeName>ovroads</sld:FeatureTypeName>
        <sld:Rule>
          <sld:Name>Interstates</sld:Name>
          <sld:Title>Interstates</sld:Title>
          <ogc:Filter>
            <ogc:PropertyIsEqualTo>
              <ogc:PropertyName>RDCLS</ogc:PropertyName>
              <ogc:Literal>010</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <sld:LineSymbolizer>
            <sld:Stroke>
              <sld:CssParameter name="stroke">#000000</sld:CssParameter>
              <sld:CssParameter name="stroke-width">0.8</sld:CssParameter>
              <sld:CssParameter name="stroke-opacity">1</sld:CssParameter>
            </sld:Stroke>
          </sld:LineSymbolizer>
          <sld:LineSymbolizer>
            <sld:Stroke>
              <sld:CssParameter name="stroke">#FFFB86</sld:CssParameter>
              <sld:CssParameter name="stroke-width">2.6</sld:CssParameter>
              <sld:CssParameter name="stroke-opacity">1</sld:CssParameter>
            </sld:Stroke>
          </sld:LineSymbolizer>
          <sld:LineSymbolizer>
            <sld:Stroke>
              <sld:CssParameter name="stroke">#000000</sld:CssParameter>
              <sld:CssParameter name="stroke-width">3.4</sld:CssParameter>
              <sld:CssParameter name="stroke-opacity">1</sld:CssParameter>
            </sld:Stroke>
          </sld:LineSymbolizer>
        </sld:Rule>
            .
            .

Copying a subset for the roads layer into GeoServer style did not produce any features the first time around.

First the Filter PropertyName is case sensitive. The loading batch I used produces lower case field names in the PostGIS tables. As a result the RDCLS property name had to be changed to rdcls through out the SLD.

If you notice in the above SLD output, the Geometry reference is missing:
    <Geometry><ogc:PropertyName>the_geom</ogc:PropertyName></Geometry>

Also missing are scale elements:
    <MinScaleDenominator>32000</MinScaleDenominator>
    <MaxScaleDenominator>20000000</MaxScaleDenominator>

After the necessary modification the SLD rules look like this:

        <Rule>
          <Name>U S Highways</Name>
          <Title>U S Highways</Title>
          <ogc:Filter>
            <ogc:PropertyIsEqualTo>
              <ogc:PropertyName>rdcls</ogc:PropertyName>
              <ogc:Literal>014</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <MinScaleDenominator>32000</MinScaleDenominator>
          <MaxScaleDenominator>20000000</MaxScaleDenominator>
          <LineSymbolizer>
            <Geometry>
              <ogc:PropertyName>the_geom</ogc:PropertyName>
            </Geometry>
            <Stroke>
              <CssParameter name="stroke">#FA3411</CssParameter>
              <CssParameter name="stroke-width">3.4</CssParameter>
              <CssParameter name="stroke-opacity">1</CssParameter>
            </Stroke>
          </LineSymbolizer>
        </Rule>

This still requires some manual intervention, but it does handle setting Filters, stroke color, stroke-width, and stroke-opacity parameters with a text editor along with cut & paste. It handled the PolygonSymbolizer in the mxd as well.

This is a time saver, however, the appearance didn’t seem to match in the resulting SLD tilesources. It appears as though some of the SLD stroke widths end up masking underlying strokes so that for example the purple highway symbolization does not appear in the SLD tilesource version.

ArcMap2SLD screenshot
Fig 2 – Silverlight tile source with new SLD

After some examination I realized that multiple LineSymbolizations were in reverse order. The thickest line width needs to be first with thinner stroke widths later in the SLD file, because rendering is in file order. Here is an example of the final resulting SLD LineSymbolization with corrected ordering:

        <Rule>
          <Name>Boulder Turnpike</Name>
          <Title>Boulder Turnpike</Title>
          <ogc:Filter>
            <ogc:PropertyIsEqualTo>
              <ogc:PropertyName>rdcls</ogc:PropertyName>
              <ogc:Literal>016</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <MinScaleDenominator>32000</MinScaleDenominator>
          <MaxScaleDenominator>20000000</MaxScaleDenominator>
          <LineSymbolizer>
            <Geometry><ogc:PropertyName>the_geom</ogc:PropertyName></Geometry>
            <Stroke>
              <CssParameter name="stroke">#000000</CssParameter>
              <CssParameter name="stroke-width">3.4</CssParameter>
              <CssParameter name="stroke-opacity">1</CssParameter>
            </Stroke>
          </LineSymbolizer>

          <LineSymbolizer>
            <Geometry><ogc:PropertyName>the_geom</ogc:PropertyName></Geometry>
            <Stroke>
              <CssParameter name="stroke">#FFFB86</CssParameter>
              <CssParameter name="stroke-width">2.6</CssParameter>
              <CssParameter name="stroke-opacity">1</CssParameter>
            </Stroke>
          </LineSymbolizer>

          <LineSymbolizer>
            <Geometry><ogc:PropertyName>the_geom</ogc:PropertyName></Geometry>
            <Stroke>
              <CssParameter name="stroke">#000000</CssParameter>
              <CssParameter name="stroke-width">0.8</CssParameter>
              <CssParameter name="stroke-opacity">1</CssParameter>
            </Stroke>
          </LineSymbolizer>
        </Rule>

ArcMap2SLD screenshot
Fig 3 – Silverlight tile source with new SLD with corrected order

Summary

This tool appears to be a helpful start for translating mxd style to sld. The result gives colors and basic widths etc but requires several modifications:

  1. match case to the PostGIS table
  2. add Geometry elements
  3. add min and max scale denominators
  4. reverse order of multiple LineSymbolizations

It’s still a good start.