3

I have a raster dataset that is classified into vegetation classes. I would like to randomly sample from the dataset based on the vegetation classes, but not sure how to do it.

Here is some example code taken from the deratify() function in the package raster.

library(raster)
r <- raster(nrow=10, ncol=10)
r[] = 1
r[51:100] = 2
r[3:6, 1:5] = 3
r <- ratify(r)
rat <- levels(r)[[1]]
rat$landcover <- c('Pine', 'Oak', 'Meadow')
rat$code <- c(12,25,30)
levels(r) <- rat
x <- deratify(r, 'landcover')

What I would like to do is to sample N number of cells (or points) from each class in x (Pine, Oak, Meadow).

Any hints on how to do this without subsetting the raster?

1

Here is a function that does just that. Some comments will follow the code.

sampleClasses <- function(r = raster, n = 8) {
 # function gets a raster object and a sample size. It samples n cells from each class
 # in the raster and return a vector with the cells indices
 vals <- unique(getValues(r)) # Get all classes
 for (val in vals) {
  cellVal <- which(t(as.matrix(r) )== val) # get All cells for each class
  if(!exists("samples")) {
    samples <- sample(cellVal, n) # sample class locations
  } else {
    samples <- c(samples, sample(cellVal, n))
  }
 }

 return(samples)
}

So basically you can run this function and call it with your raster object, and the sample size you would like to have. Note that it always uses the same n, so you can't get different samples size for different classes. Also note that the function get all classes from the raster, thus you can't choose to sample only some of them. However, both behaviors can be easily modified and adjusted to a fit other (or more general) operations.

Following is a running example on your data.

> r[sampleClasses(r = r, n = 5)]

>  [1] 1 1 1 1 1 3 3 3 3 3 2 2 2 2 2 # these are the values.


> sampleClasses(r = r, n = 5) # Function returns the cells indices

>  [1] 46 16 38 37  4 41 45 34 23 53 79 86 63 90 77

>  xyFromCell(r, sampleClasses(r = r, n = 5)) # You can also get xy coords of the cells.

x y

[1,] 126 63

[2,] -162 63

[3,] 54 9

[4,] 18 63

[5,] -126 81

[6,] -126 45

[7,] -162 -9

[8,] -54 -9

[9,] -90 27

[10,] -162 45

[11,] -162 -27

[12,] 54 -81

[13,] 54 -45

[14,] 18 -81

[15,] -162 -63

  • 1
    That is exactly what I was looking for. – user44796 Jul 27 '16 at 15:30
  • 2
    The issue with this function is that it not memory safe. This is an important consideration, not only with providing an answer but, also when writing code in general. It would be quite easy to max out available RAM on a somewhat moderate sized raster. You always have to think beyond "what has worked for you" to how may code be utilized once it is out in the world. A not so efficient but safer approach would be to recalssify each class to a single raster [1 else NA] and then use raster::sampleRandom to draw the random sample for that class. You could also coerce to an sp object and use tapply. – Jeffrey Evans Jun 18 '19 at 18:40
  • @JeffreyEvans - thank you. This is a valuable comment and lesson. – dof1985 Jun 24 '19 at 8:32
3

For a non memory-safe option you can efficiently use an apply type function on an sp SpatialPointsDataFrame object. This has the advantage of not having to use which, which replicates a vector the size of your raster (which could be huge), is much faster than a for loop and directly results in an sp point object.

Add libraries and create data

library(sp)
library(raster)

r <- raster(nrow=100, ncol=100)
  r[] <- round(runif(ncell(r), 1,5),0)

Coerce to SpatialPointsDataFrame object.

r.sp <- as(r, "SpatialPointsDataFrame")
  names(r.sp) <- "class"

Now, use tapply to pull random samples from each class. This application of the tapply function results in a list object containing the random sample indices by each class value. You can then subset or create a new object of these sample indices and the result is a SpatialPointsDataFrame representing the random samples for each class. I am doing this in one fell swoop.

p=5 # size of random sample for each class
r.sp <- r.sp[unlist(tapply(1:nrow(r.sp), list(r.sp$class),  
             FUN=function(x) x[sample(1:length(x), p)])),]

This is what the results of tapply look like outside the subset.

tapply(1:nrow(r.sp), list(r.sp$class),  
      function(x) x[sample(1:length(x), p)])

Plot resulting random sample(s) on raster.

plot(r)
  points(class.samp, pch=19)

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.