[R-sig-Geo] bootstrapping linear model and residuals kriging by RSAGA

[Vasya Pupkin]

Dear ladies and gentlemen, I don’t really know it this is a suitable platform for asking such questions, but stackoverflow could not help me, so I am trying my luck here:).
What I am doing: I am trying to make a soil erodibility map. I have soil erodibility measurements at points and I do a linear regression analysis to predict soil erodibility from predictors (morphological and remotely sensed). Then I apply the regression equation to the predictor rasters and get a raster of soil erodibility. Then I calculate residuals and extra/inter-polate with ordinary kriging (global) using my own formula discribing the variogramm. Then I summ up the map from regression equation and the interpolated rsiduals to create the final map, thus mimiking regression kriging as you have suggested. The regression was done in R.
Now I want to validate the final map with bootstrap as it is implimented in R "boot" package. In the madual it says you can wrap any complicated model in the "boot" sequence. So I wanted to write this whole regression+residual_kriging model in the boostrap in R. After googling I fould this RSAGA package that provides a direct access to SAGA's functionality. And there is this "rsaga.geoprocessor" module which works like CMD for SAGA, and will not return the results back to R environment for further calculations but will save them to a file, like SAGA's CMD normally does. And there is this "pick.from.points" module which has sort of deeper integration with SAGA's modules (not all though) and is what I think I need but I can't figure out how I can use my own variogram formula and how I can return the results back to R environment for further calculations. So in other words I need to do a regression and then predict the values and extrapolate residuals for excluded points, so the RSAGA module should return the values of extrapolated residuals back to “boot” calculation.
What I want to do is something like this (R script):
boot(mydata,
       lm( resp~predict, ...)+
       pick.from.points( lm(resp~predict)$residuals, method = "krige", model = "a+b^c",...))
Do you have any idea how it could be done? If not, could you please advise on any other ways or anybody who could know? The approach – regression + residuals kriging is critical, but not the software, so if there are any other packages for doing this, please advice. 

Best wishes,
Maksim.

	[[alternative HTML version deleted]]