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

Edzer Pebesma edzer.pebesma at uni-muenster.de
Mon Dec 14 10:13:44 CET 2015



On 11/12/15 15:29, Heuvelink, Gerard wrote:
> Hi Edzer,
> 
> I checked the Ver Hoef and Cressie (1993) paper but could not find a reference to deriving the cross-covariance from the cross-variogram as "C(h)=C(0)-gamma(h) with C(0)=max(0, sum of the positive sill values)". I did notice that their paper is mostly about the pseudo cross-variogram (defined in Eq. 9), while we and presumably gstat as well work with the 'normal' cross-variogram (Eq. 8).
> 
> I think the proper relation is Cij(h)=Cij(0)-Vij(h), where Cij(h)=E[(Zi(s)-E[Zi(s)])*(Zj(s+h)-E[Zj(s+h)])] is the cross-covariance and Vij(h)=0.5*E[(Zi(s)-Zi(s+h))*(Zj(s)-Zj(s+h))] is the normal cross-variogram. This relation is valid if Z is second-order stationary (conditions 1 and 2 on page 220 of VH&C) and satisfies the symmetric cross-covariance condition (i.e. Cij(h)=Cij(-h), condition 3 on the same page). I would recommend that you implement this relation instead of "C(h)=C(0)-gamma(h) with C(0)=max(0, sum of the positive sill values)".
> 

OK, I now changed this in 1.1-1 (r-forge); install from

https://r-forge.r-project.org/R/?group_id=378

> What I still do not understand is that the 'old' version of gstat did not complain about our Allier example, because using the alternative derivation of the cross-covariance from the cross-variogram produces a kriging system that is not positive-definite.
> 

Before 1.1-0 gstat used the LDL' decomposition, which "succeeds on some
indefinite matrices" (wikipedia).

I now kept Choleski the default (as in 1.1-0), and made LDL optional
(when set choleski=0 is passed to gstat()) for pre 1.1-0 compatibility; see

https://r-forge.r-project.org/scm/viewvc.php/pkg/tests/allier.R?view=markup&root=gstat

