Search code examples
rpolygonrastersamplingspatstat

R sample a raster with square polygons

I would like to sample a big raster by creating In small raster 100x100 cells. I don't know how to do that so any ideas are welcome

My actual lead : 

library(raster)
library(spatstat)
library(polyCub)

r <- raster(ncol=1000,nrow=1000) # create empty raster
r[] <- 1:(1000*1000)             # Raster for testing
e <- extent(r)                   # get extend
# coerce to a SpatialPolygons object
p <- as(e, 'SpatialPolygons')  


nc <- as.owin.SpatialPolygons(p) #polyCub
pts <- rpoint(50, win = nc)
plot(pts)

Now I need to generate 100x100 cell square around my 50 points and I would like to crop r using those square and stack each small raster individually ...

Solution


  • The answer by @adrian-baddeley basically has the ingredients to do what you want. If you simply want a list of small im objects that contain the 100x100 box you simply subset im objects by owin objects to extract the relevant region. Here is an example (with fewer points to avoid overplotting)

    library(raster)
library(spatstat)
library(maptools)

r <- raster(ncol=1000,nrow=1000) # create empty raster
r[] <- 1:(1000*1000)             # Raster for testing
e <- extent(r)                   # get extend
# coerce to a SpatialPolygons object
p <- as(e, 'SpatialPolygons')  

nc <- as.owin.SpatialPolygons(p)
set.seed(42)
pts <- rpoint(7, win = nc)

rim <- as.im.RasterLayer(r)
Box <- owin(c(-50,50) * rim$xstep, c(-50,50) * rim$ystep)

    The following is a list of im objects of size 100x100

    imlist <- solapply(seq_len(npoints(pts)),
                   function(i) rim[shift(Box, pts[i])])

    Here is a plot of the im objects in the region and the points on top

    plot(pts)
for(i in imlist) plot(i, add = TRUE)
plot(pts, pch = 19, add = TRUE)

    You can convert to a list of raster layers with

    rasterList <- lapply(imlist, as, Class = "RasterLayer")

    PS: The following is a list of im objects of the original size with NA outside the 100x100 box if you need that format instead

    imlist <- solapply(seq_len(npoints(pts)),
                   function(i) rim[shift(Box, pts[i]), drop = FALSE])