User:Maxbe/Messtischblätter aus OSM-Tiles

From OpenStreetMap Wiki
Jump to: navigation, search

Worum gehts?

Maxbe otm 6332 15 gk4.png

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

Maxbe mtbnummern.png

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

Maxbe mtbschiefutm32.png

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

Maxbe mtbvergleich.png

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