Bei Liip haben wir vor Kurzem eine Anwendung für Architekturstudent*innen fertiggestellt, die Shapefiles als Datengrundlage verwendet, um Informationen zu Baugrundstücken in der Stadt Zürich zu liefern. Die Informationen werden auf einer Karte angezeigt und können von den Benutzer*innen eingesehen werden. Ein ähnliches Projekt wurde von einem studentischen Team, zu dem ich gehörte, an der Hochschule für Technik der FHNW durchgeführt: eine lustige, benutzerfreundliche Webanwendung, mit der Primarschulkinder etwas über ihre Gemeinde lernen können.

Shapefiles tragen nicht nur dazu bei, die User Experience zu verbessern, sondern beinhalten auch Daten, die, wenn man sie miteinander in Zusammenhang bringt, interessante Erkenntnisse liefern. Schauen wir also mal, was ich herausfinde, wenn ich eine kleine App entwickle, um damit ein paar Shapefiles genauer unter die Lupe zu nehmen.
In diesem Blogbeitrag stelle ich einen praktischen Ansatz zur Arbeit mit Shapefiles vor und zeige euch, wie ihr auch selbst damit arbeiten könnt!

Alpenkonvention

Das erste Hindernis bei der Arbeit mit Shapefiles besteht darin, zu verstehen, was Shapefiles eigentlich sind, und sie überhaupt erst zu erwerben.

Shapefiles sind Binärdateien, die Geodaten enthalten. Sie bestehen in der Regel aus Gruppen von mindestens drei verschiedenen Dateien, die alle denselben Namen haben und sich nur durch ihre Dateiendung unterscheiden: .shp, .shx und .dbf. Ein vollständiges Datenset kann jedoch viel mehr Dateien enthalten. Das Datenformat wurde von Esri entwickelt und in den frühen 1990er Jahren eingeführt. Shapefiles enthalten sogenannte Features – eine geometrische Form mit den dazugehörigen Metadaten. Ein Feature könnte ein Punkt (einzelne XY-Koordinaten) und ein Text sein, der diesen Punkt beschreibt. Oder ein Polygon, das aus mehreren XY-Punkten besteht, eine einzelne Linie, eine Multi-Linie usw.

Um an einige Shapefiles für die Arbeit zu kommen, werfe ich einen Blick auf OpenData.swiss. OpenData stellt insgesamt 102 Shapefile-Datenssets (Stand: November 2017) zur Verfügung. Ich muss sie nur herunterladen und kann sofort anfangen, mit ihnen zu arbeiten. Für dieses Beispiel habe ich ein eher kleines Shapefile ausgewählt: Alpenkonvention enthält die Perimeter der Alpenkonvention in der Schweiz.

Da es sich hierbei um Binärdateien handelt, die ich mir nicht in einem Texteditor ansehen kann, brauche ich ein Tool, das mir zeigt, was ich gerade heruntergeladen habe.
QGIS – ein freies Open-Source-Geoinformationssystem
QGIS ist eine umfassende Lösung für alle Arten von Geoinformationssystemen (GIS). Es kann gratis heruntergeladen und installiert werden und verfügt ausserdem über zahlreiche Plugins sowie eine aktive Community. Ich werde es nun nutzen, um einen Blick auf mein zuvor heruntergeladenes Shapefile zu werfen.

Für dieses Beispiel habe ich QGIS mit dem OpenLayers-Plugin installiert. So sieht es aus, nachdem es gestartet wurde:

Um zu sehen, wo ich gerade bin und wo meine Daten hingehören, füge ich OpenStreetMap als Layer hinzu. Ausserdem verschiebe und zoome ich die Karte, damit ich die Schweiz im Blick habe.

Der nächste Schritt besteht darin, das Shapefile zu öffnen, damit QGIS es über der Karte anzeigen kann. Dazu reicht ein Doppelklick auf die .shp-Datei im Dateibrowser auf der linken Seite.

In dieser Ansicht erhalte ich bereits einige Informationen über das verwendete Shapefile, über das Gebiet, das es abdeckt, und über die Anzahl der darin enthaltenen Features. Das Shapefile der Alpenkonvention scheint aus lediglich einem einzigen Feature – einem grossen Polygon – zu bestehen, was ausreichend ist, da es nur den Perimeter der Alpenkonvention in der Schweiz abbildet.

