Mapping examples in Nextjournal

This is a reproduction on the Nextjournal platform of two R Markdown notebooks that originally were created for a Nature Toolbox article by Jeffrey Perkel, published 5 June 2018. The following maps are created using the Leaflet R plugin. The Leaflet library produces interactive maps that allow users to pan, zoom, and click on points of interest.

To modify and run this document, create an account on Nextjournal, then open this notebook and click the 'remix' button to create your own working copy. (Executable versions of this document are also available on Code Ocean and Binder.)

This Nextjournal pre-packaged environment should contain all the libraries we need to get started. Let's check...

pkgs = installed.packages()
mypkgs = c("sf","leaflet","rgdal","raster","geojsonio","netcdf","ncdf4","fields")
for (pkg in mypkgs) {
  print(paste(pkg, pkg %in% pkgs))
}

0.6s

Load libraries...

library (sf)
library (leaflet)
library (rgdal)
library (raster)
library (geojsonio)
library(ncdf4)
library(fields)
library(htmlwidgets)

# to display the Leaflet maps
options(nextjournal.display.htmlwidgetsEnabled=TRUE)

0.4s

Research sites in London

I created and read in a table of latitude/longitude values in comma-separated values (CSV) format, using values extracted from Google Maps. Then I plot those points on a topographical map of London, and overlay a second data source from Macrostrat.org, which colors the map according to its geological features.

df <- read.csv(textConnection(
"Name,Lat,Long
Nature,51.533925,-0.121553
Francis Crick Institute,51.531877,-0.128767
University College London,51.524486,-0.133997
MRC Laboratory for Molecular Cell Biology,51.524435,-0.132495
King's College London,51.511573,-0.116083
Imperial College London,51.498780,-0.174888
Cambridge University,52.206960,0.115034
Oxford University,51.754843,-1.254302
Platform 9 3/4,51.532349,-0.123806
"))

0.4s

# create a blank map
m <- leaflet() 

# mark points with blue circles, except Nature's offices in red...
m <- addCircleMarkers(m, lat=df$Lat, lng=df$Long, label=df$Name, 
               color = (ifelse (df$Name == "Nature", "red", "blue"))) 

# add an underlying basemap
m <- addProviderTiles(m, providers$OpenTopoMap) 

# center the view on London
m <- setView(m, -0.119,51.525, zoom = 8) 

# pull in MacroStrat tiles
m <- addTiles (m, 'https://tiles.macrostrat.org/carto/{z}/{x}/{y}.png', 
               options = tileOptions(opacity = 0.6))

# draw the map
m

0.7s

Hurricane Irma

In this example, we will plot the US National Hurricane Center forecast track for Atlantic hurricane Irma using publicly available data (advisory #20, 0900 GMT on 04 September 2017), overlaid with potential coastal targets. Then we will overlay sea surface temperature data for the same date, using code provided by Luke Miller at San Jose State University.

Note that Nextjournal allows you to include data files in the document itself.

al112017_5day_020.zip

First, let's create a directory to hold these shape files. Note that this is a bash cell.

mkdir -p al112017_5day_020

0.1s

Bash in R

Store the files, using R's unzip function.

unzip(al112017_5day_020.zip
, exdir="al112017_5day_020")

1.2s

Next, read in data. These latitude and longitude data were obtained from Google Maps.

df <- read.csv(textConnection(
"Name,Lat,Long
Barbuda,17.648,-61.779
Sint Maarten,18.060,-63.049
British Virgin Islands,18.349,-64.718
Santo Domingo,18.483,-69.929
San Juan,18.421,-66.059
Havana,23.111,-82.357
Key West,24.555, -81.779
Miami Beach,25.790,-80.135
New Orleans,29.999,-90.087
Houston,29.782,-95.405
"))

# Import the advisory 'shape' files that we extracted above...
irmaPolygon <- readOGR(dsn="al112017_5day_020", layer="al112017-020_5day_pgn")
irmaLine <- readOGR(dsn="al112017_5day_020", layer="al112017-020_5day_lin")
irmaPoints <- readOGR(dsn="al112017_5day_020", layer="al112017-020_5day_pts")

2.2s

Plot the map using Leaflet. Color the coastal cities in blue, the hurricane in red, and give the polygon a black border. Add text to label the figure. And center it on Havana, Cuba.

# use the DATELBL field in the irmaPoints dataset to label the map
labels <- as.data.frame(irmaPoints)$DATELBL

# create a blank map
m <- leaflet () 

# add the basemap
m <- addProviderTiles(m, providers$OpenTopoMap)  

# draw the 'cone of uncertainty'
m <- addPolygons(m, data = irmaPolygon, color = "black", weight = "2", fillColor="red", fillOpacity = 0.3) 

# draw the predicted track
m <- addPolylines(m, data = irmaLine) 

# draw the predicted track positions
m <- addCircleMarkers(m, data = irmaPoints, color = "red", radius = 2, label = labels, labelOptions = labelOptions(noHide = TRUE, textOnly = TRUE)) 

# add the coastal locations
m <- addCircleMarkers(m, lat=df$Lat,lng=df$Long,label=df$Name, radius = 2, 
labelOptions = labelOptions(noHide = TRUE, textOnly = TRUE)) 

# center the map on Havana, Cuba
m <- setView(m, -82.357,23.111,zoom=4)

# display the map
m

0.7s

Let's overlay SST data for the same day, 4 September 2017. This code comes from Luke Miller at San Jose State University in California and is based on a 2014 blog post. It retrieves sea surface temperatures for a given latitude/longitude and date.

This code requires the R script, NOAA_OISST_ncdf4.R, from https://github.com/millerlp/Misc_R_scripts/blob/master/NOAA_OISST_ncdf4.R. Download the OISST data themselves and the file lsmask.oisst.v2.nc from https://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.highres.html

NOAA_OISST_ncdf4.R

sst.day.mean.2017.nc

lsmask.oisst.v2.nc

# Basic OISST data load
source (lsmask.oisst.v2.nc
) 

ssts = extractOISSTdaily(lsmask.oisst.v2.nc
,lsmask.oisst.v2.nc
, lonW=160, lonE=359, latS=-10, latN=60, date1='2017-09-04', date2='2017-09-05')

s = brick(ssts, 
          xmn = as.numeric(attr(ssts, 'dimnames')$Long[1]),
          xmx = as.numeric(attr(ssts, 'dimnames')$Long[ncol(ssts)]),
          ymn = as.numeric(attr(ssts, 'dimnames')$Lat[nrow(ssts)]),
          ymx = as.numeric(attr(ssts, 'dimnames')$Lat[1]),
          crs = '+proj=longlat +datum=WGS84')
s = dropLayer(s, 2)

# Necessary steps to get leaflet to plot western hemi data
# See https://github.com/r-spatial/mapview/issues/6
b = shift(s, -360)
SST <- trim (crop(b, extent(-100,180,-90,90),snap='in'))
m <- addRasterImage (m, SST, colors = "Set3") 
m <- setView(m, -64.768, 32.295, zoom = 4)
m

33.2s

Note that the raster package can also plot the OISST data directly, and provide a nice legend...

plot(SST)

0.9s

Computational reproducibility

Collect session information, for computational reproducibility.

sessionInfo()

0.6s