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

[Edzer Pebesma]

Hi Zev,

There's no need to brute force, as optimize is there to help you -- my 
guess is that the function is convex. The following takes a while:

 > f = function(idp, formula, data,...) 
sum(krige.cv(formula,data,set=list(debug=0,idp=idp),...)$residual**2)
 > optimize(f, interval=c(0.01,4), formula=log(zinc)~1, data=meuse)
$minimum
[1] 3.532051

$objective
[1] 32.09126

but that is probably due to the (my?) hopelessly inefficient (though 
flexible) setup of krige.cv. Speeding up can be done by allowing a 
larger tolerance, or passing e.g. nfold=5 to optimize(). This nfold will 
also result in somewhat random output, as it randomly folds the data in 
5 partitions.
--
Edzer


Zev Ross wrote:
> 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,
>> 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
>>
>>   
>
> -- 
> Zev Ross
> ZevRoss Spatial Analysis
> 120 N Aurora, Suite 3A
> Ithaca, NY 14850
> 607-277-0004 (phone)
> 866-877-3690 (fax, toll-free)
> zev at zevross.com
> ------------------------------------------------------------------------
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>   

-- 
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/
http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de