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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Edzer Pebesma Sent: donderdag 21 augustus 2008 13:07 To: Vasat, Radim Cc: [EMAIL PROTECTED]; 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 > > _______________________________________________ > Gstat-info mailing list > [EMAIL PROTECTED] > http://mailman.geo.uu.nl/mailman/listinfo/gstat-info > -- Edzer Pebesma Institute for Geoinformatics (ifgi), University of Münster, Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/ _______________________________________________ Gstat-info mailing list [EMAIL PROTECTED] http://mailman.geo.uu.nl/mailman/listinfo/gstat-info _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo