OSM Kacheln als WMS mit MapServer: Wie Layer umschalten?

holgermappt · 2011-02-13T17:01:58.000Z

Ich möchte OSM Kacheln mittels MapServer (mapserver.org) als WMS bereitstellen. Zum Testen nehme ich JOSM. Dazu generiere ich mir Worldfiles für die Kacheln (mit lat/lon Koordinaten in Dezimalgrad), mache daraus mit gdaltindex Shapefiles, und binde die in die Mapserver-Map-Datei wms.map ein. Da die Kacheln für verschiedene Zoomstufen vorliegen, soll je nach angefragter Auflösung umgeschalten werden. Dafür gibt es eine Gruppe mit je einem Layer pro Zoomstufe.

Der erste Stolperstein ist die Angabe der Projektion der Kacheln. Da die in Merkator-Projektion sind, würde ich ja eigentlich EPSG:3857 angeben. Da JOSM aber EPSG:4326 anfragt, kommt immer nur ein weißes Bild zurück, wenn ich nicht EPSG:4326 als Quellprojektion angebe. Weiß ist der Default-Hintergrund. Geht das nicht anders oder habe ich hier etwas übersehen?

Aber jetzt zur eigentlichen Frage: Wie genau schalt ich zwischen den Zoomstufen um? Da der Zoom mit der Auflösung korreliert, sollte MINSCALEDENOM/MAXSCALEDENOM das Richtige sein. Aber welchen Wert soll ich nehmen? Wie hängt der Denominator mit der Auflösung und somit der Zoomstufe zusammen? Mit Auflösung meine ich Meter pro Pixel. Die Auflösung ändert sich auch noch vom Äquator zu den Polen hin.

Hat wer so etwas schon einmal realisiert und kann mir ein Beispiel-Mapfile zeigen? Meine Version sieht so aus:


MAP
  NAME "DEMO"
  WEB
    METADATA
      "wms_title"           "WMS Demo Server"
      "wms_onlineresource"  "http://localhost/cgi-bin/mapserv?map=wms.map&"
      "wms_srs"             "EPSG:4326 EPSG:3857"
    END
  END
  PROJECTION
    "init=epsg:4326"
  END
  IMAGECOLOR 255 255 255
  TRANSPARENT ON

  LAYER
    NAME "OSMTiles15"
    GROUP "OSMTiles"
    PROJECTION
      "init=epsg:4326"
#      "init=epsg:3395"
#      "init=epsg:3857"
#      "init=epsg:900913"
    END
#    PROJECTION
#      "proj=merc"
#      "lon_0=0"
#      "k=1"
#      "x_0=0"
#      "y_0=0"
#      "ellps=WGS84"
#      "datum=WGS84"
#      "units=m"
#      "no_defs"
#    END
    METADATA
      "wms_title"           "OSM Tiles Zoom 15"
    END
    STATUS ON
    TYPE RASTER
    TILEINDEX "osmtiles15.shp"
    TILEITEM "Location"
    MINSCALEDENOM ???
  END

  LAYER
    NAME "OSMTiles16"
    GROUP "OSMTiles"
    PROJECTION
      "init=epsg:4326"
    END
    METADATA
      "wms_title"           "OSM Tiles Zoom 16"
    END
    STATUS ON
    TYPE RASTER
    TILEINDEX "osmtiles16.shp"
    TILEITEM "Location"
    MAXSCALEDENOM ???
    #MAXGEOWIDTH ???
  END

END

Danke,
Holger

SunCobalt · 2011-02-13T17:47:04.000Z

Hallo Holger,

mir ist nicht ganz klar was Du da machst. Normalerweise würde ich die Kacheln verschmelzen und dann georeferenzieren. Dabei gibt man ja den EPSg Code vor. Einfacher ist es wahrscheinlich, sich mit Maperitive eine Grafik der Gegend mit einem hohen Zoom zu erstellen. Da ist dann netterweise gleich die World-File dabei. Die sollte man dann noch etwas bearbeiten und den extent (siehe unten) auslesen. Danach sieht die mapfile dann ungefähr so aus.
Grundsätzlich gibst Du im Layer den mit PROJECTION EPSG Codes des vorliegenden Materials an und im Kopf der Mapfile und bei METADATA den Code, der bspw von JOSM angefragt wird. gdal sollte das dann automatisch umsetzen. Um ein shapefile in ein sichtbares Bild umzuwandeln,brauchst Du shp2img und musst regeln festlegen, wie bestimmte Dinge in der Shapefile auszusehen haben. Siehe CLASS im Beispiel unten

