select raster within a specified distance from polygon boundary

Question

I want to select raster cells that are within a certain distance (for e.g. 1 km or 5 km) from the boundary of a polygon. I ultimately want to take an average of only those raster cells that are within the specified distance from the boundary of shapefile inwards.

The way I thought I would approach is to create a negative buffer inwards, and subtract the original polygon and the buffer. Then mask and crop the raster using the new polygon and take the average.

Here's sample data demonstrating what I want to do.

library(raster)

# raster
r <- raster(xmn=1035792, xmx= 1116792, ymn=825303.6, ymx=937803.6, resolution = 12.5,crs = "+init=epsg:3174")
r <- setValues(r, 0)

# polygon
x <- c(1199999, 1080000, 1093067, 1090190, 1087977, 1070419, 1180419)
y <- c(957803.6,937803.6, 894366.9, 872153.9, 853703.0, 825353.6, 805353.6)
poly.lake <- SpatialPolygons(list(Polygons(list(Polygon(data.frame(x,y))), ID = 1)))

r <- mask(r, poly.lake)
r <- crop(r, poly.lake)

plot(poly.lake)
plot(r, add = T)

Instead of taking average of the resulting raster r, I only want to average raster cells which are within a certain specified distance from the boundary.

So, create a sampling 'donut'[s]? p_buffer2 <- buffer(p_buffer, -1000), plot(r), plot(p_buffer, col = 'red', add = TRUE), plot(p_buffer2, col = 'blue', add = TRUE), of the 'red'? That might be better to do on a p[1]:p[n] basis. — Chris, Oct 09 '22 at 04:18
[sf polygon with hole](https://stackoverflow.com/questions/64055453/make-sf-object-out-of-polygons-with-holes-and-set-crs), blue is the hole. wrapped in `vect`? — Chris, Oct 09 '22 at 04:28
sorry the command `p_buffer2 <- buffer(p_buffer, -1000)` gives me `Error: [buffer] invalid index` — 89_Simple, Oct 09 '22 at 09:03
I found similar question here which is what I am also aiming to do https://gis.stackexchange.com/questions/382173/cropping-a-raster-with-a-buffered-area-around-a-polygon-line-in-r — 89_Simple, Oct 09 '22 at 22:36
Perhaps mixed with [this-gis.stackexchange.com/questions/312289](https://gis.stackexchange.com/questions/312289/r-create-multiple-linestrings-from-multiple-coordinates) — Chris, Oct 09 '22 at 22:49
Are you interested in extracting raster values which are neighbouring any of the boundaries of your polygon of just a specific one? Your illustration focusses only on the western one, so how do we now (since it seems like polygon to polyline conversion is necessary here) which side is relevant to you? — dimfalk, Oct 10 '22 at 14:34
Actually I have a country shapefile and my raster (which represents coastal flooding) covers part of the land. I am interested in looking at average flooding value within a certain distance from the shoreline. — 89_Simple, Oct 12 '22 at 11:32

Robert Hijmans · Accepted Answer · 2022-10-11T19:28:23.073

The example data but using "terra"

library(terra)
r <- rast(xmin=1035792, xmax= 1116792, ymin=825303.6, ymax=937803.6, resolution = 125, crs = "epsg:3174")
values(r) <- 1:ncell(r)

# polygon
x <- c(1199999, 1080000, 1093067, 1090190, 1087977, 1070419, 1180419)
y <- c(957803.6,937803.6, 894366.9, 872153.9, 853703.0, 825353.6, 805353.6)
p <- vect(cbind(x, y), "polygons", crs = "epsg:3174")

r <- mask(r, p)
r <- crop(r, p)

You can now take the internal buffer of p

b <- buffer(p, -10000)
x <- mask(r, b, inverse=TRUE)
global(x, mean,na.rm=T)
#          mean
#lyr.1 296549.9

Or you can take both sides like this

bb <- buffer(as.lines(p), 10000)
y <- mask(r, bb)
global(y, mean,na.rm=T)
#          mean
#lyr.1 296751.3

So there is a slight difference between these two approaches; I think because the first uses inverse=TRUE; I would go with the second approach.

Your drawing (and Chris' answer) suggests that you only want the distance to the western border. In that case, you can first find the start and end nodes you need (from 2 to 6)

plot(p)
points(p)
text(as.points(p), pos=2)

Select the segments in between these nodes and create a line type SpatVector.

g <- geom(p)
k <- vect(g[2:6,], "lines", crs=crs(p))
lines(k, col="red", lwd=2)

And now do as above.

bk <- buffer(k, 10000) 
z <- mask(r, bk)
global(z, mean,na.rm=T)
#        mean
#lyr.1 297747

If you wanted to get the part of buffer bk that is inside the original polygon p you can do

bki <- intersect(bk, p)

To complete the plot

polys(bk, lty=3, border=NA, col=adjustcolor("light blue", alpha.f = 0.4))
lines(bki, lty=3)

To be fair, I thought there was some additional 'desired' sampling area expressed through the use of the hashed line, restricting sampling to hash to interior line. — Chris, Oct 11 '22 at 19:50

Chris · Answer 2 · 2022-10-11T19:35:32.170

Finding which segments of a polygon to buffer was what puzzled me, and this seems a decent approach cast_poly_to_subsegments. Taking your poly.lake as poly_sf:

geom <- lapply(
  1:(length(st_coordinates(poly_sf)[, 1]) - 1),
  function(i) {
    rbind(
      as.numeric(st_coordinates(poly_sf)[i, 1:2]),
      as.numeric(st_coordinates(poly_sf)[i + 1, 1:2])
    )
}
+ ) |>
  st_multilinestring() |>
  st_sfc(crs=st_crs(rt)) |>
  st_cast('LINESTRING')

gives us

which is a little surprising, the 'green and red', that I assumed would be 'green'. It is wound clockwise so the desired segments to buffer are 4 & 5.

lns_buf4 <- st_buffer(st_geometry(geom)[4], 1000, singleSide = TRUE)
lns_buf5 <- st_buffer(st_geometry(geom)[5], 1000, singleSide= TRUE)
lns_buf5_neg <- st_buffer(st_geometry(geom)[5], -1000, singleSide= TRUE)

plot(st_geometry(geom), col = c('red', 'yellow', 'blue', 'green'))
plot(lns_buf4, col = 'black', add = TRUE)
plot(lns_buf5, col = 'green', add = TRUE)
plot(lns_buf5_neg, col = 'blue', add = TRUE)

Whether +/-1000 is sufficient is a further intersection test between the buffer poly(s) and the other boundary. If the desired sampling area is not rectangular, steps can be taken to construct a sampling polygon from the buffer and intersection.

#library(lwgeom)
# on poly_sf 
new_line <- draw(x = 'line', col ='blue', lwd = 2, n = 10)
lns_buf5_10k_neg <- st_buffer(st_geometry(geom)[5], -10000, singleSide= TRUE)
new_line_sf <- st_as_sf(new_line, crs = st_crs(lns_buf5_10k_neg))
buf5_nline_split <- lwgeom::st_split(lns_buf5_10k_neg, new_line_sf$geometry)
irreg_smp_area <- st_collection_extract(buf5_nline_split)[1]

Though I'm happy to see it all done in terra.

select raster within a specified distance from polygon boundary

2 Answers2