Add 900913 EPSG code to cubewerks_extra file

[hobu]

PROJCS["Google Global 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.01745329251994328,
        AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_2SP"],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1],
    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"]]

[hobu]

After some discussion with Frank, we decided it should use WGS84 for the GEOGCS

PROJCS["Google Maps Global 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.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTIONMercator_2SP?,PARAMETER["standard_parallel_1",0],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],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"]]

[hobu]

Added in r12261

[gfleming]

I suspect this WKT is incorrect and inconsistent with other definitions. Surely it should be on a sphere, not on WGS84. The WKT and proj4 definitions in cubewerx_extra 900913 do not mean the same thing.

3 scenarios with gdal 1.7:

1. gdalwarp with t_srs "+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 +over +no_defs" writes this in the geotiff header:

   PROJCS["unnamed",
       GEOGCS["unnamed ellipse",
           DATUM["unknown",
               SPHEROID["unnamed",6378137,0]],
           PRIMEM["Greenwich",0],
           UNIT["degree",0.0174532925199433]],
       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"]]]
   

This is not very 'friendly' with all the 'unnamed' and 'unknown's but at least it is correct - it is on a sphere. The incomplete WKT means that GeoServer?, for example, does not recognise the SRS.

2. gdalwarp with t_srs "EPSG:900913" writes this in the geotiff header:

   PROJCS["Google Maps Global Mercator",
       GEOGCS["WGS 84",
           DATUM["WGS_1984",
               SPHEROID["WGS 84",6378137,298.257223563,
                   AUTHORITY["EPSG","7030"]],
               AUTHORITY["EPSG","6326"]],
           PRIMEM["Greenwich",0],
           UNIT["degree",0.0174532925199433],
           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"]]]
   

This is the intended result of this ticket, but the again, the definition is wrong as it uses GCS WGS84. It is still not the same as GeoServer?'s, for example so GeoServer? still does not recognise it as being 900913.

3. Geoserver maintains yet another WKT definition of WKT in epsg.properties. I forced gdalwarp to use this by adding it to the cubewerx_extra.wkt that gdal references, as follows, copying the proj4 EXTENSION from the 900913 line:

   900914,PROJCS["WGS84 / 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], AUTHORITY["EPSG","4326"]], PROJECTION["Mercator (1SP)", AUTHORITY["EPSG","9804"]], PARAMETER["semi_major", 6378137.0], 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],  AUTHORITY["EPSG","900913"],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"]]
   

Then I ran gdalwarp with t_srs "EPSG:900914" and it wrote this in the geotiff header (complete header):

Driver: GTiff/GeoTIFF
Files: two.tif
Size is 10669, 11712
Coordinate System is:
PROJCS["WGS84_Google_Mercator",
    GEOGCS["GCS_WGS_1984",
        DATUM["D_World Geodetic System 1984",
            SPHEROID["WGS_1984",6378137.0,298.257223563]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Mercator (1SP)"],
    PARAMETER["semi_major",6378137.0],
    PARAMETER["semi_minor",6378137.0],
    PARAMETER["latitude_of_origin",0.0],
    PARAMETER["central_meridian",0.0],
    PARAMETER["false_easting",0.0],
    PARAMETER["false_northing",0.0],
    UNIT["m",1.0],
    PARAMETER["standard_parallel_1",0.0]]
Origin = (3061095.039690456818789,-2998912.070841180160642)
Pixel Size = (0.557678323349409,-0.557678323349409)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 3061095.040,-2998912.071) 
Lower Left  ( 3061095.040,-3005443.599) 
Upper Right ( 3067044.910,-2998912.071) 
Lower Right ( 3067044.910,-3005443.599) 
Center      ( 3064069.975,-3002177.835) 
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  NoData Value=0
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha
  NoData Value=0

Result: I was hoping Geoserver would finally recognise this Coordinate System as 900913 since it matches its own definition but unfortunately it doesn't get to that point because I get this error when trying to save the store, probably because the rest of the header is now incomplete?:

"Could not list layers for this store, an error occurred retrieving them: Unable to acquire a reader for this coverage with format: GeoTIFF"

GeoServer?'s 900913 WKT definition also uses WGS84 but then goes on to override that by defining both semi-minor and semi-major axes as being the same (ie. a sphere). GeoServer?'s definition also purports to be from the latest EPSG database.

So, two things I think need to be done:

1. update the cubewerx_extra to a correct WKT definition (presumably the same EPSG one that geoserver uses)
   1. NB: must be on a sphere, not WGS84
2. check that gdal operations work with this new definition.

[warmerdam]

Gavin,

Good to hear from you!

As far as I'm concerned it is not a good idea to use the 900913 code. Instead use the preferred EPSG code for this coordinate system - EPSG:3857. This can be written to GeoTIFF properly though I don't know if it will be recognised by GeoServer? or not.

Note that "Mercator (1SP)" is a geoserver/geotools only way of naming the projection method, and in my opinion flaunts any effort at interoperability with other systems. It is not supported by GDAL.

Other than that, I'm not clear on what you saw as the primary issue. Was it that the WKT declares use of WGS84 as the earth model, but the PROJ.4 string declares use of a sphere? The Google Mercator coordinate system is intended to have the mercator projection operations done on the sphere, but for the lat/long coordinates to be treated as WGS84. This is difficult to represent in conventional GIS systems and in WKT. In order to get the external behavior of geographic coordinates being treated as WGS84 we declare that as the GEOGCS and then we use a custom PROJ.4 string to force the projection to be evaluated on the sphere (with the +nadgrids=@null magic telling PROJ.4 to treat the resulting lat/long as equivelent to WGS84).

I see no need for changes in the GDAL or PROJ.4 definitions.