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
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! ;)
No comments:
Post a Comment