Tutorial: Creating buffered country POLYs for OpenStreetMap data processing


OpenStreetMap represents a lot of data. If you want to import the entire planet into a PostGIS database using osmosis, you need at least 300GB of hard disk space and, depending on how much you spent on fast processors and (more importantly) memory, a lot of patience. Chances are that you are interested in only a tiny part of the world, either to generate a map or do some data analysis. There’s several ways to get bite-sized chunks of the planet – take a look at the various planet mirrors or the cool new Extract-o-tron tool – but sometimes you may want something custom. For the data temperature analysis I did for State of the Map, I wanted city-sized extracts using a small buffer around the city border. If you want to do something similar – or are just interested in how to do basic geoprocessing on a vector file – this tutorial may be of interest to you. Instead of city borders, which I created myself from the excellent Zillow neighborhood boundary dataset, I will show you how to create a suitably generalized OSM POLY file (the de facto standard for describing polygon extracts used by various OSM tools) that is appropriate for extracting a country from the OSM planet with a nice buffer around it.

Let’s get to work.

Preparing Quantum GIS

We will need to add a plugin that allows us to export any polygon from your QGIS desktop as an OSM POLY file. We can get that OSM POLY export plugin for Quantum GIS here.

Unzip the downloaded file and copy the resulting folder into the Python plugins folder. On Windows, if you used the OSGeo installer, that might be

C:\OSGeo4W\apps\qgis\python\plugins

See here for hints where it may be for you.

The plugin should now appear in the Quantum GIS plugin manager (Plugins > Manage plugins…).
If it is not selected, do that now and exit the plugin manager.

Getting Country Borders

Easy. Download world borders from http://thematicmapping.org/downloads/world_borders.php

Unzip the downloaded file and open it in QGIS:

Geoprocessing step 1: Query

Open the Layer Query dialog by either right-clicking on the layer name or selecting Query… from the Layer menu with the TM_WORLD_BORDERS-0.3 layer selected (active).

Type “ISO2″ = “US” in the SQL where clause field and run the query by clicking OK.

Geoprocessing step 2: Buffering

The next step is to create a new polygon representing a buffer around an existing polygon. Because we already queried for the polygon(s) we want to buffer, there’s no need to select anything in the map view. Just make sure the TM_WORLD_BORDERS-0.3 layer is active and select Vector > Geoprocessing Tools > Buffer(s):

Make sure the input vector layer is TM_WORLD_BORDERS-0.3. Only the query will be affected, so we’re operating on a single country and not the entire world.

For Buffer distance, type 1. This is in map units. Because our source borders file is in EPSG:4326, this corresponds to 1 degree which is 69 miles (for the longitudinal axis, that measurement is only valid at the equator and decreases towards the poles). This is a nice size buffer for a country, you may want something larger or smaller depending on the size of the country and what you want to accomplish, so play around with the figure and compare results. Of course, if your map projection is not EPSG:4326, your map units may not be degrees and you should probably be entering much bigger values.

Select a path and filename for the output shapefile. Do not select ‘Dissolve buffer results’. The rest can be left at the default values. Push OK to run the buffer calculation. This can take a little while and the progress bar won’t move. Then you see:

Click Yes. Now we have a buffer polygon based on the US national border:

Geoprocessing step 3: Generalizing

We’re almost done, but the buffer we generated contains a lot of points, which will make the process of cutting a planet file slow. So we’re going to simplify the polygon some. This is also a QGIS built-in function.

Select Vector > Geometry tools > Simplify geometries:

Make sure your buffer layer is selected as the input. Set 0.1 (again, this is in map units) as the Simplify tolerance. This defines by how much the input features will be simplified, the higher this number, the more simplification.

Select a destination for the simplified buffer to be saved. Also select Add result to canvas. Click OK:

This dialog may not seem very promising, but it has worked. Also, I have sometimes gotten an error message after this process completes. Ignore these if you get them.

Geoprocessing step 4: resolving multipolygons

Now, if your simplified country border consists of multiple polygons (as is the case with the US) we have a slight problem. The POLY export plugin does not support multipolygons, so we need to break the multipolygon into single polygons. And even then, we will need to do some manual work if we want OSM .poly files for all the polygons. This is because the plugin relies on unique string attribute values to create different  POLY files, and we do not have those because the polygons we are using are all split from the same multipolygon. So we need to either create a new attribute field and manually enter unique string values in it, or select and export the parts to POLY files one by one and rename the files before they get overwritten.

Finale: Export as POLY

I am going to be lazy here and assume I will only need the contiguous US, so I select the corresponding polygon. After that I invoke the plugin by selecting Plugins > Export OSM Poly > Export to OSM Poly(s):

The plugin will show a list of all the fields that have string values. Select ISO2 and click Yes. Next you will need to select a destination folder for your exported POLY files. Pick or create one and push OK.

This is it! Your POLY files are finished and ready to be used in Osmosis, osmchange and other tools that use it for data processing.

By the way: you can’t load POLY files into JOSM directly, but there’s a perl script to convert POLY files to OSM files that I used in order to visualize the result.

Dataverkenning AppsForNoordHolland: INSPIRE


The Dutch province of Noord Holland and Dutch open data advocates HackDeOverheid are holding an open data app competition, with a coding event happening next week at a very special venue. I was asked to shed some light on the geospatial data sources made available, which is what this post in Dutch is about. More info at the website.

Komende zaterdag vindt aan het andere einde van de wereld – nou ja, voor mij dan – de eerste publieke bijeenkomst code event in het kader van Apps For Noord-Holland plaats. De provincie Noord-Holland en een handvol andere overheden en semi-overheden hebben databronnen opengesteld waarmee we aan de slag kunnen.

Veel van die data is geografisch – de gegevens hebben betrekking op een specifieke plek; een punt, een lijn, een vlak. De wereld van de geodata is een beetje een aparte; als je met geodata wat wil aanvangen moet je een klein beetje weten van de specifieke bestandsformaten en zaken als coordinaatsystemen en projecties. Daar zijn gelukkig uitstekende bronnen voor, waarvan ik er in de loop van deze dataverkenning een paar zal geven.

Kijkend naar de beschikbare databronnen ligt bij de meeste voor de hand wat je kunt verwachten: zwemwaterkwaliteit, bedrijventerreinen, wegwerkzaamheden. Ik neem er in deze verkenning een beetje een obscure uit: INSPIRE. Eerst maar eens proberen om in een paar zinnen uit te leggen wat INSPIRE is. Daarna gaan we naar de data kijken. En passant geef ik wat tips voor (open source) software om het bekijken en verwerken van geodata behapbaar te maken – ook voor niet-GISsers.

INSPIRE

In 2007 heeft de Europese Commissie de INSPIRE-richtlijn vastgesteld. Hierin wordt bepaald dat alle lidstaten landsdekkende, actuele data beschikbaar moeten stellen over 34 thema’s, die een breed spectrum beslaan van geografische namen tot geologie. Het doel van INSPIRE is om de beleidsvorming over mileu-gerelateerde thema’s te ondersteunen, maar de gegevens zijn in principe voor iedereen beschikbaar. Niet allemaal tegelijk trouwens; de thema’s die onder Annex I vallen moeten eind volgend jaar beschikbaar zijn, maar de thema’s die onder Annex III vallen uiterlijk pas in 2019. Daar kopen we dus nog even niks voor, maar de lidstaten beginnen al wel met Annex I-thema’s over de brug te komen. Bijvoorbeeld dus Noord-Holland, waar we nu naar aan het kijken zijn.

Kaartservice schmaartservice

Je kunt de data volgens opgave hier bekijken in een interactieve kaartservice. Een beetje geduld is wel gevraagd, ik moest ruim 10 seconden wachten op mijn kaart. Verplaatsen en inzoomen gaat daarna wel wat sneller.

Zo’n kaartservice is geinig om de data even snel te bekijken, maar er zit geen legenda bij en je kunt de data vanaf hier ook niet downloaden of verder gebruiken. Laten we deze webpagina dus maar laten voor wat-ie is en de data downloaden op de datapagina. Ongeduldig? Hier (32MB .zip) is een directe link.

Als we het ZIP-bestand uitpakken krijgen we een brij aan bestanden. Tijd voor een spoedcursus SHAPE – dat is het formaat waarin we deze data geleverd krijgen.

SHAPE

SHAPE-bestanden – formeel ESRI Shape naar de softwarebouwer die het formaat heeft ontworpen als uitwisselingsformaat voor geodata – kunnen door alle GIS-software worden ingelezen. Het is daarmee een universeel uitwisselingsformaat. Met GIS-software bedoel ik dan de software die professionals gebruiken. Google Maps en Google Earth kunnen er niet rechtstreeks wat mee. Dat is meteen ook een van de belangrijkste beperkingen. Tegenwoordig komen formaten als KML en GeoJSON meer in zwang onder de neogeografen. Omzetten is gelukkig niet zo’n probleem, dat zien we straks.

