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

user_12585098 0
Viewed 443

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

r-stars raster spatial buffer r
1 Answers

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

user_13396260
Related