[R-sig-Geo] R: GSTAT: Optimize power value for IDW (Paul Hiemstra approach)

[Alessandro]

Hi 

Normally I use the R+SAGA to calculate the IDW and create a raster, with
this follow code. I change the radius input with a loop formula to create
several raster and after check the best result (I am studying a oak forest
with low density) 

radii.list <- as.list(c(5, 10, 15, 20, 25, 30)) 
for(i in 1:length(radii.list)){ 
	rsaga.geoprocessor(lib="grid_gridding", module=0,
param=list(GRID=set.file.extension(paste("DEM_iw",radii.list[[i]],sep=""),".
sgrd"),
	SHAPES="ground.shp", FIELD=0, TARGET= 0, POWER=1.0,
RADIUS=radii.list[[i]], NPOINTS=20, 	USER_CELL_SIZE=dem.pixelsize,
USER_X_EXTENT_MIN=ground at bbox[1,1], USER_X_EXTENT_MAX=ground at bbox[1,2],
USER_Y_EXTENT_MIN=ground at bbox[2,1], USER_Y_EXTENT_MAX=ground at bbox[2,2])) 
} 


After I have 6 raster (DEM_idw_5.sgrd, DEM_idw_10.sgrd, DEM_idw_15.sgrd,
DEM_idw_20.sgrd, etc. etc.). I check hand-made the best IDW. I don't know is
It possible to improve this formula with RMSE in gstat and calculate the
best power.

Ale



-----Messaggio originale-----
Da: r-sig-geo-bounces at stat.math.ethz.ch
[mailto:r-sig-geo-bounces at stat.math.ethz.ch] Per conto di Paul Hiemstra
Inviato: martedì 2 dicembre 2008 12.07
A: Zev Ross
Cc: r-sig-geo at stat.math.ethz.ch
Oggetto: Re: [R-sig-Geo] GSTAT: Optimize power value for IDW

Zev Ross schreef:
> Hi All,
>
> ArcGIS has a nice little button that calculates the optimal power 
> value to use for inverse distance weighting based on cross-validation 
> and RMSPE. Just wondering if anyone had written something similar in R 
> -- I'm using GSTAT and I'd like to avoid back and forth with ArcGIS 
> (and obviously I'd like to avoid writing it myself as well!).
>
> Zev
>
Hi,

I don't have any code lying around, but a brute force optimization 
approach should be quite easy. Also because the range of idw-powers is 
relatively small. The speed would ofcourse depend on the size of the 
dataset. In code it would look something like (actually works :)):

library(gstat)
data(meuse)
coordinates(meuse) = ~x+y

# make function to do the cv
do_cv = function(idp) {
  meuse_idw = gstat(id = "zn", formula = zinc~1, data = meuse, set = 
list(idp = idp))
  out = gstat.cv(meuse_idw)
  return(sqrt(mean(out$residual^2))) # RMSE as optimization criterium
}

idw_pow = seq(0.1,5, by = 0.2) # the idwpower values that will be checked
cv_vals = sapply(idw_pow, do_cv) # calculate the rmse
# List of outcomes
print(data.frame(idp = idw_pow, cv_rmse = cv_vals))

After this you select the idw value with the smallest RMSE.

cheers and hth,
Paul

-- 
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:     +31302535773
Fax:    +31302531145
http://intamap.geo.uu.nl/~paul

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo