2

Is there a more efficient way to sample from a sub-region (defined by a shapefile read into a spatial object) of multiple rasters without having to use extract() for each region on every raster?

Data: I have a large study area represented by approximately 50 raster layers, one per variable (e.g. elevation, slope, vegetation, etc...) my study area is divided into 10 sub-region

Methods: I need to draw random samples (using sampleRandom() or sample() depending on object) from each region of each raster; therefore approximately 500 samples.

Issue: I can do this relatively easy by using the extract() function for small datasets, but it will take forever with these data sets. I am hoping there is a way to sample from within the area of the sub-regions without having to use extract()

Below is code (some thankfully borrowed from J. Evans) that shows the basics of what I am doing.
This works fine for small raster data sets, but trying to use extract() on large raster data sets fails.

require(raster)

r1 <- raster(ncol=10, nrow=10, xmn=0, ymn =0, xmx=100, ymx=100)
r1[] <- runif(ncell(r1),0,100)
r2 <- raster(ncol=10, nrow=10, xmn=0, ymn =0, xmx=100, ymx=100)
r2[] <- runif(ncell(r1),100,200)
r3 <- raster(ncol=10, nrow=10, xmn=0, ymn =0, xmx=100, ymx=100)
r3[] <- runif(ncell(r1),200,300)

region1 <- rbind(c(0,0), c(50,0), c(50,50), c(0,50), c(0,0))
region2 <- rbind(c(50,0), c(100,0), c(100,50), c(50,50), c(50,0))
region3 <- rbind(c(0,50), c(50,50), c(50,100), c(0,100), c(0,50))
region4 <- rbind(c(50,50), c(100,50), c(100,100), c(50,100), c(50,50))

polys <- SpatialPolygons(list(Polygons(list(Polygon(region1)), "region1"),Polygons(list(Polygon(region2)), "region2"), Polygons(list(Polygon(region3)), "region3"), Polygons(list(Polygon(region4)), "region4")))

plot(r1)
plot(polys, add=TRUE) 

r_stack <- stack(r1, r2, r3)
names(r_stack) <- c("r1", "r2", "r3")

for(i in 1:length(polys)){
    region <- polys[i]
    for(k in 1:length(names(r_stack))){
        raster_region <- extract(r_stack[[k]], region)
        sample_values <- sample(unlist(raster_region), 10)
        print(sample_values)
        ### add samples to some matrix...
    }
}
### Do analysis on values...
6

You could crop and mask the raster to each region and then use "sampleRandom" to create a random sample of the specific region(s). This should speed things up considerably.

        require(raster)

        r1 <- raster(ncol=10, nrow=10, xmn=0, ymn =0, xmx=100, ymx=100)
          r1[] <- runif(ncell(r1),0,100)
        r2 <- raster(ncol=10, nrow=10, xmn=0, ymn =0, xmx=100, ymx=100)
          r2[] <- runif(ncell(r1),100,200)
        r3 <- raster(ncol=10, nrow=10, xmn=0, ymn =0, xmx=100, ymx=100)
          r3[] <- runif(ncell(r1),200,300)

        r <- stack(r1, r2, r3)
          names(r) <- c("r1", "r2", "r3")  

        region1 <- rbind(c(0,0), c(50,0), c(50,50), c(0,50), c(0,0))
          region2 <- rbind(c(50,0), c(100,0), c(100,50), c(50,50), c(50,0))
            region3 <- rbind(c(0,50), c(50,50), c(50,100), c(0,100), c(0,50))
              region4 <- rbind(c(50,50), c(100,50), c(100,100), c(50,100), c(50,50))
        polys <- SpatialPolygons(list(Polygons(list(Polygon(region1)), "region1"),
                                 Polygons(list(Polygon(region2)), "region2"), 
                                 Polygons(list(Polygon(region3)), "region3"), 
                                 Polygons(list(Polygon(region4)), "region4")))

 rs <- as.data.frame(array(0, dim=c( 0, length(names(r))+1 )))  
      names(rs) <- c("REGION", names(r))                 
  for(i in 1:length(polys)){
      region <- polys[i]
        cr <- crop(r, extent(region), snap="out")                   
          m <- rasterize(region, cr)
           mr <- mask(x=cr, mask=m)
             rs.mat <- sampleRandom(mr, 10, na.rm=TRUE)
               rs.mat <- cbind(REGION=i, rs.mat)
        rs <- rbind(rs, as.data.frame(rs.mat) ) 
      }
    print( rs )
2
  • Thanks Jeffrey. I will try this with big rasters tomorrow, but it works great running the example. One issue I have had with using the raster stack is getting an IO error when cycling through a stack for further raster analysis. Calling the same raster by file path works fine, but called through the stack gives the generic IO error. Any idea? It is not an issue with extent and opening it through GDAL worked fine. – Mr.ecos Dec 6 '13 at 2:45
  • Implemented this logic into my functions; works great! I took the rasterize step out of the loop and made a separate function that rasterizes all of the regions before running this loop. So as not to re-run rasterize (the slowest of the functions) for each varaible being masked. Now only 9 hours of processing left... – Mr.ecos Dec 6 '13 at 21:46

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.