De naam ‘Shape-bestand’ is eigenlijk misleidend. Een ‘Shape-bestand’ bestaat uit minimaal drie bestanden die onlosmakelijk met elkaar zijn verbonden, en nog een stuk of wat optionele bestanden:

  • Data.SHP – hierin zit de eigenlijke geodata: lijnen, punten, vlakken (‘e’en soort per Shape-bestand).
  • Data.DBF – hierin vind je de ‘attributen’, de metadata voor elk object. Voor een bestand met landgebruik zou je hierin bijvoorbeeld een veld ‘type’ kunnen aantreffen met waarden als ‘industrie’, ‘winkels’, ‘wonen’. De attributen worden opgeslagen in het oeroude Dbase-formaat met alle beperkingen (veldnamen maximaal 8 karakters bijvoorbeeld) van dien.
  • Data.SHX – dit bestand is een index op de objecten in het SHP-bestand, vergelijkbaar met een index op een gewone database. Dit bestand zorgt ervoor dat het bevragen van de data efficienter gaat.

Verder is er nog een hele ris optionele bestanden behorend bij ‘een Shape-bestand’ die vooral door specialistische software worden gebruikt. Ik laat die even voor wat ze zijn, behalve eentje die belangrijk is om iets van te begrijpen:

  • Data.PRJ – dit bestandje bevat een zogenaamde projectiestring, waarin wordt gedefinieerd hoe de driedimensionale aarde op een tweedimensionaal vlak wordt geprojecteerd. Er bestaan honderden verschillende projecties en bijbehorende coordinaatstelsels. In Nederland wordt het RD-stelsel gebruikt. Wereldwijd kom je meestal ‘Plate Carree‘ of ‘Spherical Mercator‘ tegen. Het omrekenen tussen projecties is meestal een van de eerste struikelblokken voor mensen die voor het eerst met geografische data werken.

En dan nu de data…

Goed, hehe, zullen we dan nu eindelijk eens naar de data gaan kijken? De eerste vraag is hoe dan. Met Google Earth kun je alvast geen Shape-bestanden openen. Gelukkig biedt de open source-wereld een hele goede oplossing: Quantum GIS, of kortweg Qgis. Dit gratis pakket kan 80% van de dingen die grote dure GIS-pakketten van vele duizenden euroos kunnen. Bij die 80% hoort zeker de functionaliteit die wij nodig hebben. Ik ga een paar dingen laten zien:

  • Qgis downloaden en installeren
  • Data Inladen
  • Classificeren
  • Projectie instellen
  • Exporteren in een bruikbaar formaat

Daarmee zou je genoeg bagage moeten hebben om deze – en ook andere – geodatabronnen te beoordelen en naar je hand te zetten.

Qgis downloaden en installeren

Quantum GIS is er in smaken voor Windows, Mac en Linux. Je kunt de nieuwste versie oppikken van http://www.qgis.org/. Er is geen automatische update-functionaliteit, dus als je Qgis vaker gaat gebruiken moet je af en toe eens terugkomen voor de nieuwste versie. Mijn ervaring is dat dat wel de moeite waard is.

Installeer Qgis volgens de instructies en starten maar.

Als eerste zie je een leeg scherm. Geen wizard of niks. Je bent meteen in het diepe van een GIS-pakket. Links is ruimte voor de kaartlagen. In een GIS wordt een kaartbeeld altijd opgebouwd uit lagen. Achter elke laag (layer) gaat een dataset schuil. Dat kan een bestand zijn of een database-tabel. Lagen kunnen elkaar overlappen – dat is zelfs heel gebruikelijk. We zullen zo zien hoe dat uitpakt. Laten we maar wat lagen gaan toevoegen.

Data inladen

Een GIS kent twee soorten lagen: raster en vector. Rasterlagen bestaan uit plaatjes; vectorlagen bestaan uit geometrische definities van lijnen, punten en vlakken. Photoshop versus Illustrator. Shape-bestanden zijn altijd vectorlagen. Kies dus uit het menu Layer → Add Vector Layer…

Je kunt meerdere lagen ineens toevoegen, zodat je in een klap alle zestien Shapes uit de INSPIRE-dataset van Noord-Holland op je scherm hebt – wel snel maar visueel wordt het nogal een zootje. Voor nu beginnen we er dus maar eens met eentje: DPDATA_geluid_prov_wegen.shp. Door het glanzende gebrek aan databeschrijvingen moet ik wat aannames doen op basis van de bestandsnaam. Ik denk dat het gaat om geluidscontouren van de provinciale wegen.

