[GIS] Segmenting rasters into sub-regions and sample using R without using extract()

rrasterspatial-analysis

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...

Best Answer

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 )
Related Question