In mjlaine/fastgrid: Gridding point data to a regular grid

This document describes gridding of point observations to a regular grid. The code is mainly intended for interpolating station temperatures (or other meterological variables). The code can use a prior background model field. This text needs some polishing, still.

Examples

Load the fastgrid library and generate a model grid by utility function makelatlongrid. The grid can be either SpatialGrid or SpatialGridDataFrame. The latter can be used if the grid trend model uses some auxiliary variables, such as elevation. In this case these variables have to exist in both the observation data and in the background grid. Here we define a grid over Finland with 0.1° resolution in longitude direction and 0.2° in latitude direction.

library('fastgrid')
lonmin <- 19.0; lonmax <- 33.0; latmin <- 59.0; latmax <- 71.5
londx <- 0.2; latdx <- 0.1
modelgrid <- makelonlatgrid(lonmin=lonmin,lonmax=lonmax,latmin=latmin,latmax=latmax,
                            londx=londx,latdx=latdx)

Generate some random points in the defined region as a SpatialPointsDataFrame to be used as observation to be gridded. The names for the coordinates are longitude and latitude.

nobs <- 50
obs <- data.frame(longitude=runif(nobs,lonmin,lonmax),
                  latitude=runif(nobs,latmin,latmax),
                  temperature=rnorm(nobs,mean=20,sd=2))
sp::coordinates(obs) <- c('longitude','latitude')

Interpolate to model grid by Kriging. The function pointgridding is in the package fastgrid. The default trend model given here fits the mean field. It could use coordinate information and auxiliary variables defined in both observations and grid.

cov.pars <- c(2.0^2,1.5,0.1) # sigma2, correlation length, nugget
trend_model <- temperature ~ 1
t0<-proc.time()
obs.gridded <- pointgridding(pointdata=obs,modelgrid=modelgrid,
                             trend_model=trend_model,
                             cov.pars=cov.pars,
                             variable="temperature")
print(t1<-proc.time()-t0)

The output will be a sp::SpatialGridDataFrame. It can be plotted by sspplot.

library(sp)
#library(maps)
#library(maptools)
#library(ggplot2)

spplot(obs.gridded)

In another package MOSplotting there is a function MOS_plot_field.

MOSplotting::MOS_plot_field(obs.gridded,stations=obs,cmin=10,cmax=30,main='kriging interpolation')

Mathematical background

The gridding by Kriging can be described as an hierarchical linear model: \begin{equation} \label{equ:hiermodel} \begin{split} y|\eta &\sim N(H\eta,\Sigma_y), \qquad\text{observations}\ \eta|\beta &\sim N(X\beta,\Sigma_\eta), \qquad\text{spatial random field}\ \beta &\sim N(\beta_0,\Sigma_\beta) \qquad\text{hyper parameters} \end{split} \end{equation}

Here, $\eta$ is the model grid as a vector and described as a spatial random field to be estimated. Its correlation structure is defined by a spatial covariance matrix $\Sigma_\eta$ (called $B$ as the background covariance) in the code). The covariance is defined by a covariance function, which is assumed to have exponential form depending on the distance of points.

The Kriging prediction equation gives estimates for the variable $y$ at each grid location,
\begin{equation} y_\mathrm{grid} = x_\mathrm{grid}\widehat{\beta} + b \Sigma_\eta^{-1} (y - X\widehat{\beta}), \end{equation} where $b$ is a row vector of correlations between the grid location and all the observation locations and $\widehat{\beta}$ is the standard least squares estimate of ${\beta}$. The code does not calculate prediction uncertainty.

Matrix $H$ is the observation operator, which maps the model grid to observation locations. A bilinear interpolation is assumed in the code. The $\beta$ parameters are estimated using a trend model. The default trend model fits only intercept, i.e. the observation mean is returned for locations with no observation near. Currently, uniform prior for the trend model parameters is assumed. A prior mean field can be used by setting trend_model = y~-1 and substracting the mean from the observation before the fit obs <- obs - H%*%bg and then addind the prior field back to the results after fit.

Code details

The essential parts of the fortran code is the following

  call chol(b, nobs)               ! replace B with L = chol(B,’lower')
  call ldivide(b,x, nobs, npar)    ! replace x with inv(L)*x
  call ldivide(b,y, nobs, 1)       ! replace y with inv(L)*y
  call lsq(x,y,beta, nobs, npar)   ! solve x from y=x*beta
  call resid(x,y,beta, nobs, npar) ! y = y - x*beta
  call ltdivide(b,y, nobs, 1)      ! y = inv(L')*y
  ! calculate prediction over grid
  do i=1, ngrid
     x1 = grid(i,:)
     call corrdist(cgrid(i,:),cy,b1, nobs, covpars)
     ypred(i) = matmul(x1,beta) + matmul(b1,y)
  end do