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!