Jak poprawnie georeferencyjnie kafelek mercatora internetowego za pomocą gdal?

Jako przykład wezmę następujący kafelek http://a.tile.openstreetmap.org/3/4/2.png i zapiszę go jako „4_2.png”.

Współrzędne WGS84 tej płytki można obliczyć lub czytać tam klikając na odpowiedni żeton:

0 66.51326044311185 45 40.97989806962013 (West North East South)

Jak poprawnie dokonać georeferencji kafelka (używając gdal do wygenerowania geotiffu lub innego formatu georeferencji), aby:

• Mapa bitowa nie musi być rozciągana (= piksele w geotiffie są dokładnie takie same jak w oryginalnej bitmapie)
• Powstały obraz zostanie otwarty we właściwym miejscu w przeglądarce / edytorze GIS (jak na przykład w TatukGIS Free Viewer )?

(Edytowane 19 września 2011 r., Aby wyjaśnić moje pytanie i dołączyć moje wnioski)

Mój wniosek:

Najpierw pomyślałem, że trzeci pomysł (patrz poniżej) był właściwy. Otworzyłem geotiff w przeglądarce GIS i porównałem wyświetlane współrzędne z oczekiwanymi. Geotiff z drugiego pomysłu wydaje się przesunięty o 2 piksele w kierunku północnym. Dlatego uważałem pomysł 3 (lub 4) za właściwy.
Ale jeśli spróbujesz z kafelkiem na znacznie wyższym poziomie powiększenia, geotiff z pomysłu 3 zostanie ostatecznie przesunięty w kierunku południowym. Głupio było porównywać współrzędne na kafelku poziomu powiększenia 3. Granice kraju na takim poziomie powiększenia muszą być uproszczone, aby porównanie nie dawało dobrych wyników.

Dan S. miał rację, obraz kafelka jest już w EPSG: 3857. Drugi pomysł jest wtedy właściwy (i daje dobre wyniki również przy wysokich poziomach powiększenia)

Pierwszy pomysł: EPSG: 4326
Kod EPSG dla współrzędnych WGS84 to EPSG: 4326 . Po prostu używam współrzędnych WGS84 do georefence kafelka jako geotif przy użyciu 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

Powstała mapa jest wyświetlana we właściwym miejscu, ale obawiam się, że rzut nie jest prawidłowy i że może nastąpić przesunięcie na środku kafelka. Po długich próbach sprawdzenia, czy poprzez ponowne zaprojektowanie mapy za pomocą gdalwarp, pobrałem wersję demo Global Mappera i wydaje się, że tak jest (wygląda na granice jak w idei 3, ale przesunięcie wewnątrz płytki). Obraz powinien zostać rozciągnięty, aby można było użyć współrzędnych EPSG: 4326.

Drugi pomysł: EPSG: 3857
Ten kafelek korzysta z projekcji „web mercator” (alias google map map), która ma teraz kod EPSG: EPSG: 3857 (alias EPSG: 900913). Po prostu przekształcam współrzędne za pomocą 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

Moje współrzędne w metrach to:

0 10018754.1713946 5009377.08569731 5009377.08569731 (West North East South)

Teraz mogę użyć gdal_translate do wygenerowania 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

Mam wrażenie, że to nieprawda, ponieważ granice map są przesunięte w kierunku północnym. Wydaje się, że to właściwy pomysł.

Trzeci pomysł: EPSG: 3857 do EPSG: 4055
Czytałem, że „mercator sieci” używa współrzędnych WGS84, ale uważam je za współrzędne sferyczne. Ze względu na różnicę między szerokością geodezyjną i geocentryczną (patrz Wikipedia o szerokości geograficznej ), wartości szerokości geograficznej nie będą takie same na elipsoidzie ani na kuli. Odkryłem, że EPSG: 4055 jest kodem współrzędnych sferycznych na kuli opartej na WGS84.

Konwersja współrzędnych na 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

Odpowiednie współrzędne sferyczne to:

0 66.3722684317026 45 40.7894557844857 (West North East South)

Następnie robię tak, jakby te współrzędne były jeszcze na elipsoidzie (EPSG: 4326) i przekształcam je w mercator sieci:

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

Uzyskane współrzędne różnią się od współrzędnych idea2:

0 9979483.26733298 5009377.08569731 4981335.86590183 (West North East South)

Teraz muszę tylko zapisać współrzędne na mapie:

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

Ten trzeci pomysł wydaje się dawać najlepsze wyniki. Ale nie jestem pewien, czy to prawda. Jeśli pomysł 3 jest poprawny, czy istnieje kod EPSG umożliwiający wykonanie tej operacji w jednym kroku?

Czwarty pomysł: EPSG: 3857 przez towgs84 = 0,0,0,0,0,0,0,0

gdal (i najwyraźniej również epsg) definiuje EPSG: 3857 w następujący sposób:

+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

podczas gdy spatialreference.org:

+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

Jeśli użyję definicji ze spatialreference.org, w jednym kroku uzyskałem prawidłowe współrzędne (cóż, nadal nie wiem, czy są to „poprawne” współrzędne, ale przynajmniej są one samami jak w idei 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

Dlaczego jest taka różnica w definicjach EPSG: 3857?

Obraz kafelka jest już w formacie EPSG: 3857. Dlaczego nie utworzyć pliku świata, aby się do niego odwoływać?

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

W przypadku kafelka obejmującego Amerykę Północną przy powiększeniu 1 zobaczysz następującą zawartość pliku światowego:

78271.517
0
0
-78271.517
-19998372.6
19998372.6

Skąd te liczby pochodzą:

• Linia 1: szerokość piksela obrazu we współrzędnych światowych = 20037508.342789244 metrów / 256 pikseli.
• Wiersze 2 i 3: obrót, więc nie dotyczy.
• Linia 4: wysokość piksela obrazu we współrzędnych światowych. To samo co linia 1, ale ujemne, ponieważ w plikach graficznych zwiększenie y odpowiada „w dół”, podczas gdy w układzie współrzędnych, zwiększenie y odpowiada „w górę”.
• Linia 5: Współrzędna X we współrzędnych światowych środka lewego górnego piksela. Jest to -20037508.342789244, jak podają kafelki link a la carte, plus 1/2 piksela, aby przenieść go na środek.
• Wiersz 6: Ditto, tylko współrzędna Y lewego górnego rogu.

GDAL powinien odebrać twój plik świata (.pgw dla png); nadal musisz powiedzieć EPSG: 3857 w wierszu polecenia.

(Uwaga: nie miałem czasu, aby to przetestować, więc to wszystko jest poza mankietem ... ale mam nadzieję, że poprawi się za pierwszym razem!)

Funkcje niezbędne do generowania globalnych kafelków używanych w sieci. Zawiera klasy implementujące konwersje współrzędnych dla:

• GlobalMercator (na podstawie EPSG: 900913 = EPSG: 3785 ) dla Google Maps, Yahoo Maps, kafelków zgodnych z Microsoft Bing Maps

• GlobalGeodetic (na podstawie EPSG: 4326) dla mapy podstawowej OpenLayers i kafelków zgodnych z Google Earth

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

O moim drugorzędnym pytaniu w czwartych pomysłach: Dlaczego jest taka różnica w definicjach EPSG: 3857 między definicją w gdal i na spatialreference.org :

Główną różnicą jest to, że gdal używa „+ nadgrids = @ null” i odniesienia przestrzennego „+ towgs84 = 0,0,0,0,0,0,0,0”. Zgodnie z dokumentacją lub formatem PROJ.4 oba parametry są używane do transformacji układu odniesienia. Znalazłem interesujący komentarz Mikaela Rittri na serwerze list Proj4 :

"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."

Użycie „+ towgs84 = 0,0,0,0,0,0,0,0” wydaje się niewłaściwe lub przynajmniej tylko w określonych warunkach.