Um zu sehen, welche Bereiche und/oder Einzelheiten es abdeckt, kann ich den Stil des Layers ändern und auf transparent umstellen. Ausserdem zoome ich näher heran, um zu sehen, wie präzis die Form definiert ist. Seine Kanten sollten sich genau mit den Schweizer Grenzen decken.

Wunderbar. Nachdem ich mir die Form des Features angesehen habe, wende ich mich nun den Metadaten des Shapefiles zu. Diese können in der Attributtabelle des Shapefiles im Browser unten links eingesehen werden.

Auch diese Datei enthält nicht viele Metadaten, aber jetzt sehe ich, dass es darin zwei Features gibt. Wie dem auch sein, es gibt wohl noch mehr interessante Shapefiles, mit denen man arbeiten kann. Aber ich werde beim Thema Alpen bleiben.

Von Lawinen und Steinböcken

Beim weiteren Durchstöbern der Shapefiles von OpenData bin ich auf zwei weitere Shapefiles gestossen, die möglicherweise interessante Daten liefern könnten. Die Verbreitung der Steinbockkolonien und der Stand der Naturgefahrenkartierung in den Gemeinden bezüglich Lawinen.

So genial QGIS für die Analyse der vorliegenden Daten auch ist, ich möchte meinen eigenen Explorer erstellen und das vielleicht sogar auf Basis mehrerer Shapefiles und mit meinem eigenen Design. Dazu gibt es verschiedene Bibliotheken in den unterschiedlichsten Sprachen. Im Bereich Webentwicklung wären die drei interessantesten: PHP, https://pypi.python.org/pypi/pyshp teext: Python und JavaScript.

Für mein Beispiel werde ich eine kleine Frontend-Anwendung in JavaScript entwickeln, um so meine drei Shapefiles anzuzeigen: die Alpenkonvention, die Steinbockkolonien und die Lawinenkartierung. Dafür verwende ich ein JavaScript-Paket, das ganz einfach Shapefile heisst, sowie Leafletjs, um die Polygone aus den Features anzuzeigen. Aber immer eins nach dem anderen.

Auslesen der Features

Disclaimer: Die folgenden Codebeispiele verwenden ES6 und setzen ein Webpack/Babel-Setup voraus. Man kann sie nicht einfach «copy-pasten», aber sie zeigen auf, wie man mit den zur Verfügung stehenden Tools arbeiten kann.

Zuerst versuche ich das Shapefile der Alpenkonvention zu laden und in die Konsole zu importieren. Das oben erwähnte Shapefile-Paket enthält zwei Beispiele in der README-Datei, weshalb ich einfachheitshalber folgenden Ansatz anwende:

import { open } from 'shapefile'

open('assets/shapefiles/alpine-convention/shapefile.shp').then(source => {
  source.read().then(function apply (result) { // Read feature from the shapefile
    if (result.done) { // If there's no features left, abort
      return
    }

    console.log(result)

    return source.read().then(apply) // Iterate over result
  })
})

Klappt wunderbar! Das Ergebnis sind die beiden Features des Shapefiles in der Konsole. Was ich hier bereits sehen kann, ist, dass die Eigenschaften jedes Features, die normalerweise in separaten Dateien gespeichert sind, hier bereits in den Features enthalten sind. Das kommt daher, dass die Bibliothek beide Dateien abruft und zusammenführt:

Ich werfe jetzt einen Blick auf die Koordinaten, die dem ersten Feature zugeordnet sind. Als erstes fällt mir auf, dass es zwei verschiedene Datenfelder gibt, was auf zwei unterschiedliche Koordinatensätze hindeutet. Ich werde den zweiten aber vorerst ignorieren.

Und hier kommt der erste Fallstrick. Diese Koordinaten sehen überhaupt nicht wie WGS84-Koordinaten (geografische Breite / geografische Länge) aus. WGS84-Koordinaten sollten als Gleitkommazahl abgebildet sein und idealerweise mit 47 und 8 für die Schweiz beginnen. Bei diesem Shapefile sehe ich mich mit einem anderen Koordinatenreferenzsystem (CRS) konfrontiert. Um die Form auf einer Karte korrekt darzustellen oder sie mit anderen Daten zu verwenden, müsste ich diese Koordinaten zunächst in WGS84-Koordinaten umwandeln, aber dazu muss ich erst herausfinden, welches CRS in diesem Shapefile verwendet wird. Da dieses Shapefile vom Bundesamt für Raumentwicklung ARE stammt, ist das verwendete CRS höchstwahrscheinlich CH1903, mit anderen Worten das Schweizer Koordinatensystem.

