Search code examples
rrasterr-sfterra

select raster within a specified distance from polygon boundary


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)

enter image description here

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.

enter image description here


Solution

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

    enter image description here