Bagaimana cara georeferensi ubin web mercator dengan benar menggunakan gdal?

Sebagai contoh, saya akan mengambil petak berikut http://a.tile.openstreetmap.org/3/4/2.png dan menyimpannya sebagai "4_2.png".

Koordinat WGS84 ubin ini dapat dihitung atau dibaca di sana dengan mengklik ubin yang sesuai:

0 66.51326044311185 45 40.97989806962013 (West North East South)

Cara melakukan georeferensi ubin dengan benar (menggunakan gdal untuk menghasilkan geotiff atau format georeferensi lain) sehingga:

• Bitmap tidak perlu diregangkan (= piksel dalam geotiff persis sama dengan di bitmap asli)
• Gambar yang dihasilkan akan dibuka di tempat yang tepat di penampil / editor GIS (seperti misalnya di TatukGIS Free Viewer )?

(Diedit pada 19 September 2011 untuk memperjelas pertanyaan saya dan menyertakan kesimpulan saya)

Kesimpulan saya:

Saya pertama kali berpikir bahwa ide ketiga (lihat di bawah) adalah yang benar. Saya membuka geotiff di Penampil GIS dan membandingkan koordinat yang ditampilkan dengan apa yang saya harapkan. Geotiff dari ide kedua tampaknya bergeser 2 piksel ke arah utara. Itu sebabnya saya menganggap ide 3 (atau 4) sebagai yang benar.
Tetapi jika Anda mencoba dengan ubin di tingkat zoom yang jauh lebih tinggi, geotiff dari ide 3 secara definitif bergeser ke selatan. Sungguh konyol membandingkan koordinat pada ubin tingkat zoom 3. Batas negara pada tingkat zoom seperti itu disederhanakan sehingga perbandingannya tidak memberikan hasil yang baik.

Dan S. benar, gambar ubin sudah ada di EPSG: 3857. Gagasan kedua adalah yang benar (dan memberikan hasil yang baik pada level zoom tinggi juga)

Gagasan pertama: EPSG: 4326
Kode EPSG untuk koordinat WGS84 adalah EPSG: 4326 . Jadi saya hanya menggunakan koordinat WGS84 untuk georefence ubin sebagai geotif menggunakan 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

Peta yang dihasilkan ditampilkan di tempat yang tepat tetapi saya khawatir bahwa proyeksi itu tidak benar dan mungkin ada pergeseran di tengah ubin. Setelah mencoba cukup lama untuk memeriksa bahwa dengan memproyeksi ulang peta dengan gdalwarp, saya mengunduh versi demo Global Mapper dan ini tampaknya menjadi masalah (sama saja dengan ide 3 tetapi ada pergeseran di dalam ubin). Gambar harus dikencangkan agar dapat menggunakan koordinat EPSG: 4326.

Gagasan kedua: EPSG: 3857
Ubin ini menggunakan proyeksi "web mercator" (alias proyeksi google map), yang sekarang memiliki kode EPSG: EPSG: 3857 (alias EPSG: 900913). Saya cukup mengonversi koordinat menggunakan 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

Koordinat saya dalam meter adalah:

0 10018754.1713946 5009377.08569731 5009377.08569731 (West North East South)

Sekarang saya bisa menggunakan gdal_translate untuk menghasilkan 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

Kesan saya adalah ini tidak benar karena perbatasan peta bergeser ke utara. Sepertinya itu ide yang tepat.

Gagasan ketiga: EPSG: 3857 hingga EPSG: 4055
Saya membaca bahwa "web mercator" menggunakan koordinat WGS84 tetapi menganggapnya seolah-olah koordinat koordinat bola. Karena perbedaan antara garis lintang geodetik dan geosentris (Lihat Wikipedia tentang garis lintang ), nilai garis lintang tidak akan sama pada ellipsoid atau pada bola. Saya menemukan bahwa EPSG: 4055 adalah kode untuk koordinat bola pada bola berdasarkan WGS84.

Mengonversi koordinat ke 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