Um das umzurechnen, braucht es also erst etwas Mathematik. Nach etwas Recherche im Internet habe ich eine JavaScript-Lösung gefunden, mit der ich Umrechnungen zwischen CH1903 und WGS84 machen kann. Ich brauche aber nur ein paar Teile davon, weshalb ich den Code kopiere und etwas umändere.

// Inspired by https://raw.githubusercontent.com/ValentinMinder/Swisstopo-WGS84-LV03/master/scripts/js/wgs84_ch1903.js

/**
 * Converts CH1903(+) to Latitude
 * @param y
 * @param x
 * @return {number}
 * @constructor
 */
const CHtoWGSlat = (y, x) => {
  // Converts military to civil and to unit = 1000km
  // Auxiliary values (% Bern)
  const yAux = (y - 600000) / 1000000
  const xAux = (x - 200000) / 1000000

  // Process lat
  const lat = 16.9023892 +
    (3.238272 * xAux) -
    (0.270978 * Math.pow(yAux, 2)) -
    (0.002528 * Math.pow(xAux, 2)) -
    (0.0447 * Math.pow(yAux, 2) * xAux) -
    (0.0140 * Math.pow(xAux, 3))

  // Unit 10000" to 1" and converts seconds to degrees (dec)
  return lat * 100 / 36
}

/**
 * Converts CH1903(+) to Longitude
 * @param y
 * @param x
 * @return {number}
 * @constructor
 */
const CHtoWGSlng = (y, x) => {
  // Auxiliary values (% Bern)
  const yAux = (y - 600000) / 1000000
  const xAux = (x - 200000) / 1000000

  // Process lng
  const lng = 2.6779094 +
    (4.728982 * yAux) +
    (0.791484 * yAux * xAux) +
    (0.1306 * yAux * Math.pow(xAux, 2)) -
    (0.0436 * Math.pow(yAux, 3))

  // Unit 10000" to 1 " and converts seconds to degrees (dec)
  return lng * 100 / 36
}

/**
 * Convert CH1903(+) to WGS84 (Latitude/Longitude)
 * @param y
 * @param x
 */
export default (y, x) => [
  CHtoWGSlat(y, x),
  CHtoWGSlng(y, x)
]

Das Resultat kann ich nun für meine Hauptanwendung verwenden:

import { open } from 'shapefile'
import ch1903ToWgs from './js/ch1903ToWgs'

open('assets/shapefiles/alpine-convention/shapefile.shp').then(source => {
  source.read().then(function apply (result) { // Read feature from the shapefile
    if (result.done) { // If there's no features left, abort
      return
    }

    // Convert CH1903 to WGS84
    const coords = result.value.geometry.coordinates[0].map(coordPair => {
      return ch1903ToWgs(coordPair[0], coordPair[1])
    })

    console.log(coords)

    return source.read().then(apply) // Iterate over result
  })
})

Daraus ergeben sich die folgenden angepassten Koordinaten:

Das ist viel besser. Mit diesen Koordinaten überlagere ich nun die Karte, damit ich auch tatsächlich etwas sehe.

Dinge sichtbar machen

Dazu füge ich der App Leaflet und die Positron-Lite-Kacheln aus CartoDB hinzu. Das Positron-Lite-Design ist in Hellgrau/Weiss gehalten, was einen guten Kontrast zu den Features bildet, die ich darübergelagert anzeigen werde. OpenStreetMap ist grossartig, aber zu bunt, weshalb die Polygone nicht deutlich genug sichtbar sind.

import leaflet from 'leaflet'
import { open } from 'shapefile'
import ch1903ToWgs from './js/ch1903ToWgs'

/**
 * Add the map
 */
const map = new leaflet.Map(document.querySelector('#map'), {
  zoomControl: false, // Disable default one to re-add custom one
}).setView([46.8182, 8.2275], 9) // Show Switzerland by default

// Move zoom control to bottom right corner
map.addControl(leaflet.control.zoom({position: 'bottomright'}))

// Add the tiles
const tileLayer = new leaflet.TileLayer('//cartodb-basemaps-{s}.global.ssl.fastly.net/light_all/{z}/{x}/{y}.png', {
  minZoom: 9,
  maxZoom: 20,
  attribution: '© CartoDB basemaps'
}).addTo(map)

map.addLayer(tileLayer)

/**
 * Process the shapefile
 */
