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)".
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.
Gerard
-Original Message-
From: R-sig-Geo [mailto:r-sig-geo-boun...@r-project.org] On Behalf Of
r-sig-geo-requ...@r-project.org
Sent: 18 November 2015 12:00
To: r-sig-geo@r-project.org
Subject: R-sig-Geo Digest, Vol 147, Issue 17
Send R-sig-Geo mailing list submissions to
r-sig-geo@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-requ...@r-project.org
You can reach the person managing the list at
r-sig-geo-ow...@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 +
From: "Bruin, Sytze de"
To: "'r-sig-geo@r-project.org'"
Subject: Re: [R-sig-Geo] gstat now uses Lapack / failing cokriging
Message-ID: <4d812a6858d942f589a5502e117b1...@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.004130 0.0014193874 0.0008959951 0.0005655821 0.0002240733
0.000270 0.0008677227 [2,] 0.0014193874 0.004130 0.0018397575
0.0013976206 0.0008790828 0.0008677227 0.000270 [3,] 0.0008959951
0.0018397575 0.004130 0.0020030001 0.0014238096 0.0005477541 0.0011247100
[4,] 0.0005655821 0.0013976206 0.0020030001 0.004130 0.0018652970
0.0003457607 0.0008544158 [5,] 0.0002240733 0.0008790828 0.0014238096
0.0018652970 0.004130 0.0001369841 0.0005374150 [6,] 0.000270
0.0008677227 0.0005477541 0.0003457607 0.0001369841 0.003570 0.0013734154
[7,] 0.0008677227 0.000270 0.0011247100 0.0008544158 0.0005374150
0.0013734154 0.003570 [8,] 0.0005477541 0.0011247100 0.000270
0.0012245061 0.0008704261 0.0008669750 0.0017801702 [9,] 0.0003457607
0.0008544158 0.0012245061 0.000270 0.0011403233 0.0005472637 0.0013523535
[10,] 0.0001369841 0.0005374150 0.0008704261 0.0011403233 0.000270
0.0002168159 0.0008506105
[,8] [,9][,10]
[1,] 0.0005477541 0.0003457607 0.00013698