Koordinat bola yang sesuai adalah:

0 66.3722684317026 45 40.7894557844857 (West North East South)

Lalu saya melakukan seolah-olah koordinat itu di mana masih di ellipsoid (EPSG: 4326) dan mengubahnya menjadi 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

Koordinat yang dihasilkan berbeda dari yang oleh idea2:

0 9979483.26733298 5009377.08569731 4981335.86590183 (West North East South)

Sekarang saya hanya perlu menulis koordinat ke peta:

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

Gagasan ketiga ini tampaknya memberikan hasil terbaik. Tetapi saya tidak yakin apakah itu benar. Jika ide 3 benar, apakah ada kode EPSG untuk melakukan operasi ini dalam satu langkah?

Gagasan keempat: EPSG: 3857 hingga towgs84 = 0,0,0,0,0,0,0,0

gdal (dan tampaknya juga epsg) mendefinisikan EPSG: 3857 seperti itu:

+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

sedangkan spatialreference.org seperti itu:

+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

Jika saya menggunakan definisi dari spatialreference.org saya mendapatkan koordinat yang benar dalam satu langkah (Yah saya masih tidak jika mereka adalah koordinat "benar" tetapi setidaknya mereka sama seperti dengan ide 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

Mengapa ada perbedaan definisi EPSG: 3857?

Gambar ubin sudah ada di EPSG: 3857. Mengapa tidak membuat file dunia untuk referensi saja?

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

Untuk ubin yang mencakup N. America pada zoom 1, Anda akan melihat konten worldfile berikut:

78271.517
0
0
-78271.517
-19998372.6
19998372.6

Dari mana angka-angka itu berasal:

• Baris 1: lebar piksel gambar dalam koordinat dunia = 20037508.342789244 meter / 256 piksel.
• Baris 2 dan 3: rotasi, jadi n / a.
• Baris 4: ketinggian piksel gambar dalam koordinat dunia. Sama seperti baris 1 tetapi negatif, karena dalam file gambar meningkat y sesuai dengan 'turun' sementara dalam sistem koordinat, meningkatkan y sesuai dengan 'naik'.
• Jalur 5: X mengoordinasikan koordinat dunia dari pusat piksel kiri-atas. Ini adalah -20037508.342789244, seperti yang dilaporkan oleh ubin a la carte link, ditambah 1/2 piksel untuk membawanya ke tengah.
• Jalur 6: Ditto, hanya Y yang berkoordinasi dari kiri atas.

GDAL harus mengambil file worldfile Anda (.pgw untuk png); Anda masih harus mengatakannya EPSG: 3857 di baris perintah.

(Catatan: tidak punya waktu untuk menguji ini, jadi semuanya batal ... tapi semoga tetap benar pada percobaan pertama!;)

Fungsi yang diperlukan untuk pembuatan ubin global yang digunakan di web. Ini berisi kelas yang menerapkan konversi koordinat untuk:

• GlobalMercator (berdasarkan EPSG: 900913 = EPSG: 3785 ) untuk ubin yang kompatibel dengan Google Maps, Yahoo Maps, Microsoft Bing Maps

• GlobalGeodetic (berdasarkan EPSG: 4326) untuk OpenLayers Base Map dan ubin yang kompatibel dengan Google Earth

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

Tentang pertanyaan kedua saya dalam ide keempat: Mengapa ada perbedaan definisi EPSG: 3857 antara definisi di gdal dan di spatialreference.org :

Perbedaan utama adalah bahwa gdal menggunakan "+ nadgrids = @ null" dan spatialreference "+ towgs84 = 0,0,0,0,0,0,0,0". Menurut dokumentasi atau format PROJ.4 kedua parameter digunakan untuk transformasi datum. Saya menemukan komentar menarik dari Mikael Rittri pada server daftar 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."

Penggunaan "+ towgs84 = 0,0,0,0,0,0,0,0" tampaknya tidak benar atau setidaknya hanya dalam kondisi tertentu.