From static maps to dynamic maps

So, last week, we saw how we could plot shapefiles in ggplot2. This is great! But, oddly, sometimes static. What if you want to no only put shapefiles on a map, but add other information (say, sites), and then make the whole thing interactive - like a google map! You can!

Adding temperature data to SpatialPolygons

So, a shapefile has some inherent properties in and of itself so that it can be plotted.

library(sp)
library(rgdal)
## rgdal: version: 1.1-3, (SVN revision 594)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.1, released 2014/09/24
##  Path to GDAL shared files: /usr/local/share/gdal
##  Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.2-0
ecoregions <- readOGR("./join_maps/MEOW-TNC/meow_ecos.shp", layer="meow_ecos")
## OGR data source with driver: ESRI Shapefile 
## Source: "./join_maps/MEOW-TNC/meow_ecos.shp", layer: "meow_ecos"
## with 232 features
## It has 9 fields
plot(ecoregions)

And as we pointed out last week, it’s a SpatialPolygonsDataFrame, so, we can add some additional information to it, like temperature change. Let’s get the slope in temperature change by region.

library(lubridate)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
#get the linear slope of temperature change
#in degrees per decade
sst_data <- read.csv("./join_maps/hadsst_regions.csv",
                     stringsAsFactors=FALSE) %>%
  mutate(Date = as.POSIXct(parse_date_time(DateName, orders="ymd"))) %>%
  group_by(ECOREGION) %>%
  do(Temp_Slope = lm(tempC ~ scale(Date, scale=FALSE), data=.)) %>%
  mutate(Temp_Slope = coef(Temp_Slope)[2]) %>%
  #deal with scaling - second/min, min/hour, hr/day, day/year, year/decade
  mutate(Temp_Slope = Temp_Slope*60*60*24*356*10) %>%
  ungroup()

Then join it to the first SpatialPolygonsDataFrame in a new object

ecoregions_temp <- ecoregions
ecoregions_temp@data <- full_join(ecoregions@data, sst_data)
## Joining by: "ECOREGION"
## Warning in outer_join_impl(x, y, by$x, by$y): joining character vector and
## factor, coercing into character vector
#subset it down
ecoregions_temp <- subset(ecoregions_temp, !is.na(ecoregions_temp$Temp_Slope))

Interactive Maps

So, how can we visualize this in a dynamic way? Leaflet!

library(leaflet)

leaflet() %>%
  addTiles() %>%
  addPolygons(data=ecoregions_temp)

Fantastic! We now have regions on a map. The first function, leaflet, initialized the map. The second, addTiles pulls up a tileset for the map from the web. The third adds polygons from a SpatialPolygons object. And note that it works using pipes.

How can we add in the temperature information? For that, the syntax is a bit funky with the colorNumeric argument and the need to reverse our data to accomodat the colorbpalatte.

#Reverse for palette
ecoregions_temp$Temp_Slope_Rev <- -1*ecoregions_temp$Temp_Slope

#note, order of slopes is fine for RdBu
tempMap <- leaflet() %>%
  addTiles() %>%
  addPolygons(data=ecoregions_temp, stroke=FALSE,
  fillOpacity=0.5,
  color = ~colorNumeric("RdBu", ecoregions_temp$Temp_Slope_Rev)(Temp_Slope_Rev))

tempMap

Neat, no?

What else can I add to a map?

Well, you can add a lot! Let’s say, for example, I had a list of sites in these ecoregions, and maybe even some information about them - say the long-term trend in kelp populations.

kelp_slopes <- read.csv("./join_maps/kelp_slopes.csv", stringsAsFactors=FALSE)

Now, let’s add the sites onto this map so we can see if sites in warming or cooling regions are increasing or decreasing in abundance.

Let’s use BrBG as our palette so brown = decrease, green = increase. And we’ll add a popup so that you can click on points, just to help the viewer.

tempMap %>%
   addCircleMarkers(data = kelp_slopes, 
             lng = ~Longitude, lat =  ~Latitude,
             fillOpacity=1,
             color = ~colorQuantile("BrBG", kelp_slopes$mean, n = 20)(mean),
             popup = ~as.character(paste0(SiteName, "<br>Slope: ", round(mean,3))))

Now we have multiple ways to visualize what is going on with kelps and compare it to the background of change