Search code examples
rrasterr-rasterrgdal

raster in R: Subset a raster image given a pixel size

I've like to subset a raster image given the coordinates center point and the desired x and y size in pixels. In my example:

library(raster)  
set.seed(5)

##Create a raster
r <- raster(ncol=10, nrow=10, crs="+proj=utm +zone=1 +datum=WGS84", xmn=0, xmx=50, ymn=0, ymx=50)
s1 <- stack(lapply(1:4, function(i) setValues(r, runif(ncell(r)))))

# Create 10 random points in the raster
pts <- data.frame(pts=sampleRandom(s1, 10, xy=TRUE)[,1:2], status=rep(c("A","B"),5))

Now, I've like to create the 10 new images that will be a subset of size (X desired width of subset image = 5 pixels and Y desired width of subset image = 5 pixels) with center point (pts$x, pts$y), with the same pixel size. Finally, a save each image as GeoTIFF:

for(i in 1:length(pts)){
  writeRaster(r[[i]],filename=paste(pts,"_",i,sep=""),
              format="GTiff",datatype="FLT4S",overwrite=TRUE)
}

This is possible?

Solution

  • Hey I was able to do this by creating square buffer around your points, big enough to capture a 5 by 5 raster. I then used those buffers to clip your raster. Here is my code with comments:

    library(raster)  
library(rgeos)
set.seed(5)

##Create a raster
r <- raster(ncol=10, nrow=10, crs="+proj=utm +zone=1 +datum=WGS84", xmn=0, xmx=50, ymn=0, ymx=50)
s1 <- stack(lapply(1:4, function(i) setValues(r, runif(ncell(r)))))

# Create 10 random points in the raster
pts <- data.frame(pts=sampleRandom(s1, 10, xy=TRUE)[,1:2], status=rep(c("A","B"),5))
pts_points = SpatialPointsDataFrame(pts[,1:2], data.frame(pts[,3]))

## Creating square buffers around your points
pts_buffer = gBuffer(pts_points, byid = TRUE, capStyle = "SQUARE", width = 12.5)


## looping through all your points
for (i in 1:nrow(pts_buffer)){
  square = pts_buffer[i,]
  val = square@data$pts...3. ## Gets the value of the point ("A" or "B")

  ## Contains a stack of raster clipped for the current point for the raster stack
  clipped = mask(s1, square) 

  ## Export the clipped raster here
}

    A width of 12.5 was used for the buffers to make sure a 5 by 5 raster was created. To change the subset size of the rasters, just change the width of the buffer!

    This is how the buffer points look like on top of the rasters: enter image description here

    In the code, We loop over each square and clip the raster to the area, here is how the raster stack looks like after being clipped by one of the areas:

    enter image description here

    Let me know if this helps, or if anything needs to be clarified!