Temporal data processing/GRASS R raster time series processing

GRASS-R / R-GRASS for raster time series processing

This little example will guide you through the steps to export a Spatio-Temporal Raster Dataset (strds) stored in GRASS, import it into R, prepare the data properly to use the Data INterpolation Empirical Orthogonal Functions algorithm (DINEOF) and, after running it, rebuild your raster time series, export it and import the new strds into GRASS.

We will use North Carolina climatic data already available as a GRASS location. You can get it here. If not done yet, you need to unzip the file and paste into your GRASS database folder (usually named grassdata).

Shall we start?

1. Open GRASS in the location/mapset of interest

grass70 $HOME/grassdata/nc_climate_spm_2000_2012/climate_1970_2012/

• List the available raster maps

g.list type=raster
g.list type=raster pattern="*tempmean"

• Create strds with mean temperature maps

t.create output=tempmean type=strds temporaltype=absolute \
title="Average temperature" \
description="Monthly temperature average in NC [deg C]"

• Register maps in the strds, list the strds in the mapset and obtain general info for a particular space-time dataset

t.register -i input=tempmean type=raster start=2000-01-01 \
increment="1 months" \
maps=`g.list type=raster pattern="*tempmean" separator=comma` 

t.list type=strds

t.info input=tempmean

• List maps and check basic statistics of the strds

t.rast.list input=tempmean

t.rast.univar input=tempmean

• Export out of GRASS

As said before, We will use DINEOF as gap-filling technique. This method is not (yet) available in GRASS GIS, so we need to export the raster time series.

# set default region
g.region -d

t.rast.export input=tempmean output=tempmean4R.tar.gzip compression=gzip

2. Moving to R

From the GRASS console, open R (or rstudio)

rstudio $HOME/Documents/foss4g_bonn/ts_grass_r &

Note that you can directly open a project in which you were already working (as in the example above) or just type "R" or "rstudio" followed by "&" so you do not loose the command prompt in GRASS.

Already in R, we first load "rgrass7", the library that provides the interface with GRASS GIS 7, check some basic information about the session and list maps.

library(rgrass7)

# check grass environment
genv <- gmeta()
genv

# list maps
execGRASS("g.list", parameters = list(type = "raster", pattern="*temp*"))

Given that neither GRASS temporal modules nor GRASS space time data sets are enabled in R, when you need to process a time series from GRASS you need to export maps and then read them into R. For that, we need to load the following packages, too:

library(spacetime)
library(raster)
library(rgdal)

Now, we read the strds exported from GRASS into R

tempmean_in_R <- read.tgrass("tempmean4R.tar.gzip", localName = FALSE, useTempDir = FALSE)

Note that R imports our exported STRDS as a RasterStack.

Let us see if all is there:

# get information of the RasterStack
class(tempmean_in_R)
dim(tempmean_in_R)

# summary of the RasterStack
summary(tempmean_in_R)
# summary of the first layer
summary(tempmean_in_R[[1]])

# plot 3rd layer in the stack
plot(tempmean_in_R,3)

For most methods in R (theil-sen slope, fft, dineof), we will need to transform our RasterStack into a matrix or data.frame. In the particular case of this example, we need to transform our set of rasters into a matrix mxt, with m=maps and t=time, i.e.: in this matrix we have our maps as rows and time in columns. Therefore, in each column we have the time series for a given pixel.

Let us start, then.

# rasterstack to matrix --> this gives us maps as columns and time in rows
txm_tempmean <- as.matrix(tempmean_in_R) 
dim(txm_tempmean)

# transpose matrix --> we need time in columns and maps unfolded in rows
mxt_tempmean <- t(txm_tempmean)
dim(mxt_tempmean)
str(mxt_tempmean)

All the point of DINEOF is gap-filling, but as our sample dataset comes from interpolations of weather stations, we do not have such gaps... So, we will create them.