De kenner weet dat het hier om Noord-Holland gaat, maar verder zijn we nog niet zo gek veel wijzer. We kunnen een paar eenvoudige trucjes uithalen om iets meer over de data te weten te komen. Ten eerste kan je dubbelklikken op de naam van de laag in de linkerhelft van het scherm. Je krijgt dan een informatievenster met verschillende tabs, waaronder de attribuutvelden en de metadata. Ook kun je afzonderlijke objecten inspecteren; zoek naar het pijltje met ‘i’ in de knoppenbalk bovenaan en selecteer een object in de kaart. Je krijgt een venster met de attributen voor dat veld. Een addertje onder het gras is hier dat je alleen objecten kunt selecteren van de actieve laag – dat is de laag die links geselecteerd is.

Dit is allemaal leuk en aardig, maar bij een grote hoeveelheid data zoals we hier hebben krijg je er met deze tools nog geen gevoel bij. Het al genoemde gebrek aan beschrijvende gegevens bij de data noopt tot wat meer grasduinen in de data. Dat brengt ons bij de volgende stap: wat orde in de chaos aanbrengen.

Classificeren

Laten we eens gek doen en alle lagen ineens inladen. Dat duurt even – althans op mijn laptopje. Eenmaal ingeladen kun je zien dat Qgis alle lagen een willekeurige kleur heeft gegeven. Aangezien lagen elkaar overlappen zie je niet eens alle data, maar het is al een beste soep.

Met de vinkjes naast de lagen kun je individuele lagen zichtbaar en onzichtbaar maken. Je kunt lagen in de lijst ook verslepen naar een hogere of lagere positie; ze worden dan ook in de kaart verder bovenop of onderop getekend. Maar daarmee wordt het alleen maar een anders uitziende chaos. Het toverwoord om daar een einde aan te maken is classificeren!

Classificeren is het afhankelijk maken van de stijl (vulkleur, lijn, arcering) van de waarde van een van de attributen. Stel, je hebt een attribuut ‘hoogte’, dan kun je die laag classificeren zodat lager gelegen gebieden groen ingekleurd worden en hoger gelegen gebieden bruin. Je kunt op twee manieren classificeren: op individuele waarden – nuttig bij beschrijvende attributen, zoals type (denk aan wonen, industrie, winkels), of op intervallen voor numerieke attribuutwaarden zoals hoogte of bevolkingsdichtheid.

De magie schuilt in het venster met lageneigenschappen dat we al eerder tevoorschijn toverden. Hierin vind je de tab ‘Style’, die standaard op ‘Single Symbol’ staat. Als we eens kijken naar de laag met geluidscontouren die we al eerder bij de horens hadden, dan leent die zich uitstekend voor een classificatie in intervallen. Roep de eigenschappen voor die laag op, selecteer de ‘Style’ tab en selecteer ‘Graduated’ in plaats van ‘Single Symbol’. Verander ook de instelling ‘mode’ van ‘Equal Interval’ in ‘Natural Breaks’. Klik nu op ‘Classify’ en de voorgestelde classificatie op basis van de gemaakte instellingen komt tevoorschijn. Als je deze bevestigt wordt ook het kaartbeeld bijgewerkt en zie je een fraai kleurverloop. Zoom eventueel een beetje in om het beter te zien.

Misschien denk je ‘allemaal Spielerei’ – maar dat is het zeker niet! Een goede visualisatie van de data is essentieel om een idee te krijgen van het verhaal dat de data vertelt. Dat verhaal is de basis voor jouw gouden app-idee en heeft bijna altijd een verrassende plotwending. Die ontvouwt zich pas als je wat gaat spelen met de data, combineren, visualiseren. Dat is wat we nu doen. Dit is het resultaat van mijn grasduinsessie van eerder vandaag:

Hier kun je dit resultaat als Qgis-project downloaden. Als je dat bestand kopieert naar je map met INSPIRE-Shape-bestanden en opent in Qgis krijg je als het goed is precies deze visualisatie op je scherm. Als het niet werkt moet je waarschijnlijk paden aanpassen in het projectbestand. Zoek naar dergelijke regels:

<datasource>/home/mvexel/Downloads/ShapesInspire/DPDATA_CHW_ALG_VLAKKEN.shp</datasource>

Vervang waar nodig. Maar ik denk dat je die tijd beter kunt besteden aan zelf prutten met de data eerlijk gezegd.

