GTiff Output Mercator Corner Coordinates Error using 1.10.1

[itbeyond]

I am using Mapserver to generate GTiff files using the following:

PROJECTION
 "init=epsg:3857"
END

OUTPUTFORMAT
 NAME "GTiff"
 DRIVER GDAL/GTiff
 MIMETYPE "image/tiff"
 IMAGEMODE RGB
 FORMATOPTION "COMPRESS=LZW"
END

When I create a map it works as expected however upon looking at the gdalinfo I notice the coordinates for the map corners are correct in the DMS but wrong in the Mercator as per this (using gdal 1.10.1 and Mapserver 6.4.0):

Driver: GTiff/GeoTIFF
Files: mapserv.tiff
Size is 1000, 1000
Coordinate System is:
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],
        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"]],
    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"]]
Origin = (12889228.989798404000000,-3724386.063095235300000)
Pixel Size = (25.399986284008403,-25.399986284008403)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=200
  TIFFTAG_YRESOLUTION=200
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (12889228.990,-3724386.063) (115d47' 9.29"E, 31d42'10.76"S)
Lower Left  (12889228.990,-3749786.049) (115d47' 9.29"E, 31d53'48.88"S)
Upper Right (12914628.976,-3724386.063) (116d 0'50.71"E, 31d42'10.76"S)
Lower Right (12914628.976,-3749786.049) (116d 0'50.71"E, 31d53'48.88"S)
Center      (12901928.983,-3737086.056) (115d54' 0.00"E, 31d48' 0.00"S)
Band 1 Block=1000x2 Type=Byte, ColorInterp=Red
Band 2 Block=1000x2 Type=Byte, ColorInterp=Green
Band 3 Block=1000x2 Type=Byte, ColorInterp=Blue

If you take coordinate Upper Left as an example: DMS states: 115d47' 9.29"E, 31d42'10.76"S Merc states: 12889228.990,-3724386.063

The actual coordinates of the map are right in DMS but the Merc does not match the DMS after a WGS84 - Merc conversion as below:

DMS 115d47' 9.29"E, 31d42'10.76"S to DD: 115.7859 -31.70299

gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857
115.7859 -31.70299
12889227.4290409 -3724386.28296906 0

DOES NOT EQUAL
12889228.990 -3724386.063

This is not the same/correct coordinate pair and all the Merc corner coordinates in the gdalinfo are incorrect.

I believe mapserver is doing it right because if I project the map to epsg:4283 the corner coordinates are totally right, so I suspect mapserver is outputing it right.

​http://www.itbeyond.com.au/images/mapserv.tiff for reference.

[Even Rouault]

There's no problem in fact. It just depends on the number of decimals used. If you use 12889228.989798404000000,-3724386.063095235300000 ("Origin") instead of 12889228.990,-3724386.063, then you have :

$ gdaltransform -s_srs EPSG:3857 -t_srs EPSG:4326 12889228.989798404276371 -3724386.063095235265791 115.785914020523 -31.7029883195659 0

and 115.785914020523 = 115° 47' 9.2904738828

What is displayed in "Corner Coordinates" is meant as being for information purposes, so millimeters accuracy for coordinates in meters, or 1/100 of second for DMS, are most of the time sufficient. If you need more precision, you have to use the information of the geotransform ("Origin" and "Pixel Size")

And in your own conversion from DMS to decimal degrees, you have lost precision too. If you use all precision from DMS as displayed in Coorner Coordinates (which, as said above, is not full precision), you'll get a closer result to the mercator coordinates :

$ gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857
115.78591388889  -31.702988888888
12889228.9751451 -3724386.13758731 0