Hi All, Thanks to Paul and Alessandro for their suggestions. Paul's code (brute force) worked well for me and the results match up well with ArcGIS. I'm not using a large dataset so the speed isn't an issue but with a larger dataset it would be. In ArcGIS the optimization is instantaneous so clearly the software is doing something different than testing out all possible values. I used Paul's code with no substantive changes then it's straightforward to use: plot(idw_pow, cv_vals) idw_pow[which.min(cv_vals)] And pull out the min. Thanks for the advice! Zev Alessandro wrote: 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, [EMAIL PROTECTED],1], [EMAIL PROTECTED],2], [EMAIL PROTECTED],1], [EMAIL PROTECTED],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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]] Per conto di Paul Hiemstra Inviato: martedì 2 dicembre 2008 12.07 A: Zev Ross Cc: r-sig-geo@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!). ZevHi, 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 = "" 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 -- Zev Ross ZevRoss Spatial Analysis 120 N Aurora, Suite 3A Ithaca, NY 14850 607-277-0004 (phone) 866-877-3690 (fax, toll-free) [EMAIL PROTECTED] |
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo