Cholnoky Jenő hagyatékából: Dunavölgy átnézetes térképe

Egyik kedvenc blogomon, a Dunai Szigeteken volt egy blogbejegyzés Vályi Béla vízrajzi térképeiről, ezek a  térképek szabadon letölthetőek a webről viszonylag jó felbontásban. Itt található egy nagyon jó leírás arról mit és hogyan kell követni a georeferálásukhoz. A módszer lényege az, hogy a térképlapok sarokpontjainak a koordinátái meghatározhatóak, illetve a térképlapok vetülete is közelíthető, miáltal a térképlapok áttranszformálhatóak modern vetületekbe, megjeleníthetőek akár google maps összehasonlítással is.
A feladatot a PROJ4 illetve a GDAL könyvtár segítségével, python nyelv használatával oldottam meg. Mielőtt belekezdenék mutatom milyen is a végeredmény: 
Teljes méret

1. lépés: Meghatározni a térképlapok sarkainak x,y képi koordinátáit, és az össztérképen belül melyik sorban, oszlopban található.
2. lépés: Kiszámolni a gcp pontokat és meghatározni az alkalmazott vetületet
3. lépés: Átforgatni a térképlapokat, levágni a szélüket, összeilleszteni, majd csempékre vágni.
4. lépés: Ábrázolni webes felületen

1. lépés: Google docs táblázatban összegyűjtöttem a szükséges információkat, a térképek sarokpontjait, fontos, hogy a képnézegető a bal felső sarkot olvassa be 0,0 koordinátákkal.

1. ábra: jobb felső koordináta leolvasása (vajon a szélen
kála és vaníliavirág található???)

2. lépés: PROJ4 formátumban "+" kapcsolók segítségével definiáltam a  vetületet+dátumot+forgási ellipszoidot+átváltási paramétereket stb.  2. oszlop, 2. sorának térképlapjára például így néz ki ez:

+proj=sinu  -> Szinuszoidális vetület
+ellps=bessel -> Bessel 1840 (vagy 41??) ellipszoid
+lon_0=35d15' -> Kezdő meridián a térképlap közepén áthaladó hosszúsági körön kell, hogy legyen
+lat_0=0.0 -> Hossztartó paralellkör az egyenlítőn
+pm=-17.66278333 -> Ferro1
+towgs84=600,205,437 -> WGS84 hez képest az eltolási paraméterek (nem szükséges megadni)
+units=meters -> Mértékegység méter
+no_defs                           

Mivel a sarokpontok földrajzi fokban vannak megadva, a gcp pontok miatt ezeket át kell válatni:

gdaltransform.exe -s_srs "+proj=longlat +ellps=bessel +pm=-17.66278333" -t_srs "+proj=sinu  +ellps=bessel +pm=-17.66278333 units=meters"

gdal_translate -a_srs "+proj=sinu +ellps=bessel  +units=meters +lon_0=35d15' +lat_0=0.0  +pm=-17.66278333 +towgs84=600,205,437 +no_defs" -gcp 907 441 -37217.9265 5331782.6106 -gcp 8465 440 37217.9265 5331782.6106 -gcp 8556 7401 37667.2082 5262298.7502 -gcp 909 7496 -37667.2082 5262298.7502 "D:\Elkeszitett_munkak\DV\2M-VIZ (865) Duna06.jpg" "D:\Elkeszitett_munkak\DV\teszt1.tiff" 

 +lon_0 paramétert mindig a térképlap középvonalához kell igazítani (33 fok 15 perc + oszlopszám*1 fok). A hosszúság koordináták -0,5 fok és +0,5 fok.

gdalwarp -tps -t_srs "EPSG:3857"  "D:\Elkeszitett_munkak\DV\1.tiff"  "D:\Elkeszitett_munkak\DV\g1.tiff" 

Persze az egészet érdemes automatizálni. GDAL/OGR könyvtár python interfésze erre pont alkalmas.

Forrás:
Mészáros János: Vályi Béla vízrajzi térképei és geoinformatikai rendszerbe integrálásuk