2

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

khalil TEBER
  • 23
  • 4
  • 1
    Hi, sounds like an interesting question! Not sure if I understand "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" - could you please clarify what is the "considered one"? Do you mean to create a (circular?) buffer around a particular pixel, then take the mean from a subset of those pixels which are within the buffer, based on difference in elevation? – Michael Dorman Jan 19 '21 at 20:41
  • Yes exactly this is what I mean. I am trying with this function on one pixel: `temp_mean <- function(x,...){ buff_mean <- ifelse(abs(x[[2]][1] - x[[2]]) > 90,NA,x[[1]]) mean(buff_mean,na.rm=T) }` and then use it on the pixels selected with the buffer: `temp_mean(subset[st_buffer(subset2_sf,2000)])` But I can't generalize it on the raster. and I don't want to use a loop because then it will take forever with the whole rasters. – khalil TEBER Jan 19 '21 at 21:04
  • I understand, thanks! The loop approach you suggest is a good choice if the raster is not too big, so it is really a matter of data size and memory limits. This "loop" approach can be improved by (1) calculating all buffer polygons in advance, outside of the loop, and (2) using parallel processing. I have another idea, will try to develop it into an answer, meanwhile could you please share details what are the dimensions, resolution, and buffer size in the real application? Or even post an example of a raster? Thanks – Michael Dorman Jan 20 '21 at 07:34
  • The rasters I want to apply this on are the MODIS LST products for Bavaria for the years 2003 until 2019: `dimension(s): from to offset delta refsys point values x 1 408 4284425 860.124 +proj=tmerc +lat_0=0 +lon... FALSE NULL y 1 426 5604907 -860.358 +proj=tmerc +lat_0=0 +lon... FALSE NULL` so it's 173808 pixels per year/raster. I will upload the subset I am trying with now in the main post (96 pixels) – khalil TEBER Jan 20 '21 at 08:01
  • Thanks! Please see one approach below, you can test it on your rasters, first convert them to `terra` object with two bands (using function `rast` I guess), or reading directly from a file, then modify the buffer size (in the example it's `10`) and elevation threshold (in the example it's `5`). – Michael Dorman Jan 20 '21 at 08:18
  • P.S. 408x426 is not an extremely big raster, to keep the code simple I'd still try the "loop" approach just in case it does execute in a reasonable period of time... – Michael Dorman Jan 20 '21 at 08:21

1 Answers1

2

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

Michael Dorman
  • 950
  • 6
  • 8
  • 1
    Thank you very much, I think it works perfectly! I tried it on a bigger subset of the data (roughly 10% of the big raster) with a big buffer (47 km) and it worked rather fast (didn't benchmark it, but I think it's in the 20 minutes or so). I would have never had the idea of using Terra, and using it like you did :) – khalil TEBER Jan 20 '21 at 11:12