SPEI example

SPEI.R

library(raster)
library(rasterVis)
library(zoo)

## Download SPEI file
setwd(tempdir())
download.file('https://digital.csic.es/bitstream/10261/104742/5/SPEI_03.nc', 'SPEI_03.nc', method = 'wget')

## Read it with `raster`
SPEI <- brick('SPEI_03.nc')
## Use better names for layers with `yearmon`
names(SPEI) <- as.yearmon(getZ(SPEI))

## Define some toy points
lats <- seq(40, 43, by = 1)
lons <- seq(-5, -1, by = 1)
ll <- expand.grid(lons, lats)
myPoints <- SpatialPoints(ll, proj4string = CRS('+proj=longlat +datum=WGS84'))
## Extract values from the `RasterBrick` at these points
vals <- extract(SPEI, myPoints)
## Convert the result into a time series using `zoo`
## First, extract the time index
tt <- getZ(SPEI)
## Second, transpose the `vals` matrix, and assign the time index
vals <- zoo(t(vals), tt)
## Finally, column names.
names(vals) <- paste0('P', seq_len(nrow(ll)))
## Display the result with one panel for each point
xyplot(vals)
## Or superpose them in a unique panel
xyplot(vals, superpose = TRUE,
       lwd = .3, col = 'black', alpha = 0.5,
       auto.key = FALSE)

## Crop the raster using a smaller extent
SPEIc <- crop(SPEI, extent(-10, 3, 30, 50))
## Display some layers using `rasterVis::levelplot`
levelplot(SPEIc, layers = 1000:1011, par.settings = RdBuTheme)
## Superpose the points using `latticeExtra::layer`
levelplot(SPEIc, layers = 1000:1011, par.settings = RdBuTheme) +
    layer(sp.points(myPoints, pch = 19, cex = .5, col = 'black'))