[R-sig-Geo] gstat now uses Lapack / failing cokriging

[Bruin, Sytze de]

On 10/20/2015 05:58 AM, Edzer Pebesma wrote:

> gstat 1.1-0, now on CRAN, no longer comes with its own functions for 
> matrix factorization and solving systems of equations [1], but now 
> directly uses Lapack (dpotrf and dtrsm) through R's own lapack interface 
> and R_ext/Lapack.h header files. 

> For global kriging at one location from 10,000 observations, as in 

> library(sp) 
> library(gstat) 
> set.seed(1331) 
> n = 10000 
> pts = SpatialPoints(cbind(x = runif(n), y = runif(n))) 
> pts$z = runif(n) 
> k <- krige(z~1, pts, pts[1,], vgm(1, "Exp", 1)) 

> I see a speed increase from 120 (gstat 1.0-26) to 46 seconds; using 
> openblas on a 4 core laptop brings this down to 15 seconds - I expect 
> sth similar with MKL/RevoR. 

> For local kriging on large data sets with smaller neighborhoods and many 
> locations, I wouldn't expect large improvements; for global kriging of 
> large data sets to many prediction locations, krige0 may be faster when 
> you use openblas or MKL - as long as things fit in memory. 

>  I'd be happy to hear experiences (positive and negative), or otherwise 
> reactions or questions.

Hi Edzer,

With the previous version of gstat the following code for doing cokriging using a linear model of coregionalization worked fine.

library(sp)
library(gstat)

# some data
x <- c(215, 330, 410, 470, 545)
y <- c(230, 310, 330, 340, 365)
fc <- c(0.211, 0.251, 0.281, 0.262, 0.242)
por <- c(0.438, 0.457, 0.419, 0.430, 0.468)
Allier <- data.frame(x, y, fc, por)
coordinates(Allier) = ~x+y

# gstat object for co-kriging using linear model of co-regionalization
g <- gstat(id=c("fc"), formula=fc~1, data=Allier,
           model=vgm(0.00247, "Sph", 480, 0.00166))
g <- gstat(g, id="por", formula=por~1, data=Allier,
           model=vgm(0.00239, "Sph", 480, 0.00118))
g <- gstat(g, id=c("fc", "por"),
           model=vgm(0.00151, "Sph", 480, -0.00124))

# predict at single point
point.krig <- predict(g, SpatialPoints(data.frame(x=450, y=350)))

However, using the recently downloaded version of gstat, I get:

Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
Warning message:
In predict.gstat(g, SpatialPoints(data.frame(x = 450, y = 350))) :
  Warning: Covariance matrix (nearly) singular at location [450,350,0]: skipping...

Yet, the covariance matrix is not nearly singular. This appears to be caused by recent changes to the library.

Regards,

Sytze de Bruin
Wageningen University
Laboratory of Geo-Information Science and Remote Sensing 
Telephone: +31 317 481830
Mobile: +31 6 13898626