hic sunt dracones

screen-shot-2017-02-26-at-8-05-05-pm

It was quite common to decorate uncharted areas with dragons or sea monsters on medieval maps. Some came with a latin warning, HIC SUNT DRACONES. But making customised, interactive maps in R is nothing to be feared, here’s a quick guide.

First things first, we’ll need some data. abs.stat is a good place to start. Pick something of interest and export it. Keep an eye on what geographic level the data is at – Statistical Levels (SA1-SA4) are pretty common, and a few are available at the post code (POA), state suburb (SSC) or Local Government Area (LGA). For this post I’m using data from the 2011 Census by LGA.

We’ll also need to get ourselves some matching polygons. The ABS offers both shapefiles and mapinfo files, but we’ll need something that will play nicely with leaflet. Download the shapefile zip then head to mapshaper.org and drag and drop the file and export in the GeoJson. No doubt there’s a way to do this in R too, but mapshaper.org is just too darn easy.

Load it into R, along with a few packages. You’ll also want to simplify the polygons – this will make everything load faster and the final file will be much smaller.

# packages
install.packages("rgdal")
install.packages("leaflet")
install.packages("rmapshaper")

# load data
## the polygons (spatial data)
abs_lga <- readOGR("LGA_2011_AUST.json")
## data of interest (non-spatial data)
census_lga <- read.csv("ABS_CENSUS2011_LGA.csv")

# simplify
simp_abs_lga <- ms_simplify(abs_lga)

It might take a little while for the geojson to load. Once it’s done you can take a look.

head(simp_abs_lga@data)
plot(simp_abs_lga)

It should look something like this:

LGA_CODE11 LGA_NAME11 STE_CODE11 STE_NAME11 AREA_SQKM rmapshaperid
0 10050 Albury (C) 1 New South Wales 305.930832 0
1 10110 Armidale Dumaresq (A) 1 New South Wales 4230.834878 1
2 10150 Ashfield (A) 1 New South Wales 8.281188 2
3 10200 Auburn (C) 1 New South Wales 32.476919 3
4 10250 Ballina (A) 1 New South Wales 484.716361 4
5 10300 Balranald (A) 1 New South Wales 21693.124541 5

And this should pop up in the plots pane:

screen-shot-2017-02-26-at-6-11-00-pm

Because this is LGAs, there’s no polygons in the ACT or the outback.

Now we’ll merge the polygons with the other data. I’ve subsetted it as there were lots of different metrics in the file.

# subset non-spatial data
age_lga <- census_lga[census_lga$Selected.Medians.and.Averages == "Median age of persons", c(7,8,13)]
# merge the spatial and non-spatial dataframes
simp_abs_lga_age <- merge(simp_abs_lga, age_lga, by.x = 'LGA_CODE11', by.y = 'REGION')

Picking a palette is fun and easy because the leaflet package uses colorbrewer. But, if you like you can define your own palette here using hex codes. I’ve also made a custom marker with a popup for our uncharted waters (well not really but..)

# palette
bins <- c(0, 20, 25, 30, 35, 40, 45, 50, Inf)
pal <- colorBin("Reds", domain = simp_abs_lga_age$Value, bins = bins)

# custom icon
Icon <- makeIcon(iconUrl = "PKMN147b.png",
  iconAnchorX = 22, iconAnchorY = 22)

# marker locations
long <- c(144.2631, 145.0631, 143.3631)
lat <- c(-38.8136, -38.9136, -39.1136)
dragons <- as.data.frame(cbind(long,lat))

# popups (remove the spaces in the html tags)
popup <- paste(sep = "<br/ >",
               "<b >Here be dragons</b >")

Now on to the mapping. The leaflet package uses pipes to split out each layer for the map. The example below is pretty basic, there’s heaps more things you can do with this package.

You’ll also need to pick out a tile provider, there’s a number of free ones supported by leaflet. I really need to find some good data for SpinalMap, maybe the metal-ness of each postcode.

map <- leaflet(simp_abs_lga_age)%>%
  setView(144.9631, -38.5136, 7)%>%
  addProviderTiles(providers$Esri.WorldStreetMap) %>%
  addPolygons(
    fillColor = ~pal(Value),
    weight = 1,
    opacity = 1,
    color = "white",
    dashArray = "3",
    fillOpacity = 0.7
   )%>%
   addMarkers(lng = dragons$long, lat = dragons$lat,
             icon = Icon,
             popup = popup
   )%>%
   addLegend("bottomright", pal = pal, values = ~Value,
          title = "Median Age by LGA",
          opacity = 1
)
map

And that’s it. That’s all there is to it. It’ll pop up in the plot window and then up to you what you do with it.

Here it is, in all it’s interactive glory.

Rpubs seems to be the easiest way to get it out there, but you can export as html or use the htmlwidgets package too.

Now go slay some (map) dragons.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s