User:Maxbe/Messtischblätter aus OSM-Tiles
Worum gehts?
Im Forum wurde gefragt, ob man Messtischblätter aus der OpenTopoMap erzeugen kann. Freundlicherweise durfte ich einen Cache für Bayern für diese Karte aufsetzen und damit basteln.
Das Ergebnis sind Blätter wie das Bild rechts: Eine Ausschnitt aus der OpenTopoMap, in passende Projektion gestaucht und gedreht mit dem Blattschnitt des Messtischblattes Nr. 6332 und einer definierten Auflösung von 15 Meter/Pixel. Wenn das Bild kein PNG wäre, sondern ein geoTIFF, wäre es auch korrekt referenziert...
Blattschnitt
Messtischblätter haben eine feste Blatteinteilung: Das Blatt 0000 liegt irgendwo in der Nordsee mit der oberen Kante bei 56° Nord und der linken Kante bei 5° 40' Ost, das Blatt 0001 rechts daneben, das Blatt 0100 darunter. Ein Blatt ist 6 Gradminuten hoch und 10 Gradminuten breit. Das führt dann zu einem fast gleichmässigen Gitter von ca. 10x10km. Die Nummerierung ist auch einfach: die ersten beiden Ziffern geben die Zeile an, die letzten beiden Ziffern die Spalte.
Dummerweise ist diese Einteilung recht alt und die schönen glatten Gradangaben stimmen auf dem Bessel-Ellipsoid. Wenn man neuere Projektionen verwendet liegen die Breiten- und Längengrade ein paar dutzend Meter daneben. Man muss also die Einteilung erst in "alten" Graden vornehmen und dann das Gitter umwandeln. Wer sowas braucht, findet auf der unten verlinkten Seite auch ein Shapefile mit umgerechnetem Gitter und passender Nummerierung.
Ausgabe in UTM oder Gauß-Krüger
Das Bild mit dem Blattschnitt da oben ist eine Mercatorprojektion, in der so ziemlich alle frei verfügbaren Karten aus OSM ausgeliefert werden. Dummerweise mögen Landvermesser und offensichtlich auch Pilz- und Insektenkartierer lieber transversale Mercatorprojektionen wie UTM oder Gauß-Krüger. Die haben einige Vorteile auf kleinen Gebieten, wie z.B. ein rechtwinkliges Koordinatensystem und einen Maßstab, der an jeder Stelle des Kartenblattes gleich ist. Dafür decken sie auch immer nur einen einige 100km breiten Streifen auf der Erde ab, links und rechts davon muss man das System wechseln. Längen- und Breitengrade liegen in diesen Systemen näherungsweise als Trapeze vor. Sind sie in echt ja auch: Am Nordpol müssen sich die Längengrade alle treffen.
Das Gitter der Messtischblätter liegt also schief auf den UTM oder GK-Karten und auch die Karte muss entsprechend gedreht werden. Ausserdem muss sie ein bisschen gestaucht oder gedehnt werden, um bei einer festgelegten Auflösung ins Gitter zu passen. Die 17-20 Zoomstufen unserer Karten müssen stufenlos umgerechnet werden können, um auch "Zwischengrössen" zu erzeugen. Zwischengrössen kann man nicht vermeiden, weil in der Mercatorprojektion die Grösse eines Pixels (also die "Meter pro Bildpunkt") innerhalb eines Kartenausschnittes variiert. Am Südende ist die Auflösung geringer als an der oberen Kante.
Einschub: Warum das alles ein bisschen gepfuscht ist...
Wenn man das richtig machen wollte, würde man die Karte nicht in ein paar Zoomstufen rendern und dann stauchen und drehen, sondern gleich in der richtigen Auflösung, dem richtigen Ausschnitt und der richtigen Projektion erzeugen. Damit verhindert man, dass das Kartenbild beim umrechnen matschig wird, hätte alle Beschriftungen in lesbarer Grösse korrekt ausgerichtet statt mitgedreht und gestaucht und würde am Kartenrand keine halben Beschriftungen stehen lassen.
Leider steht mir aber kein Renderer zur Verfügung, der nach belieben solche Karten erzeugen kann.
Zum Nachbauen
Zutaten
Die Auswahl der Zutaten hat sich vor allem daran orientiert, was schon da war. Vermutlich kann man einiges davon auch durch anderes ersetzen.
- Ein Mapproxy ist so ziemlich das wichtigste dabei. Der kann Karten von einem Tileserver holen, zwischenspeichern und stellt damit einen WMS-Dienst zur Verfügung, an dem man Kartenausschnitte in allen möglichen Projektionen abholen kann. Er sucht sich dann aus (oder man konfiguriert es), aus welchen Kacheln welcher Zoomstufen der gewünschte Ausschnitt zusammengesetzt wird. Holt sich die, setzt das Bild zusammen, projiziert das Bild in die Wunschprojektion und liefert es aus.
- Ein Mapserver kann diesen WMS abfragen, Shapefiles (die Blattgrenzen) drüberlegen, zusätzliche Beschriftungen dazumalen, einen Maßstab ausgeben und alle möglichen Ausgabegformate wie png, jpeg, svg oder geotiff erzeugen. Der macht das fertige Bild.
- PostgreSQL und PostGIS wird zum Berechnen von Koordinaten verwendet. Z.B. müssen die Eckpunkte der schief sitzenden Blätter in die Zielprojektion umgerechnet werden. Vermutlich ist das der Teil, den man auch anders gut lösen konnte, aber wenn man schon eine Datenbank im Haus hat, kann man sie auch gut die Blattschnitte verwalten und umrechnen lassen.
Dazu kommen noch ein paar Werkzeuge und Bibliotheken, auf die diese drei zurückgreifen (proj.4, gd, gdal ...). Ich kann leider nicht sagen, was wo mit installiert wird und deshalb kein vollständiges Rezept abliefern. Ich schreibe deshalb nur, was ich für die Messtischblätter zusätzlich eingebaut habe... Ausserdem verwende ich den Tileserver für den OSM-Mapnik-Stil auf osm.org als Beispiel.
Mapproxy
Dem muss man beibringen, dass er sämtliche Projektionen ausgibt, die man möchte. Nach der Installation kann er nur wenige und z.B. GK4/EPSG:31467 ist nicht dabei (alles hinter srs:). Ausserdem muss er recht grosse Bilder ausliefern (max_output_pixels)
services:
demo:
tms:
use_grid_names: true
origin: 'nw'
wms:
versions: ['1.1.1']
srs: ['EPSG:4326', 'EPSG:31466', 'EPSG:31467', 'EPSG:900913', 'EPSG:3857']
proj_data_dir: "/usr/share/proj/"
image_formats: ['image/png', 'image/jpeg']
max_output_pixels: [20000,20000]
Ausserdem muss man den Layer für OSM hinzufügen
layers: - name: osm title: osm sources: [osm_cache]
Den Cache dafür angeben
caches:
osm_cache:
grids: [mercator]
sources: [osm_source]
max_tile_limit: 1000000
cache:
type: file
directory_layout: tc
directory: /var/spool/tilecache/osm
und ihm sagen, wo er die Kacheln abholen soll:
sources:
osm_source:
type: tile
grid: osm_grid
concurrent_requests: 1
url: http://c.tile.openstreetmap.org/%(tms_path)s.%(format)s
image:
resampling_method: bicubic
grids:
mercator:
base: GLOBAL_MERCATOR
num_levels: 20
proj.4
Bei den Services des Mapproxy steht ein "proj_data_dir". Das gibt das Verzeichnis an, in dem Mapproxy eine Datei "epsg" mit den Definitionen der Projektionen sucht und die "EPSG-Nummern" zuordnet. Dort kann man zusätzliche Wunschprojektionen eintragen und mit den korrekten proj.4-Parametern versehen. Auch bei bereits vordefinierten EPSG-Codes können dort Änderungen vorgenommen werden. Falls man z.B. auf die Verwendung von NTv2 Wert legt, zeigt man dem Programm dort den Weg zur Datei "BETA2007.gsb" mit den Shiftwerten:
<25831> +proj=utm +zone=31 +ellps=GRS80 +units=m +datum=WGS84 +no_defs <> <25832> +proj=utm +zone=32 +ellps=GRS80 +units=m +datum=WGS84 +no_defs <> <25833> +proj=utm +zone=33 +ellps=GRS80 +units=m +datum=WGS84 +no_defs <> <25834> +proj=utm +zone=34 +ellps=GRS80 +units=m +datum=WGS84 +no_defs <> <31467> +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +nadgrids=/usr/share/proj/BETA2007.gsb +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.70 +units=m <31466> +proj=tmerc +lat_0=0 +lon_0=6 +k=1 +x_0=2500000 +y_0=0 +ellps=bessel +nadgrids=/usr/share/proj/BETA2007.gsb +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.70 +units=m <31468> +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +nadgrids=/usr/share/proj/BETA2007.gsb +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.70 +units=m <31469> +proj=tmerc +lat_0=0 +lon_0=15 +k=1 +x_0=5500000 +y_0=0 +ellps=bessel +nadgrids=/usr/share/proj/BETA2007.gsb +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.70 +units=m
Mapserver
Das Mapfile ist recht schlicht: Ein Layer für die Karte selbst, ein Layer für die weissen Rahmen der Nachbarblätter, einer für die Bauchbinde mit dem Copyright und einer für den Maßstabsbalken:
#####################################################
# Mapfile fuer die Erzeugung von Blattschnittkarten #
#####################################################
MAP
#
# Alle anderen Ausgabeformate sind default. Bei gtiff kommt FORMATOPTION "COMPRESS=DEFLATE" dazu.
#
OUTPUTFORMAT
NAME "GTiff"
DRIVER GDAL/GTiff
MIMETYPE "image/tiff"
IMAGEMODE RGB
FORMATOPTION "COMPRESS=DEFLATE"
EXTENSION "tif"
END
FONTSET "osmfonts.lst"
MAXSIZE 15000
SIZE 800 800
EXTENT -180 -85 180 85
#
# Hier alle Wunschprojektionen aufzaehlen...
#
WEB
METADATA
"ows_enable_request" "*"
wms_srs "EPSG:4326 EPSG:31466 EPSG:31467 EPSG:31468 EPSG:31469 EPSG:25831 EPSG:25832 EPSG:25833 EPSG:25834 EPSG:900913 EPSG:3857"
labelcache_map_edge_buffer "-10"
wms_title "otm-messtisch"
END
IMAGEPATH "/tmp/"
IMAGEURL "/"
END
PROJECTION
"init=epsg:4326"
END
# Der Blattschnitt, wird zur Ausgabe einzelner Blaetter gebraucht. Das Blatt mit der Nummer "mtbnr"
# wird nur umrahmt, alle anderen werden weiss getoent.
LAYER
TYPE POLYGON
STATUS OFF
NAME "fademtb"
MAXSCALEDENOM 100000000
DATA "/home/www/wms/mtb/mtbschnitt.shp"
PROCESSING "LABEL_NO_CLIP=ON"
PROJECTION
"init=epsg:4326"
END
LABELCACHE off
SIZEUNITS meters
CLASSITEM 'name'
VALIDATION
"mtbnr" "[0-9]+"
END
CLASS
EXPRESSION '%mtbnr%'
MAXSCALEDENOM 100000000
STYLE
OUTLINECOLOR "#000000"
WIDTH 20
MAXWIDTH 1
OPACITY 80
END
END
CLASS
MAXSCALEDENOM 100000000
STYLE
COLOR "#ffffff"
OPACITY 60
END
END
END
# Die Rasterkarte aus dem WMS, ueber den Mapproxy eingebunden
LAYER
NAME "openstreetmap"
TYPE RASTER
STATUS OFF
CONNECTION "http://localhost/mapproxy/service?"
CONNECTIONTYPE WMS
METADATA
"wms_srs" "EPSG:4326 EPSG:31466 EPSG:31467 EPSG:31468 EPSG:31469 EPSG:25831 EPSG:25832 EPSG:25833 EPSG:25834 EPSG:900913 EPSG:3857"
"wms_name" "osm"
"wms_server_version" "1.1.1"
"wms_format" "image/png"
"wms_connectiontimeout" "180"
"wms_exceptions_format" "application/vnd.ogc.se_xml"
END
END
# Die Fussnote mit dem Copyright.
LAYER
STATUS OFF
NAME "CR"
TYPE annotation
TRANSFORM ll
UNITS PIXELS
LABELCACHE off
FEATURE
POINTS
0 0
END
TEXT "Karten: OpenStreetMap.orgXLizenz: cc-by-sa und odbl"
END
CLASS
LABEL
FONT sc
TYPE TRUETYPE
SIZE 8
WRAP "X"
BUFFER 0
POSITION ur
OUTLINEWIDTH 3
COLOR "#000000"
ANTIALIAS TRUE
OUTLINECOLOR "#ffffff"
END
END
END
# Der Massstabsbalken. Die genaue SIZE und "STATUS EMBED" wird als url-Parameter uebergeben.
SCALEBAR
UNITS kilometers
INTERVALS 3
STATUS OFF
POSITION lr
TRANSPARENT ON
STYLE 0
IMAGECOLOR "#EDEBE6"
LABEL
POSITION ul
TYPE TRUETYPE
FONT sc
SIZE 8
COLOR "#010101"
END
SIZE 250 5
COLOR "#ebebeb"
BACKGROUNDCOLOR "#646464"
OUTLINECOLOR "#000000"
END
#################################################################################################
END # of MAP
PHP-Datei
Die php-Datei wird dazu gebraucht, den schrecklich komplizierten Aufruf des Mapservers durch etwas einfacheres zu ersetzen. Statt alle möglichen Layer und Koordinaten anzugeben reicht dann die Angabe von Nummer und Projektion. An dieser Stelle erfolgt auch der Zugriff auf die Postgis-Datenbank. Der wurde vorher das Shapefile mit dem Blattschnitt verfüttert (mit shp2pgsql).
<?php
// Parameter: mtbnr Die Nummer des gewuenschten Blattes
// res Aufloesung in m/px
// proj Projektion (gk4,utm32 oder utm33)
// format Bildformat (png oder gtiff)
// URL-Parameter uebernehmen, srs als Nummer raussuchen, Header ausgeben
//
$mtbnr=$_GET["mtb"];$res=$_GET["res"];$format=$_GET["format"];
$proj=$_GET["proj"];
if($proj=="gk4"){$srs=31468;}
if($proj=="utm32"){$srs=25832;}
if($proj=="utm33"){$srs=25833;}
if($format=="png") {header("Content-Type: image/png; charset=utf-8");}
if($format=="gtiff") {header("Content-Type: image/tiff; charset=utf-8");}
// Datenbank nach der Bounding Box des Blattes im gesuchten SRS fragen und diese in x1,y1,x2,y2 stecken
$db=pg_pconnect('host=localhost port=5432 dbname=osm user=osm password=osm') or die("Keine DB, sonst alles in Ordnung");
$query="select st_box2D(st_transform(the_geom,$srs)) as bbox from mtb where name='$mtbnr';";
$re=pg_query($db,$query);
$r=pg_fetch_array($re);
$bbox=$r['bbox'];
$bbox=str_replace("BOX(","",$bbox);
$bbox=str_replace(")","",$bbox);
$bbox=str_replace(","," ",$bbox);
list($x1,$y1,$x2,$y2)=explode(" ",$bbox);
// Breite und hoehe des Bildes ausrechnen und in bw und bh stecken
$dimx=$x2-$x1;$dimy=$y2-$y1;
$bw=round($dimx/$res);$bh=round($dimy/$res);
// Nochmal 26 Pixel in der Hoehe des Bildes dazunehmen, wegen Schriftzug und Massstab unten
// Dann nochmal die neue Hoehe in der echten Welt ausrechnen
$y1-=26*$res;
$dimy=$y2-$y1;
$bh=round($dimy/$res);
// Die Breite des Massstabsbalkens
$sbx=round($bw/3+25);
// der feste Teil der URL des WMS
$url="http://localhost/cgi-bin2/mapserv6?map=..%2Fmaps%2Fmesstischblatt.map&LAYERS=openstreetmap,fademtb,cr&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&SRS=EPSG%3A$srs";
// Die Variablen Teile: BBOX, FORMAT, Blattnummer
// Bei kleineren Blaettern wird noch der Masstab zugeschaltet (STATUS+EMBED)
// und seine groesse festgelegt (SIZE)
$url=$url."&BBOX=$x1,$y1,$x2,$y2"."&WIDTH=$bw"."&HEIGHT=$bh";
$url=$url."&FORMAT=$format";
$url=$url."&MTBNR=$mtbnr";
if(($res>=4)){$url=$url."&map.scalebar=STATUS+EMBED+SIZE+$sbx+5";}
// Die URL wird an curl uebergeben und ausgegeben
$curl = curl_init($url);
curl_setopt($curl, CURLOPT_RETURNTRANSFER, FALSE);
curl_setopt($curl, CURLOPT_HEADER, FALSE);
curl_exec($curl);
exit(0);
?>
geoTIFF mit 256 Farben
Für die Pilzsammler war aus Platzgründen ein besonderes Format gewünscht: Die Bilder sollten zwar als geoTIFF vorliegen, allerdings mit 256 Farben. Soweit ich weiss, kann Mapserver das nicht erzeugen, sondern nimmt bei geoTIFF immer 24 Bit/px. Bei PNG sind 8-Bit/pixel kein Problem (einfach "png8" als Ausgabeformat nehmen), aber der gdal-Bibliothek, mit der die TIFFs erzeugt werden kann man diesen Wunsch nicht mitteilen. Solche Formate müssen also ausserhalb des Mapservers erzeugt werden.
- Man lässt sich ein 8-bit-PNG erzeugen
- lässt sich das Wordfile dazu ausgeben. Dazu schreibt man in das php-Script von oben statt der curl-Ausgabe:
echo ($x2-$x1)/$bw;echo "\n"; echo "0\n"; echo "0\n"; echo (-1)*($y2-$y1)/$bh;echo "\n"; echo $x1+($x2-$x1)/$bw/2;echo "\n"; echo $y2-($y2-$y1)/$bh/2;echo "\n";
- und verwandelt dann dieses karte.png und dem dazugehörigen worldfile karte.pgw in ein tiff:
/usr/bin/gdal_translate -co "COMPRESS=DEFLATE" -co "TFW=YES" karte.png karte.tiff
gdal_translate lässt in diesem Fall das Bild mit seiner Farbtiefe und Farbpalette in Ruhe und erzeugt ein Tiff mit 8 Bit pro Bildpunkt und einer Palette mit 256 Farben.
Zum Ansehen
Eine Seite zum bequemen Ansehen und Exportieren der bayerischen Messtischblätter aus der OpenTopoMap oder dem "osm.org Mapnik-Stil" gibt es auf geo.dianacht.de/mtb/. Dort ist auch das Shapefile mit dem Blattschnitt verlinkt und eine Liste der schon vorgerenderten Blätter in einer bestimmten Auflösung. Ausserbayerische Blätter können dort theoretisch auch erzeugt werden. Da diese aber nicht aus dem Cache kommen, ist der Zugriff stark gedrosselt. Nach vier oder fünf Blättern dürfte es keinen Spass mehr machen...