First, let’s load our friends: sf
(vector data), raster
(raster data), and dplyr
(general data tidying):
library(sf)
library(raster)
library(dplyr)
Then let’s load the mapview
package, which makes the production of web maps in R very easy:
library(mapview)
For this demonstration, we are using the cookfarm
dataset from the GSIF
package. This dataset is a list with 6 somponents storing various measurements that have been done on the farm.
data('cookfarm', package = "GSIF")
names(cookfarm)
## [1] "readings" "profiles" "bdensity" "grids" "weather"
## [6] "proj4string"
The profiles
component stores soil profile data as a data.frame
. We will turn this into a simple feature using the sf
package:
profiles <- cookfarm$profiles %>%
st_as_sf(
coords = c('Easting', 'Northing'),
crs = cookfarm$proj4string
)
# simple plot using the sf package
plot(profiles)
The grids
component is storing gridded data. We will reansform these into a RasterStack
using the rasterFromXYZ
function. This function requires the x
and y
columns to be the first column of the data.frame
(from the left) – a good opportunity to use the select
function from the dplyr
package:
grids <- cookfarm$grids %>%
select(x, y, DEM, TWI, Cook_fall_ECa, Cook_spr_ECa) %>% # Re-order the columns, and select 4 interesting variables
rasterFromXYZ(crs = cookfarm$proj4string) # The coordinate reference system is stored as one of the component in cookfarm
# simple plot using the raster package
plot(grids)
mapview
The mapview
package is basically a wrapper around the leaflet
package. It makes is very easy to create a leaflet
map. You can basically just use it as you would use plot
:
# couldn't be easier to create a map!
mapview(profiles)
But while the mapview
function can be called just as is, there’s numerous details that can be tweaked. For example the colours:
# Create palettes
library(RColorBrewer)
pal_continuous <- colorRampPalette(brewer.pal(7, "BrBG")) # For continuous data
pal_categorical <- colorRampPalette(brewer.pal(9, "Set1")) # For categorical data
# Pass the palette to mapview
mapview(
profiles,
zcol = "TAXSUSDA",
col.regions = pal_categorical,
legend = TRUE
)
In particular, a wide range of background (web) maps is available. You can pick one or several. You can find a list of these background from the Leaflet-extras project website: http://leaflet-extras.github.io/leaflet-providers/preview/
.
# Tweaking backgrounds
mapview(profiles, zcol = "BLD", col.regions = pal_continuous, legend = TRUE, map.types = "Esri.WorldImagery")
The burst
option can be interesting when you visualise categorical data:
# Burst to separate soil classes
mapview(profiles, zcol = "TAXSUSDA", col.regions = pal_categorical, legend = TRUE, burst = TRUE)
The mapview
option also work for data loaded using the raster
package. If you try to visualise a RasterStack
(as opposed to RasterLayer
), you can select and choose which layer to plot using the selection interface.
# Plot a RasterLayer
mapview(grids$TWI, col.regions = pal_continuous, na.color = "transparent", legend = TRUE)
# Plot a RasterStack
mapview(grids)
# The same sort of options are available
mapview(grids, col.regions = pal_continuous, na.color = "transparent", legend = TRUE)
You can visualise the raster extent — rather than the data itself — using viewExtent
:
viewExtent(grids)
For imagery, it is a little diffferent, because usually one wants to visualise a RGB composite of 3 layers, rather than the individual layers individually. In this case, the viewRGB
is the way to go:
# Imagery specific functions
viewRGB(poppendorf)
# The combination of bands can be changed very easily
viewRGB(poppendorf, 4, 3, 2)
mapviews
togetherConveniently, maps can be aded to each other using the +
operator:
# Create 2 maps
m1 <- mapview(grids$DEM, col.regions = pal_continuous, legend = TRUE)
m2 <- mapview(profiles, zcol = "TAXSUSDA", col.regions = pal_categorical)
# Plot both together
m1 + m2
What mapview
does behind the scenes is to change the projection system to EPSG:3857 (web mercator). In some cases, this is inconvenient and you might actually want to visualise your data in a local projected CRS. plainview
is here to help:
plainview(grids$DEM)
mapview
is a wrapper around leaflet
package, which is a R API for the popular Javascript library for web mapping called Leaflet
. Leaflet
has been designed with simplicity and rapidity in mind. For more power, you’ll have to learn a bit more about the leaflet
package itself. Their website is a great starting point.
A bunch of options can be changed using the mapviewOptions
function:
mapviewOptions(
basemaps = c("Esri.WorldImagery", "Thunderforest.Landscape"),
na.color = "transparent"
)
mapview(profiles, zcol = "BLD") + mapview(grids$DEM)
You can associate and synchronise a set of maps using the sync
function:
# Syncing several maps
m1 <- mapview(grids$DEM)
m2 <- mapview(grids$TWI)
m3 <- mapview(grids$Cook_fall_ECa)
sync(m1, m2, m3, ncol = 2, sync.cursor = TRUE)
This is an interactive analogue to the panelled graphs provided by ggplot2
or lattice
.
Another advanced visualisation tool is slideview
, which is convenient to compare two maps:
img1 <- poppendorf[[1]]
img2 <- poppendorf[[5]]
slideview(
img1,
img2,
label1 = "Poppendorf-Layer-1",
label2 = "Poppendorf-Layer-2",
legend = TRUE
)
The popups can be either a table (popupTable
, default behaviour), an image (popupImage
), or a htmlwidget
(popupGraph
).
# Table
mapview(
profiles,
popup = popupTable(profiles, zcol = 1:2)
)
# Image
mapview(
profiles,
popup = popupImage('https://www.vcard.wur.nl/WebServices/GetMedia.ashx?id=37263')
)
The leaflet
and leaflet.extras
are providing a LOT of different map widgets. The garnishMaps
function facilitates their integration with mapview
:
library(leaflet)
m <- mapview(profiles)
garnishMap(
m,
addMouseCoordinates,
addGraticule,
addScaleBar
)
# Create map
m <- mapview(profiles)
# Save interactive HTML page
mapshot(m, url = 'my_map.html')
# Save static image (PNG, JPEG, or PDF)
mapshot(m, file = 'my_image.png')
mapview
🎉library(xts)
library(dygraphs)
profiles$SOURCEID <- as.character(profiles$SOURCEID)
records <- cookfarm$readings
records$SOURCEID <- as.character(records$SOURCEID)
ids <- unique(records$SOURCEID)
# Subset sensors
ids <- sample(ids, size = 5)
idx_sensors <- which(profiles$SOURCEID %in% ids)
sensors <- profiles[idx_sensors,]
make_ts <- function(id) {
records %>%
filter(SOURCEID == id) %>%
dplyr::select(-SOURCEID) %>%
dplyr::select(Date, ends_with('VW')) %>%
xts(.$Date)
}
make_dygraph <- function(id){
ts <- make_ts(id)
dygraph(ts)
}
l_graphs <- lapply(
ids,
make_dygraph
)
make_dygraph(ids[1])
mapview(sensors, popup = popupGraph(graphs = l_graphs, width = 300, height = 300))