Creating Interactive Visualisation Tools With R

More interaction - the ggvis package

library(ggvis)
load('Rainfall.RData')
rain %>% group_by(Year) %>%
  summarise(mr=mean(Rainfall)) %>% 
  ggvis(x=~Year,y=~mr) %>% layer_lines()

Storing plots to show later

rain %>% group_by(Year) %>%
  summarise(mr=mean(Rainfall)) %>% 
  ggvis(x=~Year,y=~mr) -> rain_ggv
  • Nothing plotted yet
  • But information for plot creation is in rain_ggv
  • Can add to this later on…

Employing a stored plot

rain_ggv %>% layer_lines(stroke := 'firebrick')
  • Note := for definition
    • ordinary = selects factor
  • ggvis plots are resizable

Re-employing a stored plot

rain_ggv %>% layer_points(fill := 'navy') %>% 
  layer_lines(stroke := 'firebrick')

Smoothing in ggvis

rain_ggv %>% layer_smooths(span=1,stroke := 'dodgerblue')
  • span is the degree of smoothing
  • span=1 covers whole dataset, but is ‘tapered’

Smooth plus points

rain_ggv %>% layer_smooths(span=1,stroke := 'dodgerblue') %>% layer_points(fill := 'firebrick')
  • Demonstrates signal vs. noise comparison

Interactive smoothing

rain_ggv %>% layer_smooths(span=input_slider(0.1,1.5),stroke := 'dodgerblue') 
  %>% layer_points(fill := 'firebrick')

  • As yet, ggvis can’t do this in a simple presentation
  • But you can in RStudio

Interactive smoothing 2

rain_ggv %>% layer_points(fill := 'firebrick') %>% 
  layer_smooths(span=waggle(0.1,1.5),stroke := 'dodgerblue')

  • As before, can’t yet embed in a presentation
  • waggle should be last item in pipeline
    • When R sees it, it starts waggling and ignores remainder

Neater Labelling

rain_ggv %>% layer_smooths(span=1,stroke := 'dodgerblue') %>% layer_points(fill := 'firebrick') %>%
  add_axis(type='x',format='4d',subdivide=10) %>% add_axis(type='y',title='Rainfall (cm per year)')

How format works in ggvis

  • Quite complex, but
    • 4d,5d etc means decimal occupying 4, 5 spaces
    • 4,d etc as above with commas after 1000’s
    • 6.2f display decimals to 2 places
    • 03d pad with leading zeroes - ie 007
    • See here for full information

Spatial Interpolation

Basic Idea

  • We have some observations at a number of irregular points
    • Each point has
      • a geographical location
      • An observed measurement (eg Rainfall)
    • Task is to estimate measurement levels at other locations
    • Basically to ‘fill in the gaps’

Method 1 - Voronoi Tesselation

  • Sometimes also called ‘Thiessen Polygons’
    • Each polygon is nearer to one point than the others
    • Assign point value to whole polygon

How to do this in R - 1

  • A lot of this is down to data manipulation
    • Create the geographical point data set
    • Create the Thiessen polygons
    • Join the point date onto the Theissen polygons
    • Map the result
# Step 1 - create a January rainfall point data set
# Create the January rainfall data
ms_rain <- rain %>% 
  group_by(Station,Month) %>% 
  summarise(rain = mean(Rainfall))
ja_rain <- ms_rain %>% filter(Month == 'Jan')

# Create an sf point object for the stations, then join the january averages
stat_dat <- st_as_sf(stations,coords=c('Long','Lat'),crs=4326) %>% 
  st_transform(29902) %>% 
  left_join(ja_rain)

How to do this in R - Steps 2 + 3

# First off create the voronoi tesselation (Theissen polys) - Step 2
vor <- stat_dat %>% 
  st_union() %>% 
  st_voronoi() %>% 
  st_collection_extract() %>% 
  st_sf() 

# Next join back the January rainfall data
vor_jan <- vor %>% st_join(stat_dat)

How to do this in R - Clip to Irish coastline

coast <- st_read('ireland_coast.geojson',quiet=TRUE) # Stops printing info
vor_jan <- vor_jan %>% st_intersection(coast)
tm_shape(vor_jan) + tm_borders() + tm_shape(stat_dat) + tm_dots(size=0.3)

How to do this in R - The map

tm_shape(vor_jan) + tm_polygons(col='rain') + tm_legend(legend.outside=TRUE)

Field-based Smoothing

  • A different approach
    • Models the observations as part of a field
      • That is, a continous ‘surface’ over geographical space
      • Typically detects trends

Interpolation Vs Trend Smoothing

  • Interpolation (1d) \(\Leftrightarrow\) Thiessen/Voronoi Polygons (2d)
  • Trend Smoothing (1d) \(\Leftrightarrow\) Field modelling (2d)

How this works

  • Fit a trend (field) model to the data
    • eg Thin plate spline, Kriging
    • Evaluate on a raster (pixel-based) grid
    • Map it

How to do it: 1. Fit a thin-plate spline

library(fields) # You may need to install this package
tps_model <- Tps(st_coordinates(stat_dat), stat_dat$rain,method='REML')
tps_model
## Call:
## Tps(x = st_coordinates(stat_dat), Y = stat_dat$rain, method = "REML")
##                                               
##  Number of Observations:                25    
##  Number of parameters in the null space 3     
##  Parameters for fixed spatial drift     3     
##  Model degrees of freedom:              9.3   
##  Residual degrees of freedom:           15.7  
##  GCV estimate for sigma:                19.23 
##  MLE for sigma:                         18.04 
##  MLE for rho:                           63790 
##  lambda                                 0.0051
##  User supplied rho                      NA    
##  User supplied sigma^2                  NA    
## Summary of estimates: 
##                  lambda       trA      GCV      shat -lnLike Prof converge
## GCV        3.059061e-05 23.750015 209.1276  3.233613     101.4826       NA
## GCV.model            NA        NA       NA        NA           NA       NA
## GCV.one    3.059061e-05 23.750015 209.1276  3.233613           NA       NA
## RMSE                 NA        NA       NA        NA           NA       NA
## pure error           NA        NA       NA        NA           NA       NA
## REML       5.101444e-03  9.295568 588.6009 19.228773     101.2859        3

Raster evaluation

  • raster a package to handle pixel-based geographical objects
    • extent a bounding box for the pixels
    • resolution the size of the pixels (here 1km)
    • interpolate fits the field model to the center of each pixel
    • We need to set the coordinate system for the raster - crs
      • it works differently from sf
    • Use the coast as a mask - block out pixels not on the land mass
      • rasterize turns the coast sf object into a raster
library(raster) # You may need to install this
rex <- extent(st_bbox(coast))
cover_grid <- raster(rex,resolution=1000)
rain_grid <- interpolate(cover_grid,tps_model)
crs(rain_grid) <- CRS("+init=epsg:29902")
m <- rasterize(coast,rain_grid)
rain_grid <- mask(rain_grid,m)

Map of Field

tm_shape(rain_grid) +
  tm_raster(style='cont',title='Monthly Rainfall (cc)') + 
  tm_shape(coast) + 
  tm_borders()

Conclusion

💡 New ideas

  • New R ideas

    • raster objects

    • Thiessen / Voronoi Polygons

    • Interpolation

    • Smooth Trend Modelling/Fields

    • ggviz and interactive exploration

    • Next lecture - Reproducible research