WUR Geoscripting

Extra: Geoscripting in R: demos

This page is intended as a demonstration of geoscripting potential. This is done by providing geo-spatial showcases (e.g. data manipulation, visualization, analysis). Every showcase is followed by:

A reflection of what is the interest of using scripting for that particular problem.
An assessment of the R knowledge required to perform that particular problem.

Note that most of the examples given here are using packages that are now deprecated. They are remains from old tutorials and exercises from this course in the past. Nonetheless, they are still useful for demonstration and inspiration!

Demo 1: Integrated analysis chain from download to analysis

In the script below we do the following:

Download SRTM data for a country e.g. Belgium.
Randomly sample height data.
Visualise results

# Load libraries
library(raster)
library(spatstat)
library(rgeos)

# You can choose your own country here
bel <- raster::getData('alt', country = 'BEL', mask = TRUE) # SRTM 90m height data
belshp <- raster::getData('GADM', country = 'BEL', level = 2) # administrative boundaries

# Create random points
dran <- runifpoint(500, win = as.vector(extent(bel)))

# Make the random point spatial points
S <- SpatialPoints(data.frame(x = dran$x, y = dran$y),
                   proj4string = CRS(proj4string(belshp)))

# Select only the ones within belgium
Sint <- gIntersection(S, belshp)

# Create a map
plot(bel)
plot(belshp, add = TRUE)
plot(
  Sint,
  add = TRUE,
  col = "red",
  pch = 19,
  cex = 0.2
)

Now we can sample the height data, and plot following random point id:

out <- extract(bel, Sint, df = TRUE)
colnames(out) <- c("id", "height")
head(out)

##   id height
## 1  1      8
## 2  2      0
## 3  3     31
## 4  4      0
## 5  5      0
## 6  6     10

plot(out, type = "p", pch = 19)

This example can easily be run for every country of the world with minimal human intervention, and therefore also constitute an other case of automation.

The contribution of scripting

Ability to automatically download data from any place in the world.
Quick and efficient analysis by sampling data randomly.
Can easiliy be repeated for other countries, other datasets (reproducible).
Fully integrated analysis chain, from data download to output.

Basic scripting skills required

String handling (generate automatic file names, recognize patterns, etc …).
Querying lists and data.frames.
Error handling.
Handling raster and vector data.
Understanding and using basic spatial analysis functions.

Demo 2: Visualize the temporal profile of a pixel

The example below illustrates how a small function allows to quickly visualize temporal profile of a raster time-series.

# Download a spatio-temporal dataset
download.file(url = 'https://raw.githubusercontent.com/loicdtx/bfastSpatial/master/data/tura.rda', destfile = 'tura.rda', method = 'auto')
load('tura.rda')

There is now a variable named tura in your environment, it is a time-series of NDVI over an area of Ethiopia.

library(zoo)

# Define the function to extract and plot the time series
click2ts <- function(x) {
  val <- click(x, n = 1)
  z <- getZ(x)
  plot(
    zoo(t(val), z),
    type = 'p',
    pch = 20,
    xlab = 'Time',
    ylab = 'NDVI (-)'
  )
}

Then run the code block below. R will be waiting for you to click on the map, and quickly after the time-series will appear.

plot(tura, 1)
click2ts(tura)

The contribution of scripting

Exploring the temporal dimension of a dataset is often not possible in classical GIS environments. The flexibility of scripting allows to do that with minimal efforts.

Basic scripting skills required

Understanding the structure of spatial objects.
Writing functions.

Demo 3: Implement a custom made function to a raster time-series

The example below uses a time-series of MODIS VCF data, which corresponds to percentage tree cover, to investigate temporal trends of tree cover over the Netherlands.

# Load libraries
library(zoo)
library(lubridate)
library(raster)

# Download data (if it doesn't work, try with method='wget')
download.file(url = 'https://raw.githubusercontent.com/GeoScripting-WUR/Scripting4Geo/gh-pages/data/MODIS_VCF_2000-2010_NL.rds',
              destfile = 'MODIS_VCF_2000-2010_NL.rds',
              method = 'auto')

# Read the data
modis <- readRDS('MODIS_VCF_2000-2010_NL.rds')

# Clean data (values > 100 correspond to water)
modis[modis > 100] <- NA

# Visualize
plot(modis, 1)

# For the example we will reduce the extent so that processing does not take too long
e <- extent(340101, 370323, 5756221, 5787772)
plot(e, add = TRUE)
modis_sub <- crop(modis, e)
plot(modis_sub, 1)

# Define function to calculate temporal trends
fun <- function(x) {
  ts <- zoo(x, time)
  df <- data.frame(t = decimal_date(index(ts)), vcf = c(ts))
  out <- try(lm(vcf ~ t, data = df)$coefficients[2], silent = T)
  if (class(out) == 'try-error')
    out <- NA
  return(out)
}

# Run the function spatially
time <- getZ(modis)
out <- calc(x = modis_sub, fun = fun)

# Visualize output
plot(out)
hist(out, main = 'Tree cover change at 250 m resolution (2000-2010)', xlab = 'Percentage change')

You should see a map displaying the tree cover percentage temporal trends for the 2000-2010 period.

Contribution of scripting

Scripting, thanks to its flexibility, offers the possibility to implement any custom made function to any dataset.

Basic scripting skills required

Writing functions.
Manipulating spatial objects.
Knowing the details of the functions/algorithms to be implemented spatially.

Demo 4: Interactive maps

Check out the nice new leaflet package for R, that creates interactive web maps with the JavaScript ‘Leaflet’ Library.

library(leaflet)
m <- leaflet()
m <- addTiles(m)
m <- addMarkers(m, lng = 5.665349, lat = 51.987870, popup = "Wageningen University")
m

See here for the result of the above code: Wageningen Leaflet Demo

For more information: http://www.htmlwidgets.org/showcase_leaflet.html

The above map is nice, but still pretty static; it only displays two variables, which we did not choose. The year the first park was created for the dot color, and the total number of parks in the country for the dot size. To make that map a bit more interactive we can embed it in a shiny App. Shiny is among the latest developments of the R world, supported by RStudio; it provides the possibility to create interactive plots. This app can either be run in an RStudio session or online.

Advanced examples are possible, as shown here, with an integration of shiny and leaflet. (Example provided by RStudio)

Contribution of scripting

“User friendly” interface to JavaScript libraries.

Basic scipting skills required

Dataframe manipulation.
String manipulation (subsetting, pattern matching).
Function writing.

Reproducibility

Note that all examples shown in this demo have at least one characteristic in common; they are all reproducible. While you may not always remember the parameters you have selected when analyzing or creating map in a GIS environment, scripting leaves a fully reproducible mark of what you have done.

Summary

Geoscripting allows us to:

Automate steps and deal with large amounts of data.
Solve complex problems (run custom made functions, while traditional GIS environments are limited to the features they provide).
Make work reproducible (also beneficial for collaboration).
Visualize data and build interactive geo-data visualization tools.

WUR Geoscripting

Extra: Geoscripting in R: demos

Demo 1: Integrated analysis chain from download to analysis

The contribution of scripting

Basic scripting skills required

Demo 2: Visualize the temporal profile of a pixel

The contribution of scripting

Basic scripting skills required

Demo 3: Implement a custom made function to a raster time-series

Contribution of scripting

Basic scripting skills required

Demo 4: Interactive maps

Contribution of scripting

Basic scipting skills required

Reproducibility

Summary

Inspiration and more info