2

I have numerous large stacked raster file with multiple layers and I would like to select 250 random points and extract the data values from all layers for these points.

I had been using the raster package sampleRandom(raster,n, na.rm=TRUE, xy=TRUE) to select random points before extracting the data from my raster object but this works without replacement. I, however, would like to sample random points from the raster files with replacement as can be done with the r base package sample(x, size, replace = TRUE, prob = NULL).

I had considered just running a loop using sampleRandom so that each point has a chance of being sampled for each of the 250 runs;

samp_all <- matrix(, nrow = 250, ncol = 10)
for (i in 1:250)
{samp <- sampleRandom(layer1, 1, xy = TRUE, na.rm = TRUE)[, -3]
samp_all[i,] <- samp}
samp2 <- extract(stack_all, samp_all)

But this seems to run pretty slowly and I have quite a lot of files I would like to run this on.

Is there a simpler way to randomly sample points from raster layers with replacement that would allow me to extract data from my raster?

  • Would you like to get the same data points for all rasters? E.g. if a data-point at (row == 1; col == 5) was sampled, would you like to get its values from all raster layers? or would you like a different sample for each layer? – dof1985 Aug 16 '16 at 9:44
  • I'd like to get the same data points from all raster layers, just with the added possibility that some data points might get selected more than once. – J. Cee Aug 16 '16 at 9:49
  • You could just modify the original raster package function and source it back to R. You can always download source code from the packages CRAN website. There is only one line of code in the sampleRandom function that would need to be changed to add replacement. 'cells <- sample(ncell(r), size = size, replace = TRUE)' Once you have done this just use the source() function or copy-and-paste the function into the R console. – Jeffrey Evans Jul 6 '18 at 16:44
0

This function also does the trick, and on top of that samples only from cells that are not NA. It mimics dismo::randomPoints() but runs much faster. It returns a matrix, just like the dismo function.

# r is a raster, n is the number of points requested
fastRandomPoints <- function(r, n) {
  if(raster::nlayers(r) > 1) r <- r[[1]]
  v <- raster::getValues(r)
  v.notNA <- which(!is.na(v))
  x <- sample(v.notNA, n)
  pts <- raster::xyFromCell(r, x)
  return(pts)
}
  • The reason that it is faster is that you have removed the memory-safe aspects of the function. There is no real problem with that excepting when somebody wants to sample a raster with billions of cells. – Jeffrey Evans Jul 6 '18 at 16:39
  • Thanks for your reply. Can you explain what you mean by memory-safe, and why the function above would be a problem for a very large raster? – ndimhypervol Jul 7 '18 at 20:54
  • On large rasters the raster package uses gdal to read and write blocks for processing. In this way the raster is never fully read into memory. Certain function calls (eg., getValues) force the raster to be read. If you look at @dof1985's answer, it is using the row/column index so, never reads a vector of raster values. However, it also does not account for NA's. – Jeffrey Evans Jul 8 '18 at 3:21
  • Thanks for the clarification. So if that is true, can you explain why this function runs much faster than the regular randomPoints() even though getValues() reads the raster into memory? Can't one just rm() the raster out of memory when the function is done? – ndimhypervol Jul 9 '18 at 16:05
  • When you can get away with it, it is faster, processing wise, to operate on objects in memory. If raster is swapping blocks of data it is bound to slow things down. However, this is also a memory asymptote where processing really slows down on very large data. It is not a matter of removing the object, that is just good practice, but rather having a means of processing say, 10 billion observations. This is where the raster package advantage of reading/writing blocks of digestible data becomes very desirable. – Jeffrey Evans Jul 9 '18 at 16:53
1

Try the following code. The idea is to create a set of [row, col] indicators based on the RasterStack dimensions (rows and cols). Than you can easily used these indicators on the stack to subset all values.

The function below sampleStack gets a RasterStack and n number of values to sample, and gives a data frame with [row, col] positions and values extracted by layer.

library(raster)
# Generate raster layers

r <- raster(matrix(rnorm(100, 0, 1), nrow = 10))
for (i in 1:5) {r <- stack(r, raster(matrix(rnorm(100, 0, 1), nrow = 10)))} # run 5 times

sampleStack <- function(r, n) {
  rowSample <- sample(1:r@nrows, size = n, replace = TRUE)
  colSample <- sample(1:r@ncols, size = n, replace = TRUE)
  pairs <- data.frame("rowInd" = rowSample, "colInd" = colSample)
  out <- as.data.frame(cbind(pairs, as.data.frame(t(apply(pairs, MARGIN = 1, function(x) {return(r[x[1],x[2]])})))))
  colnames(out)[3:ncol(out)] <- names(r)
  return(out)
}

# Example Run
sampleStack(r = r, n = 14)

   rowInd colInd   layer.1.1  layer.2.1   layer.1.2  layer.2.2    layer.1    layer.2
1       5      4 -0.09678111  1.6844843  0.82574090  0.5328165  0.66721846  0.1936958
2       4      8  0.64982724  0.5467126  1.59975344  0.2757094  0.94797866 -0.1798319
3       3      4 -1.35927393 -1.2774878  0.77616160  0.1429519  1.10396643  1.1444793
4       7      9  1.45380719 -0.6128730 -0.53011041 -0.3138787 -0.86586255  0.9056694
5       6      2 -0.49808353 -0.1272448 -1.96004940 -0.5663870 -0.07217682 -1.7568981
6       4     10  0.01607546 -0.5113896  1.19713933 -1.6322803 -1.04051134  0.7135125
7       2      4  1.10798593 -0.2610036  0.56009222  2.4618433 -0.44356484  1.0332427
8       3      1  1.70168812  1.7643488 -0.09976064 -0.2386893 -1.04266622  0.2019014
9      10      3 -1.82791351  1.5666126 -1.79275437  0.2946007  0.96467732 -0.6951626
10      6      4  0.54795051  0.1378088 -1.53793046 -0.5989934 -1.64424273 -0.3463153
11      1      2 -0.86287672 -0.2408750 -0.81438516 -2.0200205  1.16523355  0.4052408
12      2      2 -0.47916254 -0.6778470  0.79086436 -0.5692255  0.96205715 -0.5146865
13      5      4 -0.09678111  1.6844843  0.82574090  0.5328165  0.66721846  0.1936958
14      6      1 -1.04973148 -0.3973457 -0.24445969  0.4061588 -1.50143806  0.3896232
  • 1
    Thanks that works great pulling everything out together. I did also manage to overcome the non replacement issue of sampleRandom by sampling based on cell and then extracting xy values. cell<-ncell(layer1) cells<-sample.int(ce‌​ll, 250, replace = TRUE) samp_all<-xyFromCell‌​(layer1, cells, spatial=FALSE) samp2 <- extract(stack_all, samp_all) – J. Cee Aug 16 '16 at 14:25
  • I really like this approach however, it does not account for NA's. Often rasters are irregular, with a considerable amount of nodata along the edges of the raster. If this is not accounted for you can plausibly end up with a large portion of your random same being NA. This also returns values but not the coordinates. Often you want a spatial object (points) returned for additional analysis, which is an option in the raster::sampleRandom function with sp=TRUE. – Jeffrey Evans Jul 6 '18 at 16:49
  • 1
    You are right, this does not account for NAs. @ndimhypervol have provided a good approach to handle NAs in is answer. Regarding coordiantes, I dont think it was specified by the OP. Anyhow the function could be re-arranged to include both. – dof1985 Jul 7 '18 at 20:44

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.