OSM: Zugriff auf die Attribute des Admin_Centre beim Rendern von Verwaltungsgrenzen

In der OpenStreetMap-Datenbank sind die Informationen über Verwaltungseinheiten (Gemeinden, Kreise, Länder) meist auf die Grenzrelation und den Knoten mit der Rolle admin_centre verteilt. Will man zum Rendern der Grenze oder Fläche einer Verwaltungseinheit auf Informationen zurückgreifen, die beim Knoten gespeichert sind, z.B. das Attribut population oder name:xx, dann stellt man fest, dass üblichen Werkzeuge diese Informationen nicht bereit stellen. Nachfolgend werden Methoden beschrieben, wie man dieses Problem bei mit osm2pgsql erstellten Datenbanken und bei aus OSM erzeugten Shape-Dateien lösen kann. PostGis / osm2pgsql osm2pgsql überträgt zwar Attribute der äußeren Wege, also der Wege mit Rolle outer, auf das für die Relation erzeugte Polygon, für den Knoten mit der Rolle admin_centre wird nichts vergleichbares gemacht. Es wird auch kein Link vom Polygon zum Knoten erzeugt. Die im Slim-Modus von osm2pgsql erzeugte Tabelle <prefix>_rels enthält zwar in der Spalte members sämtliche Mitglieder der Relation und deren OSM-ID, aber leider in einer nicht direkt und effizient mit SQL auswertbaren Form. Beispiel: psql -d osm -c "select members from planet_osm_rels\ where hstore(tags)->'name'='Corneilla-la-Rivière';" members ---------------------------------------------------------------------------- {w121772197,outer,w100006391,outer,w121772204,outer,n287557758,admin_centre} Die Member der Relation werden also in einem Array abgelegt, zwei aufeinander folgende Array-Elemente pro Member. Mit einer kleinen Funktion kann man die Node-Id des Admin_Centre extrahieren: CREATE OR REPLACE FUNCTION admin_centre(anyarray) RETURNS bigint LANGUAGE sql AS $function$ SELECT cast(substring($1[i] from 2) as bigint) FROM generate_series(array_lower($1,1),array_upper($1,1),2) i WHERE $1[i+1]='admin_centre';$function$ Damit kann man nun die Polygone der Verwaltungsgrenzen zusammen mit den Attributen des Admin_Centre abfragen: select poly.name, point.name, point."name:ca", point.population from planet_osm_rels r, planet_osm_polygon poly, planet_osm_point point where poly.admin_level='8' and admin_centre(r.members) = point.osm_id and r.id = -poly.osm_id ; Je nach Anwendungsfall und Performanceanforderung kann man damit auch die benötigten Attribute in die jeweils andere Tabelle übertragen oder eine Tabelle mit der Zuordnung von Knoten zu Relationen erstellen: create table admin_rel as select -id as rel_id, admin_centre(members) as admin_centre from planet_osm_rels where hstore(tags)->'boundary'='administrative' and hstore(tags)->'admin_level' in ('8','6','4'); Shape-Dateien Hierfür wird die Osmium-Bibliothek, das darauf basierende Javascript-Tool osmjs und QGis benutzt. Mit osmjs können auf einfache Weise Shape-Dateien aus den Objekten einer OSM-Datei zu erzeugt werden. Wie bei osm2pgsql werden dabei zwar Polygone aus Grenzrelationen erstellt, aber ohne Bezug zum admin_centre. Es ist jedoch möglich, eine CSV-Datei zu erstellen, welche diesen Bezug enthält: /* Osmium Javascript Example: shape_export.js run with: osmjs -2 -m -l sparsetable -j shape_export.js OSMFILE */ var shp_places = Osmium.Output.Shapefile.open('./places', 'point'); shp_places.add_field('id', 'string', 12); shp_places.add_field('type', 'string', 32); shp_places.add_field('name', 'string', 32); var shp_boundaries = Osmium.Output.Shapefile.open('./boundaries', 'polygon'); shp_boundaries.add_field('id', 'integer', 10); shp_boundaries.add_field('type', 'string', 32); shp_boundaries.add_field('boundary', 'string', 255); shp_boundaries.add_field('admin_level', 'string', 255); var csv_file = Osmium.Output.CSV.open("./rels.csv"); var node_tags = { place: { village: 'villaget', city: 'city', town: 'town'} //, //shop: { supermarket: 'supermarket' } } Osmium.Callbacks.init = function() { print("Init"); } Osmium.Callbacks.node = function() { for (var key in this.tags) { if (node_tags[key]) { var type = node_tags[key][this.tags[key]]; if (type) { shp_places.add(this.geom, { id: this.id, type: type, name: this.tags.name }); } } } } Osmium.Callbacks.area = function() { if (this.tags.boundary) { shp_boundaries.add(this.geom, { id: this.id, type: this.tags.boundary, admin_level: this.tags.admin_level}); } } Osmium.Callbacks.relation = function() { for (var key in this.members) { for (var key2 in this.members[key]) { if (key2=='role' && this.members[key].role=='admin_centre'){ csv_file.print(this.id*2+1, this.members[key].ref, this.members[key].role, this.members[key].type); } } } } Osmium.Callbacks.end = function() { shp_places.close(); shp_boundaries.close(); csv_file.close(); print("Done"); } Das Skript erstellt eine Shape-Datei mit den Grenzpolygonen, eine Shape-Datei mit den Place-Knoten und eine CSV-Datei, welche die Ids der Grenzpolygone und der zugehörigen Place-Knoten enthält. Dabei wird berücksichtigt, dass osmjs den Polygonen eine Pseudo-ID mit dem Wert [2 * OSM-ID der Relation + 1 zuweist. Diese drei Dateien lädt man als Layer in QGis. Dann verknüpft man den CSV-Daten-Layer über die Knoten-ID mit dem Knoten-Layer und anschliessend den Polygon-Layer über die Pseudo-ID mit dem Knoten-Layer. Will man die Daten mit Tilemill / Mapnik verenden, kann man die Polygon-Layer als Shape-Datei speichern.

4. Oktober 2013

Version V027 des Skripts rel2gpx.pl

In diese Version sind zwei Patches von Andreas Tille eingearbeitet. Zur Download-Seite

30. September 2013

QGIS mit Python-Unterstützung unter Debian Jessie

Derzeit funktioniert die Python-Einbindung von QGis in Debian Jessie/Testing nicht. Das liegt daran, dass QGIS 1.7 die Version 8.1 des SIP-Interfaces unterstützt, in Jessie aber bereits auf SIP Version 10 umgestellt wurde. Versuche mit den Debian-Paketen für QGIS 1.8 von https://qgis.org/debian/ haben das Problem nicht behoben. Letztendlich habe ich QGIS 2.0 aus den Quellen erzeugt. Quellen von qgis.org herunterladen. Quellarchiv in ein Verzeichnis entpacken Compilieren entsprechend der Anleitung in der Datei INSTALL. Bei mir waren nur die Schritte aus Abschnitt 3.3. und 3.7. notwendig. Beim Konfigurieren mit ccmake wird der Pfad PYTHON_LIBRARY nicht automatisch gesetzt und muss daher von Hand eingetragen werden (in ccmake), z.B. /usr/lib/x86_64-linux-gnu/libpython2.7.so ldconfig ausführen. (Stand: 14.09.2013)

9. Juli 2013

GPX-Track auf OpenStreetMap-Karte anzeigen und ausdrucken

Es gibt eine Reihe von Tools und Online-Diensten, mit denen man GPX-Tracks auf einer OpenStreetMap-Karte anzeigen kann: QLandkarteGT, GpsPrune, gpsies.com,... Beim Ausdrucken bzw. Erzeugen einer Vektorgrafik stößt man allerdings schnell an Grenzen. Die Darstellung der Tracks (Farbe, Breite, Transparenz) ist nicht konfigurierbar und der Zoomlevel wird anhand des Kartenausschnitts gewählt, usw. Dieses Video hat mir einen Weg zur Lösung aufgezeigt: Quantum GIS (QGIS) mit OpenLayers-Plugin. OpenLayers-Plugin Installation Das OpenLayers-Plugin ermöglicht es, OSM-Kacheln und Kacheln aus anderen mit OpenLayers darstellbaren Quellen (Google, Bing) in einem QGIS-Layer anzuzeigen. Dazu muss zunächst das Plugin installiert werden. Bis QGIS 1.8 funktioniert das problemlos mit mit dem Plugin Installer. Zunächst muss im Menü "Erweiterungen" der Plugin Installer installiert werden, dann kann mit diesem das OpenLayers-Plugin installiert werden. QGIS 1.9/2.0 ist mit dem Plugin nicht kompatibel. Es existiert aber ein Workaround, der hier beschrieben ist: [node:254]. Layer hinzufügen Danach erscheint im Menü "Erweiterungen" ein Untermenu "Openlayers Plugin", über das eine Kachelquelle ausgewählt und mit einem Layer verknüpft werden kann. Die Kacheln werden entsprechend dem Zoomlevel geschrumpft bzw, gedehnt, was zu einer unscharfen Darstellung führt. Ein Workaround besteht darin, den Maßstab explizit so einzustellen, dass die Kacheln in Original-Auflösung angezeigt werden. Stellt man den Maßstab 1:2257 ein, erhält man die Kacheln des Zoomlevels 18. Durch Zoomen mit der Maus bekommt man dann die Kacheln in den OSM-Zoomleveln angezeigt. Weitere Kachelquelle hinzufügen Eine zusätzliche Kachequelle kann auf einfache Weise Hinzugefügt werden. Im Verzeichnis ~/.qgis2/python/plugins/openlayers_plugin/html existiert für jede Kachelquelle eine HTML-Datei, die den Javascript-Code enthält, mit dem der entsprechende OpenLayers-Layer erzeugt wird. Man erstellt nun eine Kopie einer dieser Dateien und passt den Code entsprechend an. Für OSM-kompatible Kachelquellen genügt es in der Regel, die Datei osm.html zu kopieren und den URL anzupassen. Anschließend muss noch in ~/.qgis2/python/plugins/openlayers_plugin/openlayers_plugin.py, etwa bei Zeile 120, ein Eintrag für die neue Kachelquelle nach diesem Muster eingefügt werden: self.olLayerTypeRegistry.add( OlLayerType(self, 'OpenStreetMap', 'osm_icon.png', 'osm.html', True) ) nach einem Neustart von QGIS ist die neue Kachelquelle verfügbar. Alternative Methode Ein anderer Weg zur Anzeige der OSM-Kacheln in QGIS ist hier beschrieben. GPX-Track Ein GPX-Track kann als Vektor-Layer angezeigt werden: Menüpunkt Vektorlayer hinzufügen anwählen im Dialog die GPX-Datei auswählen und öffen im anschliessenden Auswahlfenster den Layer tracks und ggf. weitere Layer auswählen und mit OK bestätigen im Layer-Fenster die Eigenschaften der Track-Linie festlegen: Farbe, Breite,... Karte drucken / Grafik erzeugen Hierzu wird der Dialog Neue Druckzusammenstellung geöffnet und mit neue Karte hinzufügen ein Kartenauschnit in das Druckblatt eingefügt. Dabei wird der Kartenausschnitt aus dem Hauptfenster übernommen. Dieser kann im Modus Inhalt verschieben nach Belieben gezoomt und verschoben werden Die so erzeugten Ausdrucke haben eine miserable Qualität. Der Zoomlevel der Kacheln wird entsprechend Auflösung des Ausdrucks und der Größe des Kartenauschnitts automatisch gewählt. Dann werden die Kacheln skaliert um den gewählten Kartenausschnitt zu erhalten. Das Ergebnis ist in den meisten Fällen ziemlich unscharf und für den Druck in der Regel nicht verwendbar. Durch geeignete Wahl der Druckgröße und des Maßstabs des Kartenausschnitts kann man trotzdem eine akzeptable Bildqualität erreichen. Beispiele: Format Auflösung (dpi) Maßstab Zoomlevel A43005725015 A41505725014 A4725725013 A17257.25013 A172144.44812 A1300144.44814 A13002.311.16810 A172577.79210 Das beste Ergebnis erziele ich allerdings mit einer Bildschirmkopie, da die Kacheln auf dem Bildschirm unskaliert dargestellt werden.

9. Juli 2013

Mapping Party in Collioure

Am 1. Juni 2013 organisieren Annexe21, Perpinux und die OpenStreetMap-Mapper des Departements Pyrénées-Orientales eine Mapping Party in à Collioure. Alles weitere auf https://annexe21.net/cartopartie2013

28. Mai 2013

Echoes of Robert Johnson, Prades 02/03/2013

4. März 2013

Big Daddy Wilson au Mediator de Perpignan (17/02/2013)

18. Februar 2013

PostGIS: Wege ausserhalb von Landuse-Polygonen ermitteln

Für meine Velocarte wollte ich einen Overlay-Layer erstellen, in dem Straßen, die sich wegen ihres geringen Verkehrsaufkommens besonders für Alltagsradler eignen, farblich hervorgehoben werden. In einem ersten Schritt habe ich dafür alle Unclassified- und Service-Wege selektiert. Es hat sich aber gezeigt, dass darunter auch viele stark befahrene Straßen im innerstädtischen Bereich und in Industrie- und Gewerbegebieten waren. Um diese Straßen auszuschließen, habe ich zunächst die PostGIS-Funktion ST_WITHIN verwendet: select osm_id,name,ref,highway,tracktype,way from planet_osm_line as l where highway in ('unclassified','service') and not exists (select * from planet_osm_polygon as p where p.landuse in ('residential','industrial','commercial') and st_isvalid(p.way) and st_within(l.way,p.way)); Da ST_Within auf vollständiges Enthaltensein prüft, werden von dieser Abfrage auch Straßen erfasst, die aus dem Umland in ein Wohn- oder Gewerbegebiet hineinführen. An dieser Stelle bin ich auf die Funktion ST_Difference gestoßen, welche, angewendet auf einen Weg und ein Polygon, den Teil des Wegs liefert, der nicht im Polygon liegt. Der erste Ansatz war: ...

2. Februar 2013

PostGIS: Daten für den Renderingprozess vorbereiten

Die nachfolgende Anleitung gilt für Debian Wheezy (Stand 01/2013) Anlegen der PostGIS-Datenbank Installation der PostgreSQL-Pakete: aptitude install postgresql-9.1 postgresql-9.1-postgis postgresql-contrib-9.1 postgis Installation des Pakets postgis. Wheezy enthält PostGIS 1.5. Wenn, wie in meinem Fall, Funktionen benötigt werden, die erst ab PostGIS 2.0 verfügbar sind, kann diese Version problemlos aus den Quellen installiert werden. Anlegen eines Users und einer Datenbank. Als Superuser: su - postgres createuser <username> Soll die neue Rolle ein Superuser sein? (j/n) j #psql --command "ALTER USER <username> WITH SUPERUSER"; psql --command "ALTER USER <username> WITH ENCRYPTED PASSWORD 'osm'"; Als User <username>: createdb osm psql -d osm --command "CREATE EXTENSION hstore;" #createlang plpgsql osm psql -d osm -f /usr/share/postgresql/9.1/contrib/postgis-1.5/postgis.sql psql -d osm -f /usr/share/postgresql/9.1/contrib/postgis-1.5/spatial_ref_sys.sql Importieren der OpenStreetMap-Daten mit osm2pgsql Mit Osm2pgsql können OpenStreetMap-Daten in eine PostGIS-Datenbank importiert werden. Im Wesentlichen werden dabei die OSM-Objekte (Ways, Nodes, Relations), ihre Attribute (Tags) und Geometrie in Datenbanktabellen übernommen. Nodes werden nach planet_osm_nodes, Ways nach planet_osm_line und geschlossene Ways nach planet_osm_polygon importiert. Welche Attribute übernommen werden sollen, wird in einer Konfigurationsdatei (Import Style) festgelegt. Relationen vom Typ Multipolygon, Boundary und Route werden nach planet_osm_rels importiert. In diese Tabelle werden nur die Attribute der Relation und Verweise auf ihre Members eingetragen. Über die Konfigurationsdatei kann festgelegt werden, dass für Relationen mit bestimmten Attributen auch die Geometrie nach planet_osm_polygon bzw. planet_osm_line importiert werden. Daten importieren: osm2pgsql -c -m -s -d osm -U <username> -W -H localhost --bbox <bbox> <osmdatei> <osmdatei> ist die Datei OSM-Format (.osm, .osm.bz2 ou .pbf), welche die zu importierende Daten enthält. Nachbearbeiten der importierten Daten Die mit osm2pgsql und dem Default-Importstyle erzeugte Datenbank enthält alle Daten, die man für das Rendern einer Karte im OSM-Mapnik-Standard benötigt. Für einen Karte in einem davon wesentlich abweichenden Stil kann es sinnvoll bzw. notwendig sein, die importierten Daten mit SQL nachzubearbeiten. Ein einfaches Anwendungsbeispiel ist das Ersetzen der Werte "yes", "true", "1" durch "yes", um Abfragen in Mapnik oder Tilemill zu vereinfachen: UPDATE planet_osm_line SET tunnel = (CASE WHEN tunnel IN ('yes','true','1') THEN 'yes'::text ELSE tunnel::TEXT END); Dies kann natürlich auch in der Datenbankabfrage in Mapnik/Tilemill gemacht werden. Erledigt man das bereits im Vorfeld, dann werden die Abfragen in Mapnik/Tilemill kompakter und lesbarer. Komplexere Anwendungsfälle werden in einem separaten Beitrag behandelt.

30. Januar 2013

Tilemill unter Debian Wheezy installieren

node.js & npm Um Tilemill übersetzen zu können, muss erst node.js installiert werden. Quellen von https://nodejs.org/download/ herunterladen Übersetzen und installieren: tar xvzf tar xvzf node-v<version>.tar.gz cd node-v<version> ./configure make su -c 'make install' Tilemill Tilemill übersetzen: git clone https://github.com/mapbox/tilemill.git cd tilemill npm install Anschließend kann Tilemill mit ./node_modules/tilemill/index.js gestartet werden.

30. Januar 2013