I'm experimenting with cokriging in Gstat and I'm getting some weird results I cannot explain. Your help will be highly appreciated.
I'm randomly sampling (I keep 10% of the image) an RGB 128x128 (8 bits per band) color image and I'm interpolating it from the samples using simple kriging, simple cokriging, ordinary kriging and ordinary cokriging. Local neighborhoods of radius 5 are used in all the cases. I fitted experimentally simple direct and cross variogram exponential models, assuring that a linear model of corregionalization could be found by gstat.
When I use ordinary kriging I get better results than when I use simple kriging, and when I use ordinary cokriging the results are also better than when I use simple cokriging. Good, this is what I expected but...
I get EXACTLY the same results using kriging and cokriging in both cases (simple and ordinary). I mean that simple kriging and cokriging results are identical and ordinary kriging and cogriging results are identical too. Essentially this means that adding the cross-variogram does not improve the interpolation.
Gstat reports to find the linear model of coregionalization in both cokriging (simple and ordinary) cases. Nevertheless the results are identical to their kriging counterparts.
I can admit that adding channels to the predictor does not improve significantly the result (due to the specific inter- and intra- channel correlation chracteristics), but I cannot understand they're virtually identical (point to point) since they're different algorithms. Moreover, this happens all the time (I've run several times the experiments, changing the sampling location and its density). I've tried to set a debug option to see the kriging weights and some other internal info, but I couldn't manage to understand the log.
What is shocking to me is that there is a significant cross-correlation among channels. As an example you can see below the gst command used for simple cokriging. For simple kriging it is the same, without the cross-variogram lines. For ordinary kriging just the sk_means are removed.
If I'm right, the correlation index function between bands i and j is
(I'm assuming that SEMIvariograms are the inputs to Gstat)
r_ij(h) = (\gamma_ij (0) / (sqrt(\gamma_i(0))*sqrt(\gamma_j(0))) *
exp(-h/30)
and when one puts numbers there is a significant cross correlation specially
between channel 1 and 2: r_12(h) = 0.9081 exp(-h/30)
It is true that inside a neigborhood of radious 5 all the samples you
get in a single channel are highly correlated (eg. r_1(5) = 0.8465), but
in my opinion this is not reason enough to get exactly the same result
when I add the other channels (note the r_12(h) above is quite high).
Now, if I add white gaussian noise to each band (sigma=30) in the RGB
image and I modify the direct and cross variogram definitions in the file
below by adding
900 Err(0) to the direct variograms and 0 Err(0) to the cross-variograms
(in order to make Gstat recognize the linear coregionalization model) I
get a significant improvement when using cokriging over kriging (both for
the ordinary and simple cases. Note that this is not interpolation anymore
but approximation -or smoothing-). Nevertheless, the bad news here are
that now ordinary cokriging is worse than simple cokriging (I've averaged
the MSE after several independent runs), something I cannot understand
either.
I'm running gstat 2.4.0 on Linux
Thanks in advance and best regards.
Juan Ruiz-Alzola
'data(f1): 'gsttemp-filterdata.eas', x=1, y=2, v=3,
radius=5, sk_mean=0;'
'data(f2): 'gsttemp-filterdata.eas', x=1, y=2, v=4,
radius=5, sk_mean=0;'
'data(f3): 'gsttemp-filterdata.eas', x=1, y=2, v=5,
radius=5, sk_mean=0;'
'variogram(f1): 8191.78 Exp(30);'
'variogram(f2): 6003.86 Exp(30);'
'variogram(f1,f2): 6368.51 Exp(30);'
'variogram(f3): 2535.61 Exp(30);'
'variogram(f2,f3): 1133.12 Exp(30);'
'variogram(f1,f3): 1316.44 Exp(30);'
'data(): 'gsttemp-locations.eas', x=1,y=2;'
'set output = 'gsttemp-ok.out';'
-- Juan Ruiz-Alzola, PhD, Assoc. Prof. (1) - (2) - (3) (1) Dep. Signals and Communications, Univ. Las Palmas GC, Spain (2) Research Unit, Dr. Negrin hospital, Las Palmas GC, Spain (3) Dep. of Radiology, Harvard Medical School, Boston, MA, USA Phone/Fax (Univ.): +34 928 452980 / 1243 , (hosp.): +34 928 449289