# Check for gaps in data
sum(is.na(mxt_tempmean))

# Create some holes. 
# Total pixel counts: 14456*156=2255136
# NULL values must appear as NaN
set.seed(46)
n=400000
mxt_tempmean[mysamples<-sample(length(mxt_tempmean), n)]<-NaN
sum(is.na(mxt_tempmean))

Now we are ready, to use DINEOF to fill our gappy data. You can find more info about DINEOF here:

LINKS!!!

and the package is here: https://github.com/marchtaylor/sinkr

1. Load library and other required libraries

library(sinkr) library(irlba) library(Matrix)

1. or, download the dineof function and source it:
2. source('~/Documents/foss4g_bonn/ts_grass_r/dineof.r')

1. Run the algorithm - default settings

result_tempmean <- dineof(mxt_tempmean)

The result_tempmean object is a list with all results. You may investigate yourself, try with other settings and compare.

1. Explore what's inside

names(result_tempmean)

For the purpose of this example, we will only extract the reconstructed spatio-temporal matrix and display it:

1. Extract gap-filled matrix

tempmean_new <- result_tempmean$Xa

1. Set color palette

library(RColorBrewer) pal <- colorRampPalette(c("blue", "cyan", "yellow", "red"))

1. Display the gappy and gap-filled spatio-temporal matrixes

par(mfrow=c(2,1)) image(mxt_tempmean, col=pal(100), main="Original data with gaps") image(tempmean_new, col=pal(100), main="Gap-filled data")

1. Plot the time series of a single pixel from both matrixes

pixel<-10000 ylim<-range(c(mxt_tempmean[,pixel], tempmean_new[,pixel]),na.rm=TRUE) plot(mxt_tempmean[,pixel], t="l", ylim=ylim, lwd=2, ylab="Temperature", xlab="Time") lines(tempmean_new[,pixel], col=2, lty=3) legend(123,13,c("original","gap-filled"),lty = c(1,3),

      lwd = c(2,1),col = c(1,2), cex=0.8)

abline(h=0, v=0, col="gray60")

Let us assume that we are happy with the result, we need now to go back to a STRDS format. We need to rebuild the raster time series starting from a matrix mxt. We first create empty rasters of the same dimensions of our original data and then, we fill them with each row of the gap-filled matrix. What we obtain is a list of rasters.

1. Create raster objects to fill with data from Xa

tempmean_2_grass <- raster(nrows=104, ncols=139,

           xmn=624500, xmx=694000, 
           ymn=208500, ymx=260500, 
           crs=tempmean_in_R)

dim(tempmean_2_grass)

1. Extract data from each row of Xa and build up a raster layer
2. Output: list of rasters

tempmean_new_rl <- lapply(1:nrow(tempmean_new), function(i) {

 setValues(tempmean_2_grass, tempmean_new[i,])

} ) length(tempmean_new_rl) dim(tempmean_new_rl3)

1. Plot a raster layer in the list

plot(tempmean_new_rl3)

1. Compare with layer 3 in the original data

plot(tempmean_in_R,3)

So far we have a list of rasters. But to export them and import them back into GRASS, we need again a RasterStack.

1. rebuild RasterStack

tempmean_new_rs <- stack(tempmean_new_rl) class(tempmean_new_rs)

1. we need to add time to the RasterStack

time_4_new_rs <- seq(as.Date("2001-01-01"),as.Date("2012-12-01")) tempmean_new_rs_and_time <- setZ(tempmean_new_rs,time_4_new_rs,name="time")

1. export RasterStack to GRASS with write.tgrass

write.tgrass(tempmean_new_rs, "tempmean_new_from_R.tar.gzip")

Now, switch back to GRASS console and import the strds with t.rast.import. This command will create a strds and register all maps in it.

t.rast.import input=tempmean_new_from_R.tar.gz output=tempmean_dineof base=tempmean_dineof extrdir=/tmp

Enjoy!