Reading and plotting OSM data in QGIS

This post will be a little different than the previous posts. R gets a rest for the time being, with this post focusing on some GIS tools.

GIS, or Geographic Information System seems to fallen out of flavour amongst those with 'data science' leanings. R is a fantastic toolset, with some great geospatial packages that are built on solid GIS libraries, but sometimes I feel it just doesn't compare to using programs made to fit for geo work. This is why this post will focus on visualising data using QGIS (http://www.qgis.org).

Don't be alarmed by this. QGIS is Open Source and well supported. It has a python console that allows the user to script things as required, and also hooks up to other GIS software such as GRASS, GDAL and SAGA. Goto the link above to check more out about it.

The purpose of this post will be to show how easy it is to download OpenStreetmap, or OSM, data into QGIS and have a look around. For those unaware of what OSM is, check it out here.

Firstly, I'm on Ubuntu. Which means scripting things is a breeze, we can download the OSM data easily from geofabrik in a one liner. We'll get data for The Netherlands:

 wget -c http://download.geofabrik.de/europe/netherlands-latest.osm  

This .osm file can be quite large, and compressed versions are available. But we consider it as almost a source data file that we will import into a database. I've used PostgreSQL along with the GIS extension PostGIS for the last two years now. It works great and can do many great things.

In order to import the .osm file into PostgreSQL + PostGIS we require a library, GDAL. This is a well known library that operates across multiple platforms - you may have used this package in R. Again, this is open source. We're going to use the ogr2ogr program from GDAL. The manual describes the function of this program quite well:

This program can be used to convert simple features data between file formats performing various operations during the process such as spatial or attribute selections, reducing the set of attributes, setting the output coordinate system or even reprojecting the features during translation.

Assuming you've installed things correctly, and you have access working PostgreSQL + PostGIS database, you can simply insert the .osm data directly in.

 ogr2ogr -f "PostgreSQL" PG:dbname=data netherlands-latest.osm -gt 65536 -overwrite --config PG_USE_COPY YES --config OVERWRITE YES --config SPATIAL_INDEX NO -progress  

The command instructs ogr2ogr to insert data from 'netherlands-latest.osm' into 'PostgreSQL' database 'data'. I'll leave the additional parameters for you to investigate, but quickly mention '--config PG_USE_COPY YES'    as it greatly enhances the speed of the operation on certain systems.

Now we can fire-up QGIS and investigate what we've done. You'll need to click Layer>Add Layer>Add PostGIS Layers. Upon connecting to the database in QGIS, we should see five tables with vector geometries. These are:

1. points
2. lines
3. multilinestrings
4. multipolygons
5. other_relations

We're going to add roads to our map, so we select lines and filter with the following code:

 "highway" in ('motorway','motorway_junction','motorway_link','primary','primary_link','trunk','trunk_link')  

This allows us to add only roads in the lines table that are of type motorway, primary or trunk. This filter is important on low-spec'd machines are often it's possible to import more data than your machine can handle. Let's see the outcome.

Image 1. Adding a PostGIS geometry table with a filter in QGIS

Image 2. Motorways, Primary and Trunk roads of The Netherlands via OSM

Lovely, that's a lot of data. You can now add in other tables, shapefiles, text files - whatever you can find in ogr2ogr. I prefer uploading everything to PostGIS so that I can easily reproduce my work.

With a little work, you can achieve some beautiful maps and unlock some incredibly useful data from OSM. Below is just a small sample of the potential of OSM data with QGIS, some additional SQL and creative styling.

Image 3. Bank branch locations in Amsterdam and surrounding region

This post barely touches on any of the fantastic things you can do with QGIS, PostGIS or OSM data, and doesn't even mention the interaction with R. Hopefully I'll have time to show more cool things this software combination has to offer.

Cheers

Plot contour polygons in Leaflet using R

cb

08/10/2014

This post illustrates how easy it is to visualise geospatial data using R. In particular, we will calculate a 2d density estimate of our geo data using the KernSmooth package, transform the data using SP, then finally visualise in Leaflet using the LeafletR and RColorBrewer packages.
The data I’m using is the Melbourne UrbanForest dataset and has been released by the Victorian Government. You can check it out here. Let’s load it into R.

# Load packages
library(leafletR)
library(sp)
library(KernSmooth)
library(RColorBrewer)

# Read in the data and remove trees of unknown age
data = read.csv("~/Downloads/Melbourne_s_Urban_Forest_Tree_data.csv")
# Remove observations with no Age Description
data = subset(data,data$Age.Description != '')

Next, we view the estimated 2D Density. I’m effectively binning over a gridcell of 0.00225 longitude by 0.00225 latitude. This roughly corresponds to a 250m x 250m grid. Normally I’d reproject data into a projection that is in meters (GDA94 / Australian Albers works well) but I’m more interested in just visualising at the moment.

# Apply 
d2d = bkde2D(cbind(data$Longitude,data$Latitude),bandwidth=c(0.00225,0.00225))
# Visualise
contour(d2d$x1,d2d$x2,d2d$fhat)

Cool! That looks kinda nice. We’re on the right track. The result of d2d is gridpoints however, we still need to bin these points into buckets.
The code below creates a set of contour lines from the grid data, then transforms the data into an SP object. Quite a bit is going on here; it’s important to think about the underlying datasets we’re working with, we move from grid points in 3d, to linestrings in 2d, ‘fill’ the lines to form independent polygons and then finally collect these polygons in one object. It is important to follow this procedure - polygons within polygons can form islands and islands within islands (sounds like Inception the movie right?) if collected incorrectly. So yes take care.

# Create linestrings
lines = contourLines(x=d2d$x1,y=d2d$x2,z = d2d$fhat,nlevels = 8)

# Create independent polygons within a list
dd1 = sapply(1:length(lines),function(i) Polygon(as.matrix(cbind(lines[[i]]$x,lines[[i]]$y))))

# Merge all independent polygons into a Polygons object (this contains multiple polygons)
dd2 = sapply(1:length(lines),function(i) Polygons(list(dd1[[i]]),i))

# Don't forget to remember the contour value for each polygon - we store it into a dataframe for use in the next step
poly_data = data.frame(Value = sapply(1:length(lines),function(i) lines[[i]]$level))

# Merge the Polygons object dd2 with the dataframe containing the contour level data, poly_data.
dd3 = SpatialPolygonsDataFrame(SpatialPolygons(dd2),data = poly_data)

Phew, that’s over. This step is a little confusing at first, there is some wrangling around to Let’s store our contour values (we will use them in the styling) and convert the polygons to JSON:

# Convert our dd3 SpatialPolygonDataFrame object to JSON
dd_json = toGeoJSON(dd3,name="MelbourneTree")

# Store the unique levels of the contours, this will come in handy for colouring
values = unique(sapply(1:length(lines),function(i) lines[[i]]$level))

Almost there... The LeafletR package in R is very powerful. Just one line creates a beautiful Leaflet .html file that can be uploaded to the web or shared around work. Styling, as you can see is also straightforward when you use RColorBrewer to generate a nice palette.

# Create a style for the Leaflet map
sty = styleCat(prop="Value",val=values,style.val=brewer.pal(length(values),"Greens"),leg = "Tree Cover")

# Create the map object. This will automatically create a .html file on your machine
map = leaflet(dd_json,base.map = "osm",popup="Value",style=sty)

It’s a little effort to get the data into the right place, but the results are definitely worth it =).
Cheers


EDIT: Removed some typos!