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