[R-sig-Geo] Spatial interpolation using thin plate spline [fields] through time

[Robert Way]

Hello,
I am trying to interpolate air temperature surfaces through time and space
in my study region. After testing a number of different techniques (e.g.
co-kriging) I have evaluated that for my test data the Thin Plate Spline
implementation from the *fields *package gives the best results against
out-of-sample data over a test period.

I am looking to apply the TPS interpolation to cover the whole time period
of study (e.g. 50 years) at the seasonal timescale with .geotiffs of
interpolated air temperatures being output at each step meaning that I am
looking to iteratively feed the TPS function the same X-inputs
(xposition,yposition, zposition) from the station data but different
columns of my Y-input with the temperature of each season through the time
interval.

I am hoping to get some advice on the best way of solving this problem
which efficiently uses computing resources considering this project will be
a huge drain on computing time for the more powerful computers in our Lab.
For instance, is it best if I find a way to build a function and feed it
with a loop or alternatively other methods such as with apply(). Another
question is whether it is best to output each raster individually or to
create a rasterstack or rasterbrick and split it afterwards. Overall i've
found outputting the .tif files and running calcs on large rasters in R to
be time consuming so any advice would be helpful. Thank you for your time.


Although I have not written a test example the code below is very simple
and shows what I have done:

### Load required packages ###
library(fields)
library(raster)
library(spatial.tools)
library(gdalUtils)
library(mapplots)
library(rgdal)
library(gstat)

### Use more computing resources ###
sfQuickInit()

### Load X-Data - in this case three columns with x,y and z positions per
station ###
xdata <- file.choose()
x <- read.csv(xdata, TRUE, ",")

### Load Y-Data - in this case data.frame of air temperatures through time
per station ###
### Columns = time intervals ### Rows = individual stations
ydata <- file.choose()
y <- read.csv(ydata, TRUE, ",")

### Load elevation data ###
elev1k <- raster("dem1k.tif")

### Thin plate spline fit for first time interval ###
fit <- Tps(x, y[,1])

## Output diagnostic plots ###
set.panel(2,2)
plot(fit)
dev.off()

### Creating interpolated surface for first time interval ###
result <- interpolate(elev1k,fit, xyOnly=FALSE)

### Output geotiff ###
writeRaster(result, "TPSresult.tif", overwrite=TRUE)

-- 

*Robert*

	[[alternative HTML version deleted]]