open('assets/shapefiles/alpine-convention/shapefile.shp').then(source => {
  source.read().then(function apply (result) { // Read feature from the shapefile
    if (result.done) { // If there's no features left, abort
      return
    }

    // Convert CH1903 to WGS84
    const coords = result.value.geometry.coordinates[0].map(coordPair => {
      return ch1903ToWgs(coordPair[0], coordPair[1])
    })

    console.log(coords)

    return source.read().then(apply) // Iterate over result
  })
})

Daraus ergibt sich bereits eine schöne Karte, mit der ich arbeiten kann:

Der nächste Schritt besteht darin, die Karte mit dem Tool zum Laden von Shapefiles zu verbinden. Zu diesem Zweck erstelle ich aus den Koordinaten, die ich vorhin in WGS84-Koordinaten umgewandelt habe, Leaflet-Polygone.

// ...

// Convert CH1903 to WGS84
const coords = result.value.geometry.coordinates[0].map(coordPair => {
  return ch1903ToWgs(coordPair[0], coordPair[1])
})

const leafletPolygon = new leaflet.Polygon(coords, {
  weight: 0.5,
  color: '#757575',
  fillOpacity: 0.3,
  opcacity: 0.3,
})

leafletPolygon.addTo(map)

// ...

Jetzt kann ich mir das Ergebnis ansehen:

Schön!

Das Endergebnis

Durch ein leicht erweitertes User Interface und das Hinzufügen eines Umschalters für die verschiedenen Shapefiles, die ich vorhin heruntergeladen habe, kann ich eine kleine App entwickeln, die mir sowohl die Gefahrenzonen für Lawinen als auch die Lebensräume der Steinböcke in der Schweiz zeigt:

Hier sehen wir die beiden Shapefiles im Einsatz: Die hellgrünen, grünen, gelben und roten Polygone sind Gemeinden, die einer gewissen Lawinengefahr ausgesetzt sind (hellgrün = geringe Gefahr, grün = mittel bis geringe Gefahr, gelb = mittlere Gefahr, rot = hohe Gefahr), die darübergelegten dunkleren Polygone zeigen die Steinbockkolonien.

Diese Karte liefert bereits gute Erkenntnisse über das Verhalten von Steinböcken: Offenbar leben sie nur in Zonen mit geringer bis mittlerer Lawinengefährdung. Sie neigen auch dazu, mittel- und hochgefährdete Gebiete zu meiden. Da Steinböcke alpines Gelände bevorzugen und dort in der Regel eine gewisse Lawinengefahr besteht, ist dies plausibel. Aber jetzt habe ich einige handfeste Daten, die diese Behauptung tatsächlich untermauern!

Die fertige Anwendung präsentiert sich hier in Aktion. Den Quellcode findet man in diesem Archiv.

Die Fallstricke

Obwohl Shapefiles eine gute Möglichkeit darstellen, um GIS-Daten zu bearbeiten, sollte man doch einige Fallstricke vermeiden. Einer, den ich bereits erwähnt habe, ist ein unübliches CRS. In den meisten Fällen erkennt man das bei der Anwendung, aber es empfiehlt sich zu prüfen, welches CRS im Shapefile verwendet wird. Der zweite grosse Fallstrick ist die Grösse der Shapefiles. Wenn Shapefiles direkt im Browser mit JavaScript verwendet werden, kann es vorkommen, dass der Browser abstürzt, weil er die riesige Menge an Polygonen nicht verarbeiten kann. Um das zu vermeiden, gibt es aber verschiedene Lösungen: Man kann das Shapefile entweder vereinfachen, indem man unnötige Polygone entfernt, oder es vorbearbeiten und die Polygone in einer Art Datenbank speichern, während man nur die aktuell sichtbaren Polygone abfragt.

Denkanstösse

Ich persönlich liebe es, mit Shapefiles zu arbeiten. Es gibt viele sinnvolle Anwendungsfälle für Shapefiles und es ist äusserst spannend, sich mit den darin enthaltenen Daten auseinanderzusetzen. Auch wenn es den einen oder anderen Fallstrick zu beachten gibt, ist es im Allgemeinen sehr praktisch mit Shapefiles zu arbeiten, da man mit wenig Aufwand etwas erfolgreich entwickeln und zum Laufen bringen kann. Die OpenData-Community verwendet sehr häufig Shapefiles, weshalb sie eine etablierte Norm sind. Ausserdem gibt es zahlreiche Bibliotheken und Tools, die die Arbeit mit Shapefiles so fantastisch machen, wie sie ist.