mapfile header


MAP
    NAME OSM-WMS
    FONTSET "/osm/mapserver/fonts.list"
    STATUS ON
    UNITS METERS


    WEB         
        METADATA
            "wms_title"     "Thomas WMS-Server"
            "labelcache_map_edge_buffer" "-10"
            "wms_onlineresource" "to_be_replaced_by_wms.py"
            "wms_srs"       "EPSG:4326 EPSG:900913"
            "wms_connectiontimeout" "90"
        END
    IMAGEPATH    "/osm/mapserver/imagepath/"

    END
    CONFIG "PROJ_LIB" "/usr/share/proj"
    CONFIG "MS_ERRORFILE" "stderr"
    EXTENT 0 0 180 180
    MAXSIZE 4000
    SIZE 800 800
    UNITS dd
    IMAGECOLOR "#B3C6D4" 
    PROJECTION
        "init=epsg:4326"
    END

Layer eines georeferenzierten Bildes


 LAYER
    NAME "Hiddensee"
    STATUS ON
    TYPE RASTER

    DATA "/osm/mapserver/hiddensee/Hiddensee_tiled.tif"

    PROJECTION
          "init=epsg:4326"
    END

    METADATA
     "wms_title"  "Hiddensee"
     "wms_extent" "13.065000 54.504001 13.132001 54.605000"
     "wms_srs"    "EPSG:4326"
     "copyright"    ""
    END
  END

Beispiel einer gerenderten Shapefile


  LAYER
    NAME "RoadsVmap0"
    TYPE LINE
    METADATA
    "wms_title" "RoadsVmap0"        
    "wms_srs"    "EPSG:4326"  #copied from other layer
    "wms_extend" "17.071282 41.521486 24.781496 46.529345"
    END
    STATUS ON
    TRANSPARENCY 100
    PROJECTION
    "init=epsg:4326" 
    END
    CLASS
       NAME "trans-road-I" 
       STYLE
         WIDTH 2        
         COLOR 37 213 202     
       END
    END
  END

Den extent bekomme ich so “python print_wms_extent.py meinbild.tif”

Das Bild in Tiles aufteilen “gdal_translate -of GTiff -co “TILED=YES” vorher.tif nachher.tif”

Und schon vorab verschiedene Zoomlevel bereitstellen “gdaladdo -ro --config USE_RRD YES -r gauss meinbild.tif 2 4 8 16 32 64 128 256”

MINSCALEDENOM/MAXSCALEDENOM sind afaik dafür, wenn Du Rendering Regeln brauchst, also wenn der Mapserver als Renderer arbeitet wenn Du abhängig vom Zoomlevel unterschiedliche Daten/Shapefiles nehmen willst. Das sieht dann bspw so aus


LAYER
    STATUS ON 
    GROUP "default" 
    PROJECTION "init=epsg:4326" 
    END 
    CONNECTIONTYPE postgis CONNECTION "host=localhost dbname=osm user=www-data password=www-data port=5432" 
    PROCESSING "CLOSE_CONNECTION=DEFER" 
    maxscaledenom 25010 
    minscaledenom 10010
    TYPE LINE
    NAME "roadsclose"
    DATA "way from (select way, osm_id ,highway,ref, namesr, tunnel from osm_line where highway is not null order by priority desc) as foo using unique osm_id using srid=4326"
    LABELITEM "namesr"
    CLASSITEM "highway"
    CLASS
        EXPRESSION "motorway"
        TEXT ([ref])
        STYLE
            WIDTH 10
            OUTLINEWIDTH 1
            OUTLINECOLOR "#24277A"
        END
        STYLE
            WIDTH 10
            COLOR "#4F719D"
        END
        LABEL
            TYPE TRUETYPE FONT sc PARTIALS FALSE OUTLINEWIDTH 3
            PRIORITY 1
            FONT scb
            SIZE 9
            MINDISTANCE 150
            MAXLENGTH 6
            COLOR 255 255 255
            BACKGROUNDCOLOR "#FF5C1D"
            BACKGROUNDSHADOWCOLOR "#FDB462"
        END
    END

