How to georeference a web mercator tile correctly using gdal?

Question

As an example I will take the following tile http://a.tile.openstreetmap.org/3/4/2.png and save it as "4_2.png".

The WGS84 coordinates of this tile can be calculated or read there by clicking the corresponding tile:

0 66.51326044311185 45 40.97989806962013 (West North East South)

How to georeference the tile correctly (using gdal to generate a geotiff or another georeferenced format) so that:

• The bitmap don't need to be stretched (= the pixels in the geotiff are exactly the same as in the original bitmap)
• The resulting image will be opened at the right place in a GIS viewer/editor (like for example in TatukGIS Free Viewer)?

(Edited at Sept 19 2011 to make my question clearer and include my conclusions)

---

My conclusion:

I first though that the third idea (see below) was the right one. I opened the geotiff in a GIS Viewer and compared the displayed coordinates with what I expected. The geotiff out of the second idea seems to be shifted 2 pixels towards north. That is why I considered idea 3 (or 4) as the right one.
But if you try with a tile at a much higher zoom level, the geotiff out of idea 3 is definitively shifted towards south. It was silly to compare coordinates on a tile of zoom level 3. The country boundaries at such zoom level are mus simplified so that the comparison doesn't give good results.

Dan S. was right, the tile image is already in EPSG:3857. The second idea is then the right one (and gives good result at high zoom levels too)

---

First idea: EPSG:4326
The EPSG code for WGS84 coordinates is EPSG:4326. So I simply use the WGS84 coordinates to georefence the tile as geotif using gdal_translate:

gdal_translate -of Gtiff -co tfw=yes -a_ullr 0 66.51326044311185 45 40.97989806962013 -a_srs EPSG:4326 4_2.png t4326.tif

The resulting map is displayed at the right place but I fear that the projection is not correct and that there might be a shift in the middle of the tile. After trying a long time to check that by reprojecting the map with gdalwarp, I downloaded a demo version of Global Mapper and this seems to be the case (sames borders as by idea 3 but a shift inside the tile). The image should be streched to be able to use EPSG:4326 coordinates .

---

Second idea: EPSG:3857
This tile use a "web mercator" projection (alias google map projection), which now has an EPSG code: EPSG:3857 (alias EPSG:900913). I simply convert the coordinates using gdaltransform:

gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857
0 66.51326044311185
0 10018754.1713946 0
45 40.97989806962013
5009377.08569731 5009377.08569731 0

My coordinates in meters are:

0 10018754.1713946 5009377.08569731 5009377.08569731 (West North East South)

Now I can use gdal_translate to generate a geotiff:

gdal_translate -of Gtiff -co tfw=yes -a_ullr 0 10018754.1713946 5009377.08569731 5009377.08569731 -a_srs EPSG:3857 4_2.png t3857.tif

My impression is that this is not correct because the borders of the maps are shifted toward north. It seems to be the right idea.

---

Third idea: EPSG:3857 through EPSG:4055
I read that "web mercator" uses WGS84 coordinates but consider them as if they where spherical coordinates. Due to the difference between a geodetic and a geocentric latitude (See Wikipedia about the latitude), the latitude values will not be the same on an ellipsoid or on a sphere. I found that EPSG:4055 is the code for spherical coordinates on a sphere based on WGS84.

Converting the coordinates to EPSG:4055:

gdaltransform -s_srs EPSG:4326 -t_srs EPSG:4055
0 66.51326044311185
0 66.3722684317026 -17964.0621483233
45 40.97989806962013
45 40.7894557844857 -9152.84527519904

The corresponding spherical coordinates are then:

0 66.3722684317026 45 40.7894557844857 (West North East South)

Then I do as if those coordinates where still on the ellipsoid (EPSG:4326) and convert them to web mercator:

gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857
0 66.3722684317026
0 9979483.26733298 0
45 40.7894557844857
5009377.08569731 4981335.86590183 0

The resulting coordinates differ from the one by idea2:

0 9979483.26733298 5009377.08569731 4981335.86590183 (West North East South)

Now I just have to write the coordinates into the map:

gdal_translate -of Gtiff -co tfw=yes -a_ullr 0 9979483.26733298 5009377.08569731 4981335.86590183 -a_srs EPSG:3857 4_2.png t3857_through_4055.tif

This third idea seems to give the best results. But I am not sure if it is correct. If idea 3 is correct, is there an EPSG code to do this operation in one step?

---

Fourth idea: EPSG:3857 through towgs84=0,0,0,0,0,0,0

gdal (and apparently epsg too) defines EPSG:3857 like that:

+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

whereas spatialreference.org like that:

+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

If I use the definition from spatialreference.org I got the correct coordinates in one step (Well I still don't if they are the "correct" coordinates but at least they are the sames as by idea 3):

gdaltransform -s_srs EPSG:4326 -t_srs "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
0 66.51326044311185
0 9979483.26733298 -17964.0621483233
45 40.97989806962013
5009377.08569731 4981335.86590183 -9152.84527519904

Why is there such a difference in the definitions of EPSG:3857?

Answer 1

The tile image is already in EPSG:3857. Why not just create a world file to reference it?

http://en.wikipedia.org/wiki/World_file

For the tile that covers N. America at zoom 1, you'd be looking at the following worldfile contents:

78271.517
0
0
-78271.517
-19998372.6
19998372.6

Where those numbers came from:

• Line 1: width of an image pixel in world coordinates = 20037508.342789244 meters / 256 pixels.
• Line 2 and 3: rotation, so n/a.
• Line 4: height of an image pixel in world coordinates. Same as line 1 but negative, because in image files increasing y corresponds to 'down' while in the coordinate system, increasing y corresponds to 'up'.
• Line 5: X coordinate in world coordinates of the center of the top-left pixel. This is -20037508.342789244, as reported by the tiles a la carte link, plus 1/2 of a pixel to bring it to the center.
• Line 6: Ditto, only Y coordinate of the top left.

GDAL ought to pick up your worldfile (.pgw for the png); you'll still have to tell it EPSG:3857 on the command line.

(Note: didn't have time to test this, so it's all off the cuff... but hopefully correct on the first try anyway! ;)

Answer 2

Functions necessary for generation of global tiles used on the web. It contains classes implementing coordinate conversions for:

• GlobalMercator (based on EPSG:900913 = EPSG:3785) for Google Maps, Yahoo Maps, Microsoft Bing Maps compatible tiles

• GlobalGeodetic (based on EPSG:4326) for OpenLayers Base Map and Google Earth compatible tiles

http://svn.osgeo.org/gdal/sandbox/klokan/gdal2tiles.py

Answer 3

About my secondary question in the fourth ideas: Why is there such a difference in the definitions of EPSG:3857 between the definition in gdal and on spatialreference.org:

The main difference is that gdal uses "+nadgrids=@null" and spatialreference "+towgs84=0,0,0,0,0,0,0". According to the documentation or the PROJ.4 format both parameter are used for datum transforms. I found an interesting comment from Mikael Rittri on the Proj4 listserver:

"The +to definition you suggest is not correct for WGS84 Web Mercator, though.
This is because the datum shift you use, 

   +towgs84=0,0,0,0,0,0,0

is not the same as unmodified lat/long, or "without any datum shift" as you say.
It means "no datum shift" only if the +from and +to definitions use the same ellipsoid,
and in your case the +from definition uses the WGS84 ellipsoid, while the +to definition
uses a sphere instead.  

So, you have to use a trick: use 

   +nadgrids=@null 

instead."

The use of "+towgs84=0,0,0,0,0,0,0" doesn't seems to be right or at least only under specific conditions.