> Gerard
> 
> -----Original Message-----
> From: R-sig-Geo [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of r-sig-geo-request at r-project.org
> Sent: 18 November 2015 12:00
> To: r-sig-geo at r-project.org
> Subject: R-sig-Geo Digest, Vol 147, Issue 17
> 
> Send R-sig-Geo mailing list submissions to
> 	r-sig-geo at r-project.org
> 
> To subscribe or unsubscribe via the World Wide Web, visit
> 	https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> or, via email, send a message with subject or body 'help' to
> 	r-sig-geo-request at r-project.org
> 
> You can reach the person managing the list at
> 	r-sig-geo-owner at r-project.org
> 
> When replying, please edit your Subject line so it is more specific than "Re: Contents of R-sig-Geo digest..."
> 
> 
> Today's Topics:
> 
>    1. Re: gstat now uses Lapack / failing cokriging (Bruin, Sytze de)
>    2. Re: gstat now uses Lapack / failing cokriging (Edzer Pebesma)
>    3. issues with NAs when using resample function (Beth Forrestel)
>    4. Re: issues with NAs when using resample function (Stephen Stewart)
>    5. grass ascii support in the raster package (Cornel Pop)
> 
> 
> ----------------------------------------------------------------------
> 
> Message: 1
> Date: Tue, 17 Nov 2015 17:14:22 +0000
> From: "Bruin, Sytze de" <sytze.debruin at wur.nl>
> To: "'r-sig-geo at r-project.org'" <r-sig-geo at r-project.org>
> Subject: Re: [R-sig-Geo] gstat now uses Lapack / failing cokriging
> Message-ID: <4d812a6858d942f589a5502e117b1234 at scomp5295.wurnet.nl>
> Content-Type: text/plain; charset="us-ascii"
> 
> Hi Edzer, thanks for your prompt reply. I tried reproducing your matrix r but got a different result, i.e. a valid covariance matrix. The cross-covariances are different. Using the same example as above, my code is as follows:
> 
> 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))
> 
> 
> dists <- spDists(Allier)
> r11 <- variogramLine(g$model$fc, dist_vector = dists, covariance = T)
> r12 <- variogramLine(g$model$fc.por, dist_vector = dists, covariance = T)
> r22 <- variogramLine(g$model$por, dist_vector = dists, covariance = T) r <- cbind(r11, r12) r <- rbind(r, cbind(r12, r22))
> 
>> r
>               [,1]         [,2]         [,3]         [,4]         [,5]         [,6]         [,7]
>  [1,] 0.0041300000 0.0014193874 0.0008959951 0.0005655821 0.0002240733 0.0002700000 0.0008677227  [2,] 0.0014193874 0.0041300000 0.0018397575 0.0013976206 0.0008790828 0.0008677227 0.0002700000  [3,] 0.0008959951 0.0018397575 0.0041300000 0.0020030001 0.0014238096 0.0005477541 0.0011247100  [4,] 0.0005655821 0.0013976206 0.0020030001 0.0041300000 0.0018652970 0.0003457607 0.0008544158  [5,] 0.0002240733 0.0008790828 0.0014238096 0.0018652970 0.0041300000 0.0001369841 0.0005374150  [6,] 0.0002700000 0.0008677227 0.0005477541 0.0003457607 0.0001369841 0.0035700000 0.0013734154  [7,] 0.0008677227 0.0002700000 0.0011247100 0.0008544158 0.0005374150 0.0013734154 0.0035700000  [8,] 0.0005477541 0.0011247100 0.0002700000 0.0012245061 0.0008704261 0.0008669750 0.0017801702  [9,] 0.0003457607 0.0008544158 0.0012245061 0.0002700000 0.0011403233 0.0005472637 0.0013523535 [10,] 0.0001369841 0.0005374150 0.0008704261 0.0011403233 0.0002700000 0.0002168159 0.0008506105
>               [,8]         [,9]        [,10]
>  [1,] 0.0005477541 0.0003457607 0.0001369841  [2,] 0.0011247100 0.0008544158 0.0005374150  [3,] 0.0002700000 0.0012245061 0.0008704261  [4,] 0.0012245061 0.0002700000 0.0011403233  [5,] 0.0008704261 0.0011403233 0.0002700000  [6,] 0.0008669750 0.0005472637 0.0002168159  [7,] 0.0017801702 0.0013523535 0.0008506105  [8,] 0.0035700000 0.0019381256 0.0013776943  [9,] 0.0019381256 0.0035700000 0.0018048825 [10,] 0.0013776943 0.0018048825 0.0035700000
> 
>> eigen(r)$values
>  [1] 0.0124609506 0.0055040132 0.0046425761 0.0035980843 0.0031064016 0.0028989486 0.0028042439  [8] 0.0018107335 0.0010254131 0.0006486352
> 
> 
> 
> ------------------------------
> 
> Message: 2
> Date: Tue, 17 Nov 2015 20:38:33 +0100
> From: Edzer Pebesma <edzer.pebesma at uni-muenster.de>
> To: r-sig-geo at r-project.org
> Cc: "Heuvelink, Gerard" <Gerard.Heuvelink at wur.nl>
> Subject: Re: [R-sig-Geo] gstat now uses Lapack / failing cokriging
> Message-ID: <564B8239.6080401 at uni-muenster.de>
> Content-Type: text/plain; charset="windows-1252"
> 
> 
> On 17/11/15 18:14, Bruin, Sytze de wrote:
>> Hi Edzer, thanks for your prompt reply. I tried reproducing your matrix r but got a different result, i.e. a valid covariance matrix. The cross-covariances are different. Using the same example as above, my code is as follows:
>>
>> 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))
>>
>>
>> dists <- spDists(Allier)
>> r11 <- variogramLine(g$model$fc, dist_vector = dists, covariance = T)
>> r12 <- variogramLine(g$model$fc.por, dist_vector = dists, covariance =
>> T)
>> r22 <- variogramLine(g$model$por, dist_vector = dists, covariance = T) 
>> r <- cbind(r11, r12) r <- rbind(r, cbind(r12, r22))
>>
>>> r
>>               [,1]         [,2]         [,3]         [,4]         [,5]         [,6]         [,7]
>>  [1,] 0.0041300000 0.0014193874 0.0008959951 0.0005655821 0.0002240733
>> 0.0002700000 0.0008677227  [2,] 0.0014193874 0.0041300000 0.0018397575
>> 0.0013976206 0.0008790828 0.0008677227 0.0002700000  [3,] 0.0008959951
>> 0.0018397575 0.0041300000 0.0020030001 0.0014238096 0.0005477541
>> 0.0011247100  [4,] 0.0005655821 0.0013976206 0.0020030001 0.0041300000
>> 0.0018652970 0.0003457607 0.0008544158  [5,] 0.0002240733 0.0008790828
>> 0.0014238096 0.0018652970 0.0041300000 0.0001369841 0.0005374150  [6,]
>> 0.0002700000 0.0008677227 0.0005477541 0.0003457607 0.0001369841
>> 0.0035700000 0.0013734154  [7,] 0.0008677227 0.0002700000 0.0011247100
>> 0.0008544158 0.0005374150 0.0013734154 0.0035700000  [8,] 0.0005477541 0.0011247100 0.0002700000 0.0012245061 0.0008704261 0.0008669750 0.0017801702  [9,] 0.0003457607 0.0008544158 0.0012245061 0.0002700000 0.0011403233 0.0005472637 0.0013523535 [10,] 0.0001369841 0.0005374150 0.0008704261 0.0011403233 0.0002700000 0.0002168159 0.0008506105
>>               [,8]         [,9]        [,10]
>>  [1,] 0.0005477541 0.0003457607 0.0001369841  [2,] 0.0011247100
>> 0.0008544158 0.0005374150  [3,] 0.0002700000 0.0012245061 0.0008704261 
>> [4,] 0.0012245061 0.0002700000 0.0011403233  [5,] 0.0008704261
>> 0.0011403233 0.0002700000  [6,] 0.0008669750 0.0005472637 0.0002168159 
>> [7,] 0.0017801702 0.0013523535 0.0008506105  [8,] 0.0035700000
>> 0.0019381256 0.0013776943  [9,] 0.0019381256 0.0035700000 0.0018048825 
>> [10,] 0.0013776943 0.0018048825 0.0035700000
>>
>>> eigen(r)$values
>>  [1] 0.0124609506 0.0055040132 0.0046425761 0.0035980843 0.0031064016
>> 0.0028989486 0.0028042439  [8] 0.0018107335 0.0010254131 0.0006486352
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
> 
> Building upon your earlier example,
> 
> predict(g, SpatialPoints(data.frame(x=450, y=350)), debug = 32)
> 
> gives you the generalized covariance matrix that is used for the cokriging, which I looked at. gstat computes generalized covariances as
> C(0)-gamma(h) with C(0) = max(0, sum of the positive sill values), instead of the sill of all sill values. If one of the sill components is negative, this matters.
> 
> I looked in the Ver Hoef & Cressie 1993 paper, but couldn't find out which one is right. Maybe Gerard can also take a look at it; the fix would be trivial.
> 
> --
> Edzer Pebesma
> Institute for Geoinformatics (ifgi),  University of M?nster, Heisenbergstra?e 2, 48149 M?nster, Germany; +49 251 83 33081
> Journal of Statistical Software:   http://www.jstatsoft.org/
> Computers & Geosciences:   http://elsevier.com/locate/cageo/
> Spatial Statistics Society http://www.spatialstatistics.info
> 
> -------------- next part --------------
> A non-text attachment was scrubbed...
> Name: signature.asc
> Type: application/pgp-signature
> Size: 490 bytes
> Desc: OpenPGP digital signature
> URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20151117/d7d023d9/attachment-0001.bin>
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi),  University of Münster,
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/
Spatial Statistics Society http://www.spatialstatistics.info

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20151214/453c92a1/attachment.bin>


More information about the R-sig-Geo mailing list