[Gstat-info] LOCAL universal kriging with GLOBAL reg. coeff.estimation

A late follow up question to this thread....

What if the local neighbor hood was restricted to the range of 
autocorrelation? My impression was that values beyond that have little 
weight in the interpolation anyways.

I realize that this is of course not statistically optimal, but for 
those with large data sets (and lacking the computing power) would 
this not give a reasonable approximation?

Dear Radim, Edzer,

I was thinking about the same problem few years ago (I assume that you 
work with auxiliary maps and
not only coordinates).

I think (have a feeling) that local and global Universal kriging 
should be treated as two things
(especially if you put a very narrow search radius). This is because, 
in the case of KED algorithm,
both regression and kriging part of predictions are solved at the same 
time. Hence, if we limit the
search window, but keep a constant variogram model, we could obtain 
very different predictions then
if we use a global (regression-kriging) model. Only if the variogram 
of residuals if absolutely
stationary (constant), then we can limit the search window to fit the 
KED weights. In fact, I am
confident that software packages should not allow UK/KED predictions 
with a limited search radius,
because this conflicts the model assumptions (a bit deeper discussion 
about this problem can be
found in my lecture notes section 2.2 "Localized versus local models").

In the case of regression-kriging, this is less problematic because it 
is computationally cheap to
fit the regression part, and then you are allowed to limit the search 
radius to interpolate the
residuals (we still cheat but not so much). But then the problem is 
that we do not have an estimate
of the RK kriging variance (you could estimate the OLS prediction 
variance and then OK-residuals
prediction variance and then add them together - both can be done in 
gstat - but then you completely
ignore correlation of the two terms; although this is often 
insignificant).

If gstat already does something like this (global GLS simulations of 
trend with local kriging), then
this would be a legitime solution. But I still can not run much 
geostatistics with big datasets
(>>1000 points, >>10^6 grids) in R+gstat (on standard PCs), not to 
mention geostatistical
simulations (that are computationally much more demanding than 
predictions)... if you would have to
do it 100 times, this would really take a lot of processing.

One alternative is to run UK/KED in SAGA software (I have tested it - 
it runs about 5-10 times
faster than R+gstat and has no memory limit problems), which is 
possible through R also, e.g.:

rsaga.get.usage("geostatistics_kriging", 3)
rsaga.geoprocessor(lib="geostatistics_kriging", module=3,
param=list(GRIDS="DEM.sgrd;SLOPE.sgrd;PLANC.sgrd;TWI.sgrd;SINS.sgrd;SMU1.sgrd;SMU3.sgrd;SMU4.sgrd;SM 

U9.sgrd", GRID="SOLUM_rk.sgrd", SHAPES="baranja.shp", FIELD=0, 
MODEL=1, NUGGET=0, SILL=200,
RANGE=500, INTERPOL=0))

see also http://geomorphometry.org/view_scripts.asp?id=24

But the last version of SAGA had some bug in this module, so I 
remember that I did not get correct
results (I did not test the most recent version of SAGA).

I personally think that we should simply think about implementing 
local RK (regressions/variograms
in moving window; maps of regression coefficients, variogram 
parameters, R-square etc.), and this
would then solve this problem (+give us much more insight into the data).


Tom Hengl
http://spatial-analyst.net


-----Original Message-----
From: gstat-info-bounces at geo.uu.nl [mailto:gstat-info-bounces at 
geo.uu.nl] On Behalf Of Edzer Pebesma
Sent: donderdag 21 augustus 2008 13:07
To: Vasat, Radim
Cc: Gstat-info at geo.uu.nl; sig-geo
Subject: Re: [Gstat-info] LOCAL universal kriging with GLOBAL reg. 
coeff.estimation

Good question, I'll include r-sig-geo as well.

actually the prediction equations you end up with are kind of funny, and
I've never seen them written out. Two different covariance matrices
being inversed.

Package gstat can do the two-step approach: global BLUE, then simple
kriging of residual, but that ignores the correlation of both terms,
when added.

I'm quite sure that if you do Gaussian simulation however with gstat,
the trend is fitted globally (and simulated from the corresponding GLS
mean & covariance), whereas the kriging is done locally, so that should
result in realisations that have the variance you want. You just have to
compute the average and variance, to get estimates of the kriging mean
and variance you're looking for. Also, it might be computationally
demanding, both in time and space (memory).
-- 
Edzer

Vasat, Radim wrote:

Dear all,

My question is how to performe LOCAL universal kriging with
estimation of regression coefficients from the whole area
(GLOBAL estimation) at once.

For large data (I have more than 8000 records) is computationaly too
demanding to handle with all the data for UK prediction. Hence,
I prefer to do it LOCALLY (namx=100), but the reg. coeff.
would be estimated only with limited data in this case (not all
data are included).

To estimate the reg. coeff. first (with BLUE=TRUE) and then perform
the simple kriging with beta equals to the reg. coeff. and nmax=100
might be the
solution. But the reg. coeff. estimation error is not included into the
final kriging variance in this case.

Anybody knows how to join these two requirements (LOCAL UK and GLOBAL
estinmation of reg. coeff.) into one prediction???

I will appreciate any coments!
Radim