Projectie instellen

Goed, terug naar het saaie werk. Als je goed hebt opgelet heb je misschien gezien dat de Shape-bestanden van de provincie Noord-Holland geen .PRJ-bestanden hebben. Foei Provincie! Hoe kunnen we zonder berschrijvende gegevens of .PRJ-bestand nu weten welke projectie de bestanden hebben? Dit is essentiele informatie bij een geodatabestand.

De context en mijn ervaring vertellen me dat de provincie het Nederlandse RD-stelsel heeft gebruikt, maar dat had er best ergens bij mogen staan. We moeten dit nu eerst aan Qgis vertellen, voordat we opdracht kunnen geven om de data naar een andere projectie om te rekenen en te exporteren.

Ook dit gaat in het venster van de laag-eigenschappen – ja, we moeten dit per laag instellen. Dit gaat met de instelling ‘Specify CRS’ in de ‘General’ tab. Stel daar in dat de laag in ‘EPSG:28992 – Amersfoort / RD New’ (sic) is gedefinieerd.

Exporteren in een bruikbaar formaat

Hier gaat het uiteindelijk allemaal om! We willen de data in zo’n formaat dat we er in onze webapplicatie, Android- of iPhone-app wat mee kunnen aanvangen. Mijn geinformeerde gok zegt dat het Shape-formaat niet voldoet, maar bijvoorbeeld KML wel. Laten we dus proberen de data naar dit formaat te exporteren.

Er is nog meer aan de hand: naast het juiste bestandsformaat is ook de juiste projectie essentieel! Als je jezelf geen gedonder op de hals wilt halen later bij het coden van je app – en gedonder met projecties, dat wil je niet – dan zorg je nu voor de juiste projectie van je data. De brave borsten van de provincie hebben keurig de Nederlandse RD-standaard aangehoden voor deze dataset – overigens niet in overeenstemming met de INSPIRE-richtlijn, maar dat vlug terzijde. Wij willen voor onze app graag de Google-standaard Spherical Mercator (of de oude standaard Plate Carree, de procedure is hetzelfde).

Selecteer eerst de laag die je verder wilt gaan gebruiken en vervolgens uit het menu Layer → Save as…

Hier kiezen we een doelbestand, KML als uitvoerformaat, en we specificeren de projectie als ‘Google Mercator’. Gebruik de zoekfunctie – het is een lange lijst met projecties! Geen wonder dat het, ook voor ervaren GISsers, een bron van verwarring blijft.

Druk nu op OK en je KML wordt gebakken!

Open het resultaat in een ander programma – bijvoorbeeld Google Earth – om te testen of het allemaal een beetje klopt:

Let niet op mij – ik moet nog iets instellen met fonts geloof ik. Maar het gaat erom dat de data klopt! Merk op dat de fraaie kleurtjes etc. niet mee zijn gekomen. Dat kan misschien wel maar daar gaat het niet om. Jouw app zal de data misschien helemaal niet willen laten zien, maar bijvoorbeeld gebruiken als onzichtbare drempel terwijl je je met je telefoon door het land verplaatst, of om een waarde te berekenen op basis van je huidige lokatie. En als je wil visualiseren, zal jouw kaartframework daar weer eigen voorzieningen voor hebben.

Mooi werk! We hebben nu KML-data in een projectie waar we allemaal wat aan hebben.

Conclusie

Pfoe. Het is uitgedraaid op een hele GIS-tutorial. Dat was helemaal niet de bedoeling maar wel leuk om te doen! En volgens mij ook nuttig. Met deze skills kun je een heleboel veel voorkomende uitdagingen met geodata te lijf. Veel succes komende zaterdag en daarna met de App-competitie!

Meer informatie

  • Een handleiding voor Qgis vind je hier.
  • Als je eens een projectie moet opzoeken, dan doe je dat op SpatialReference.org
  • Andere open source GISsen die ik ken zijn GRASS, gvSIG en uDig, maar ik gebruik zelf het liefste Qgis vanwege stabiliteit en toegankelijkheid. Een nog grotere lijst vind je hier.
  • Wil je geodata snel en makkelijk op het web beschikbaar maken? Kijk dan zeker eens naar de OpenGeo Suite. Deze link is naar de gratis Community-editie maar er is ook een Enterprise-editie met support enzo.
  • Over projecties kun je eindeloos lezen. De documentatie van Google Earth bevat een toegankelijk verhaal over projecties en datums (daar hadden we het nog niet eens over gehad…)