Search code examples
rbufferrasterspatialr-stars

extract mean value of raster with buffer condition on second layer/attribute

enter image description here enter image description here

I have this stars object (could be also formatted to raster):

stars object with 2 dimensions and 2 attributes
attribute(s):
   LST_mean        elevation    
 Min.   :14.98   Min.   :296.0  
 1st Qu.:16.89   1st Qu.:346.9  
 Median :17.64   Median :389.3  
 Mean   :17.52   Mean   :389.2  
 3rd Qu.:18.18   3rd Qu.:428.3  
 Max.   :20.11   Max.   :521.6  
dimension(s):
  from to  offset    delta                       refsys point values    
x   71 83 4387654  860.241 DHDN / 3-degree Gauss-Kru... FALSE   NULL [x]
y   33 41 5598885 -860.241 DHDN / 3-degree Gauss-Kru... FALSE   NULL [y]

Which has 2 attributes (layers in the case of raster): temperature and elevation. Using temperature, I would like to select the pixels that fall within a buffer and return the mean, only for the pixels whose difference in elevation with the considered one everytime is less than 90 meters.

Any ideas how to do this? Calculating the averages of the pixels that fall within the buffer is very easy, but I couldn't find a way to set any condition on them.

I will be immensly grateful for your help and suggestions. Approaches using other packages than satrs are also very welcome :)

Solution

  • Please see below a solution using terra. The code uses terra::extract to create two corresponding lists:

    • The pixel values
    • The surrounding buffer values

    Subsequently the values are processed pairwise, using mapply, with a function similar to the one you suggested.

    It's the first time I'm using terra but seems like terra::extract is much faster than raster::extract, therefore this solution may be feasible even for a large raster.

    Creating sample data:

    library(sf)
library(terra)

r = rast(ncol = ncol(volcano), nrow = nrow(volcano), xmin = 0, xmax = ncol(volcano), ymin = 0, ymax = nrow(volcano))
values(r) = volcano
s = r
s[] = rnorm(ncell(s))
r = c(r, s)
crs(r) = ""
plot(r)

    input

    Calculating buffers:

    pnt = as.points(r, values = FALSE)
pol = buffer(pnt, 10)

    Extracting raster values from points:

    x = extract(r, pnt)
head(x)

    ##      ID lyr.1       lyr.1
## [1,]  1   100 -0.03525223
## [2,]  2   100  0.31525467
## [3,]  3   101  0.94054608
## [4,]  4   101  0.37209238
## [5,]  5   101 -0.38388234
## [6,]  6   101 -0.03120593

    Extracting raster values from buffers:

    y = extract(r, pol)
head(y)

    ##      ID lyr.1       lyr.1
## [1,]  1   100 -0.03525223
## [2,]  1   100  0.31525467
## [3,]  1   101  0.94054608
## [4,]  1   101  0.37209238
## [5,]  1   101 -0.38388234
## [6,]  1   101 -0.03120593

    Now, the extracted values can be processed sequentially using mapply. First, we convert the objects to lists:

    x = as.data.frame(x)
x = split(x, x$ID)
y = as.data.frame(y)
y = split(y, y$ID)

    Next, we use mapply to make the necessary calculation, each time considering the current focal point value x and surrounding buffer values y:

    f = function(x, y) {
    d = abs(x[, 2] - y[, 2])  ## differences
    values = y[, 3]  ## values
    mean(values[d < 5], na.rm = TRUE)  ## Mean of subset
}
result = mapply(f, x, y)

    Finally, putting the results back into the raster template:

    u = r[[1]]
values(u) = result
plot(u)

    output