2

I want to perform IDW interpolation using R, and more specifically, using the idw command from the gstat package. Also, I want to use the optimal number of neighbors (point locations) and distance power, which will be determined by what combination of "idp" and "nmax" produces the lowest RMSE in leave-one-out-cross validation. For this reason, after each run of IDW with different "nmax" and "idp" parameters, I perform the leave-one-out cross validation. In order to do this, I have written a double for loop in R that tests each "idp" parameter for each number of "nmax" neighbors, but my script can't produce the result I want.

power = seq(from = 1, to = 1.6, by = 0.1)
neigh = seq(from = 4, to = 5, by = 1)
# January
for(i in power) {
  for(j in 1:length(neigh)
    jan_idw = idw(Jan ~ 1, st, grd, nmax = power, idp = neigh)
    list_jan_idw = c(list_jan_idw,jan_idw )
    krige.cv(Jan ~ 1, st, nfold = 52)
  }
} 

For each run, I want the results to be stored in a data frame or matrix.

  • doesn't krige.cv already give you a data.frame? If you want to save all the data.frames from your loops, use a list. Then you will have a list with length(power)*length(neigh) data.frames. Very simple. – aaryno Jun 17 '15 at 15:29
4

You could use a list to save all the krige.cv data.frames. I changed your example to the meuse data to make it reproducible. I'm assuming your data is similarly formatted.

require(gstat)
require(sp)
data(meuse)
data(meuse.grid)
coordinates(meuse) = ~x+y
coordinates(meuse.grid) = ~x+y 
gridded(meuse.grid) = TRUE

power = seq(from = 1, to = 1.6, by = 0.1)
neigh = seq(from = 4, to = 5, by = 1)
# January
results=list()
results.cv=list()
for (i in power) {
  for (j in neigh) {
    results[[paste0(i,"_",j)]] = idw (zinc ~ 1, meuse, meuse.grid, nmax = i, idp = j)
    results.cv[[paste0(i,"_",j)]] = krige.cv (zinc ~ 1, meuse, nfold = 52)
  }
} 

to get access to a specific list element, for example for power=1, neigh=4:

results[["1_4"]] # idw model

To see how the results are stored in the lists:

> str(results)
List of 14
 $ 1_4  :Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
  .. ..@ data       :'data.frame':  3103 obs. of  2 variables:
  .. .. ..$ var1.pred: num [1:3103] 634 713 654 604 857 ...
  .. .. ..$ var1.var : num [1:3103] NA NA NA NA NA NA NA NA NA NA ...
  .. ..@ coords.nrs : int [1:2] 1 2
  .. ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
  .. .. .. ..@ cellcentre.offset: Named num [1:2] 178460 329620
  .. .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
  .. .. .. ..@ cellsize         : Named num [1:2] 40 40
  .. .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
  .. .. .. ..@ cells.dim        : Named int [1:2] 78 104
  .. .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
  .. ..@ grid.index : int [1:3103] 69 146 147 148 223 224 225 226 300 301 ...
   ...

    >  str(results.cv)
List of 14
 $ 1_4  :Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  .. ..@ data       :'data.frame':  155 obs. of  6 variables:
  .. .. ..$ var1.pred: num [1:155] 796 738 606 481 369 ...
  .. .. ..$ var1.var : num [1:155] NA NA NA NA NA NA NA NA NA NA ...
  .. .. ..$ observed : num [1:155] 1022 1141 640 257 269 ...
  .. .. ..$ residual : num [1:155] 225.5 402.8 33.9 -223.6 -99.7 ...
  .. .. ..$ zscore   : num [1:155] NA NA NA NA NA NA NA NA NA NA ...
  .. .. ..$ fold     : int [1:155] 1 36 21 40 48 36 39 18 18 32 ...
  .. ..@ coords.nrs : num(0) 
  .. ..@ coords     : num [1:155, 1:2] 181072 181025 181165 181298 181307 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:155] "1" "2" "3" "4" ...
  .. .. .. ..$ : chr [1:2] "x" "y"
  .. ..@ bbox       : num [1:2, 1:2] 178605 329714 181390 333611
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:2] "x" "y"
  .. .. .. ..$ : chr [1:2] "min" "max"
  .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. .. ..@ projargs: chr NA

You can be more creative about how you name your list indices. You can also make a list of lists (e.g., results[[i]][[j]], in which case you would have to initialize each results list as a new list outside the inner loop, as in results[[i]]=list())

  • Thank you very much @aaryno. I run it and it came up with the following warnings: "In if (d$nmax == Inf) nmax = as.integer(-1) else nmax = as.integer(d$nmax) : the condition has length > 1 and only the first element will be used" Also, when I do result[[power[i]_neigh[j]] it always returns NULL. – Maria Karypidou Jun 17 '15 at 16:04
  • for loops were not formulated correctly. I changed the code to use the example meuse so I could run it through. – aaryno Jun 17 '15 at 16:48
  • Thank you very much @aaryno. This is exactly what I wanted the script to do! – Maria Karypidou Jun 17 '15 at 17:59
  • 1
    I fail to see how the power and neighbors are passed to the krige.cv call. If I don't specify a model in krige.cv it defaults to IDW? – Daniel Mar 14 '17 at 13:32
  • 2
    Just found a way to set idp in krige.cv call (gis.stackexchange.com/questions/190558/…) – Daniel Mar 14 '17 at 14:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.