holgermappt · 2011-02-13T19:58:33.000Z

Mein Ausgangspunkt sind normale Kacheln, die von tile.openstreetmap.org kommen. Da ich manchmal keinen Onlinezugang habe, speichere ich die Kacheln und nutze sie über einen lokal laufenden MapServer als Hintergrund in JOSM. Damit kann man schon mal gut die Tracks anschauen, auch wenn man erst später etwas an den OSM-Daten machen kann. Dabei stapele ich die Kacheln im Shapefile übereinander, wobei die mit der höchsten Auflösung ganz oben liegen. Wenn also Kacheln in hoher Auflösung fehlen, dann werden dort entsprechend gröbere Bildbereiche vom MapServer eingebaut, weil durch fehlende Kacheln in hoher Auflösung die darunter liegenden sichtbar werden. Das Shapefile wird nicht direkt angezeigt, sondern dient als Information, welche Kacheln an einer bestimmten Stelle liegen. Das funktioniert so auch ganz gut. Problem ist nur, dass es zu viel Arbeit für den MapServer gibt, wenn zu viele Kacheln in die Berechnung eines angeforderten Bildes einfließen. Daher generiere ich Shapefiles, die nur Kacheln bis zu einer bestimmten Zoomstufe enthalten. Auf die soll dann umgeschaltet werden, wenn ich im JOSM herauszoome. Und um diese Umschalterei geht es. Ich will keine Kacheln aus Rohdaten erzeugen, sondern die bereits vorhandenen anzeigen. Ich will auch keine extra Grafiken in anderen Zoomstufen generieren, weil das der MapServer schon alles macht.

Was mir hauptsächlich fehlt ist ein Verständnis für die Bedeutung von Denominator und die Einbindung von Bildern in Merkator-Projektion. Ich könnte den Denominator auch empirisch ermitteln, hätte aber gerne eine Regel zur Herleitung. Die MapServer Doku verwendet leider nur Begriffe, mit denen ich nichts anfangen kann.

Holger

SunCobalt · 2011-02-13T20:52:35.000Z

MINSCALEDENOM/MAXSCALEDENOM ist der minimale/maximale Maßstab für ein Layer, also 1:5000 bis 1:50000 wäre “MINSCALEDENOM 5000” + “MAXSCALEDENOM 50000”. Wenn Du für jede OpenLayers Zoomstufe den Maßstab rausfindest, hast Du die Daten schnell. Du darfst nur keine Lücke lassen, sprich ein WMS Layer (minscale…) 5000 bis (maxscale…) 9999 und das nächste dann ab (minscale) 10000, da JOSm auch zwischenrein zoomen kann. Für jedes Layer müsstest Du beides angeben, nur für das kleinste Layer nur MAXSCALEDENOM und für das größte Layer nur MINSCALEDENOM. Wie das mit dem Stapeln funktionieren soll, weiß ich gerade nicht. Evtl von oben immer mit der schlechtesten Zoomstufe anfangen und MAXSCALEDENOM weglassen. Dann klatscht er vielleicht solange Tiles übereinander, bis er bei höheren Zoomstufen nichts mehr findet. Aber ob Du bei einer Deutschlandansicht wirklich Kacheln mit einem Zoomlevel 18 sehen möchtest, weiß ich nicht. Daher vielleicht doch ein großzügiges MAXSCALEDENOM setzen.
Ich finde trotzdem, dass Du Dir das viel zu kompliziert machst.

Die Antwort auf die EInbindung von Bilder in 900913 hatte ich Dir oben schon gegeben

SunCobalt · 2011-02-13T21:31:49.000Z

achja, gdal hast Du installiert? Du brauchst gdalwarp um 900913 in 4326 zu übersetzen. Es muss auch für den Mapserver auffindbar sein. Tippe mal irgendwo gdalinfo in die Kommandozeile und schaue, was passiert. Und in der (ebenfalls auffindbaren) epsg Datei muss ein Eintrag für 900913 existieren, der so aussieht “<900913> +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 +no_defs <>” …ohne die “”

Dann immer "“wms_srs” “EPSG:4326” unter METADATA im mapfile Header und auch im Layer (beim zweiten weiß ich nicht obs nötig ist, schadet aber nicht) und


PROJECTION
          "init=epsg:900913"
END

in die Layers, deren Daten in der Merkator Projektion vorliegen.

holgermappt · 2011-02-14T17:03:11.000Z

SunCobalt:

MINSCALEDENOM/MAXSCALEDENOM ist der minimale/maximale Maßstab für ein Layer, also 1:5000 bis 1:50000 wäre “MINSCALEDENOM 5000” + “MAXSCALEDENOM 50000”. Wenn Du für jede OpenLayers Zoomstufe den Maßstab rausfindest, hast Du die Daten schnell.

Genau da komme ich nicht weiter. Was besagt ein Maßstab von 1:5000? Auf einer Papierkarte würde es bedeuten, dass 1 cm in der Karte 50 m in der Wirklichkeit entspricht. Aber wir haben hier Bilder mit Pixeln. 1 Pixel im Bild entspricht 5000 Pixeln in der Wirklichkeit? Wie viele Pixel hat dann die Welt? Wie viele Pixel hat der Äquator? Wie viele sind entlang des 40. Breitengrades? Wie viel “ganz oben” (85,051 °N)? Welchen Maßstab hat Zoom 15?

gdal ist installiert, gdalwarp auf der Kommandozeile verfügbar. Da libgdal.so.1 als Abhängigkeit des MapServers in der Paketverwaltung angezeigt wird, sollte auch die Umprojektion funktionieren. Hat wer ein Worldfile für eine der Merkator-Projektionen als Beispiel? Oder Koordinaten aus einen georeferenzierten Bild? Ich glaube fast die Koordinaten in meinen Worldfiles sind irgendwie verkehrt. Oder kann ich selbst Bild+Worldfile+EPSG-Kode unter Angabe eines anderen EPSG-Kodes in Bild+Worldfile umwandeln?

holgermappt · 2011-03-06T21:11:10.000Z

Wie hängt der Denominator mit der Auflösung und somit der Zoomstufe zusammen?

Den entscheidenden Hinweis habe ich auf http://www.geog.fu-berlin.de/de/Karto/umn_karten/legenden.shtml gefunden: “Der UMN MapServer geht bei der Berechnung des numerischen Maßstabs von einer Bildschirmauflösung von 72 ppi (pixel per inch) aus.”. Damit ergibt sich, dass Zoom 15 am Äquator einen Maßstab von 1:13542 hat, bei 85,051 °N ist er 1:1168. Wenn die Auflösung im Vordergrund steht, muss also der MINSCALEDENOM für Zoom 15 auf 13542 gesetzt werden. Das hat natürlich den Nachteil, dass nahe der Pole ein Layer mit viel zu hoher Auflösung selektiert und die Kacheln runterskaliert werden. Der MAXSCALEDENOM ist entsprechend des nächsten Layers zu wählen. Die 72 ppi ist der Default von RESOLUTION innerhalb von MAP.

Die Projektion habe ich auch hinbekommen. Der richtige Code ist EPSG:3857, die Koordinaten müssen in Metern angegeben. Die Rechenvorschrift dazu habe ich auf http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/ gefunden. EPSG:3857 ist das Gleiche wie EPSG:900913, wobei 900913 für google steht. EPSG:900913 ist ein mehr oder weniger inoffizieller Code der Projektion, ist aber oft zu finden.

Zu erwähnen ist noch, dass die Eigenschaft von JOSM, Merkator (EPSG:3857) Bilder als WGS84 (EPSG:4326) anzufordern, durch die dabei stattfindende Umprojektion zu geringen Qualitätsverlusten führt. Vor allem Schrift wird dabei unleserlich.

Vielen Dank an SunCobalt für einige entscheidende Hinweise.