Tilecache positions my gdal_translated map 20km south of where I want it

Question

I have created a geotiff file through gdal_translate:

gdal_translate -a_ullr 547812.2384341886 6834734.336439242 547967.5011478929 6834622.35081697 -a_srs EPSG:900913 pdz3.jpg pdz.tif

When I load this tif file in qgis, I see my map on the right location. Also, when I ask the gdalinfo from this file I get the following information:

Driver: GTiff/GeoTIFF
Files: pdz.tif
Size is 1024, 724
Coordinate System is:
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"]]]
Origin = (547812.238434188649990,6834734.336439241655171)
Pixel Size = (0.151623743851815,-0.154676273856209)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  547812.238, 6834734.336) (  4d55'15.89"E, 52d22'37.58"N)
Lower Left  (  547812.238, 6834622.351) (  4d55'15.89"E, 52d22'35.36"N)
Upper Right (  547967.501, 6834734.336) (  4d55'20.91"E, 52d22'37.58"N)
Lower Right (  547967.501, 6834622.351) (  4d55'20.91"E, 52d22'35.36"N)
Center      (  547889.870, 6834678.344) (  4d55'18.40"E, 52d22'36.47"N)
Band 1 Block=1024x2 Type=Byte, ColorInterp=Red
Band 2 Block=1024x2 Type=Byte, ColorInterp=Green
Band 3 Block=1024x2 Type=Byte, ColorInterp=Blue

I tested the 4d55'15.89"E, 52d22'37.58"N values in google map, and they also point to the proper location.

However, when I want this geotiff tiled - through tilecache - the map is placed about 20km south of the point I want it to appear (and appears in qgis). When I georeference the map in qgis self, and send it to tilecache I get the map on the right location. Calling gdalinfo on the qgis referenced tif I get roughly the same information (a couple of meters difference, due to mismatch clicks ;-))

The tilecache configuration file of the map is

[cache]
type=GoogleDisk
base=/home/gworx/data/tilecache

[jvv-chan-5-pdz]
type=GDAL
file=/home/gworx/data/layer-files/1320850804986-pdz.tif
bbox=547812.238434,6834734.33644,547967.501148,6834622.35082
spherical_mercator=true
srs=EPSG:900913
metaTile=true
extension=jpeg
mime_type=image/jpeg
tms_type=google

The bbox values are in line with the values I've set when calling gdal_translate. I am pretty much clueless at this moment on how to get this to work. If anyone can point me in the right direction on what I could be doing wrong, I would greatly appreciate it!

Answer 1

I suggest to use EPSG:3857 as CRS for the raster during the tiling process.

EPSG:900913 is not included in current GDAL and QGIS anymore, thus the special spherical mercator treatment is not applied. The proj definition for EPSG:3857 is

+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

which is a bit contradictional to

SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]

typically causing the offset of 20 to 30 km to the North.

Answer 2

20 km wrong in north/south direction is very probable to be caused by mixing spherical and ellipsoid representations of Earth.

Google Maps and Open Street Maps (OSM) use a tweaked Mercator (which seem to be called "Web Mercator") that is faster to calculate. The tile-software behind OSM likely also uses the spherical faster model, while your input data is probably using an ellipsoid earth representation.

The problem can possibly be fixed by adding a parameter semi_minor = 6378137 and re-projecting all data to a single common coordinate system.

For a long discussion, see (somewhat heavy to digest) http://alastaira.wordpress.com/2011/01/23/the-google-maps-bing-maps-spherical-mercator-projection/