R sample a raster with square polygons

Question

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

Your objective is unclear to me: Are the given inputs a big raster in the format of the `raster` package and some (50) given points in the `ppp` format of `spatstat` and then you want to extract 100x100 raster cells around each point? What do you mean by stacking each small raster individually? — Ege Rubak, Nov 16 '17 at 21:42
Also there must be edge cases to think of: E.g. what should happen when a point is close to the raster border and the raster doesn't contain a 100x100 grid around it? — Ege Rubak, Nov 16 '17 at 23:07
Sorry for the delay, _i)_ You had understood exactly what I'm looking for! extract a 100x100 raster cell around each point. For me the best output will be a raster stack. _ii)_ I don't know (yet) own to manage point near the raster border. — delaye, Nov 27 '17 at 11:13

score 2 · Answer 1 · answered Nov 17 '17 at 05:03

If you want to use spatstat then you need to convert the raster object r into an object of class im supported by spatstat. You can do this conversion in the maptools package. Call this image object rim. Then you can do as follows

Box <- owin(c(-50,50) * rim$xstep, c(-50,50) * rim$ystep)
BoxesUnion <- MinkowskiSum(pts, Box)
W <- intersect.owin(as.mask(rim), BoxesUnion)

This would give you the subset of the raster that is covered by the squares. If you want to keep the squares separate, do something like

M <- as.mask(rim)
BoxList <- solapply(seq_len(npoints(pts)), 
                      function(i) intersect.owin(M, shift(Box, pts[i])))

Then BoxList is a list of the individual sub-rasters.

If I understand your proposition, `BoxList` is a List of windows ! It sounds great ! Now I don't understand how to use those windows to mask the initial raster or the `im`... I'm not comfortable with `owin` object. — delaye, Nov 27 '17 at 12:27

score 2 · Accepted Answer · answered Nov 28 '17 at 09:55

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])

R sample a raster with square polygons

2 Answers2