Bardzo lubię mapy. Pewnego dnia wcześniej w tym roku wpadłem na pomysł, aby ozdobić ścianę drogową mapą Polski. W pomyśle był jeden mały haczyk: nie chciałem byle jakiej mapy, chciałem mapę Google. Szybki przegląd internetu nie przyniósł żadnych konkretnych rozwiązań – super, będę pierwszy!

Mapy Google wyświetlają różną ilość szczegółów w zależności od stopnia przybliżenia. Gdy na ekranie mieści się cała Polska, szczegółów nie ma wcale. Moja mapa ścienna miała być duża (planowałem coś w okolicach metr na metr), więc potrzebowałem szczegółów. Zdecydowałem się na napisanie aplikacji w C#, której zadaniem było wyświetlanie kolejnych fragmentów mapy i robienie screenshota. Rozwój aplikacji zakończył się w momencie, gdy zadziałała, więc nawet guziki nie były podpisane:

Dla średniej szczegółowości wystarczyło pokrycie Polski za pomocą ~120 screenów, a dla wysokiej zrobiłem ich ponad 1200.

Samo robienie screenów to łatwy fragment. Trudniej było skleić je w całość. Pamiętajmy, że ziemia jest kulą, więc samo sklejenie obrazków spowoduje powstanie dziur, ponieważ relacja pomiędzy stopniami długości geograficznej a kilometrami zależy od szerokości geograficznej. W sieci łatwo znalazłem relację pomiędzy pixelami mapy a fizycznymi odległościami w zależności od szerokości geograficznej:

Meters per pixel = 156543.03392 * Math.cos(latLng.lat() * Math.PI / 180) / Math.pow(2, zoom)

W kolejnym kroku wycinałem z moich screenów fragment mapy bez elementów interfejsu – dla ułatwienia zawsze taki sam prostokąt na środku ekranu. Robiłem to oczywiście automatycznie za pomocą polecenia convert. Przy okazji przycinania konwertowałem wycięty fragment na format tiff. Mając obliczone z wcześniejszego wzoru współrzędne geograficzne rogów fragmentu mapy konwertowałem plik graficzny tiff do formatu geotiff za pomocą polecenia gdal_translate. Geotiff to nic innego jak zwykły tiff z dodatkowymi informacjami dotyczącymi jego orientacji w przestrzeni. Ten proces był sterowany skryptem w MATLABie:

lat_min = 48.7;
lat_max = 55;
lon_min = 13.8;
lon_max = 24.6;

lat_step = .6;
lon_step = 1.2;


xof = 100;
yof = 100;
x = 1920 - 2*xof;
y = 982 - 2*yof;

for lat = lat_min:lat_step:lat_max
	for lon = lon_min:lon_step:lon_max
		im_width = 1720;
		im_height = 782;
		zoom = 10;

		meters_vertical = (156543.03392 * cos(deg2rad(lat))  / (2^zoom));

		[lat_tm, lon_tm] = reckon(lat, lon, km2deg(meters_vertical*(im_height/2)/1000), 0); 
		[lat_bm, lon_bm] = reckon(lat, lon, km2deg(meters_vertical*(im_height/2)/1000), 180); 

		meters_horizontal_top = 156543.03392 * cos(deg2rad(lat_tm)) / (2^zoom);
		meters_horizontal_bottom = 156543.03392 * cos(deg2rad(lat_bm)) / (2^zoom);

		[lat_tr, lon_tr] = reckon(lat_tm, lon_tm, km2deg(meters_horizontal_top*(im_width/2 - .5)/1000), 90); 
		[lat_tl, lon_tl] = reckon(lat_tm, lon_tm, km2deg(meters_horizontal_top*(im_width/2 - .5)/1000), 270); 

		[lat_br, lon_br] = reckon(lat_bm, lon_bm, km2deg(meters_horizontal_bottom*(im_width/2 - .5)/1000), 90); 
		[lat_bl, lon_bl] = reckon(lat_bm, lon_bm, km2deg(meters_horizontal_bottom*(im_width/2 - .5)/1000), 270); 

		system(sprintf('convert source/mapa_%8.5f_%8.5f.png -crop %dx%d+%d+%d -resize %dx%d -transparent white crop/mapa_%8.5f_%8.5f.tiff', lat, lon, x, y, xof, yof, 2*x, 2*y, lat, lon));
		system(sprintf('gdal_translate -of GTiff -a_srs WGS84 -a_ullr  %8.5f %8.5f %8.5f %8.5f crop/mapa_%8.5f_%8.5f.tiff crop/geomapa_%8.5f_%8.5f.tiff', lon_tl, lat_tl, lon_br, lat_br, lat, lon, lat, lon));
	end
end

Następnie za pomocą skryptu gdal_merge.py połączyłem wszystkie geotiffy w jednego wielkiego geotiffa. Tutaj praca mogłaby się zakończyć, ale chciałem, żeby moja mapa miała jakieś fajne odwzorowanie, ramkę itp. Najlepszym znanym mi programem do rysowania map jest pakiet GMT (http://gmt.soest.hawaii.edu/). Za pomocą prostego skryptu przygotowałem pdf (ps to prawie pdf) z mapą do wydruku:

gmt gmtset PS_MEDIA=Custom_60ix80i
gmt gmtset MAP_FRAME_WIDTH=.25i
gmt gmtset FONT_ANNOT_PRIMARY=32p,Helvetica,black
gmt gmtset COLOR_NAN=white
gmt psbasemap -R13.8/24.5/48.7/55 -Jl19.20/19.30/52/60/1:600000 -B1g.5:."": -K -P > mapa_all.ps
gmt grdimage ALL.tiff -D -R13.8/24.5/48.7/55 -Jl19.20/19.30/52/60/1:600000 -Gf -Q -P -O -K >> mapa_all.ps
gmt pscoast -R -Df -J -N1/3p,- -A0/0/1  -O -P >> mapa_all.ps

Tak wyglądał efekt końcowy po wydrukowaniu:

Nie dzielę się plikiem z mapą, bo jest bardzo duży i nie jestem pewnie czy google nie chciałoby mi urwać za to głowy.