0

I would like to find the cells, and their corresponding centroid, that hold the maximum value for each band in a multiband raster.

Converting the stars-object to a sf-dataframe makes this easy somehow. But I also think this is super unefficient. And I have the strong feeling that this is an easy task which I can't figure out how to do... (using subsetting, st_apply, aggregate, ....)

This is the code (using the conversion to a sf-table)

library(stars)
library(sf)
tif = system.file("tif/L7_ETMs.tif", package = "stars")
x = read_stars(tif)
# for each cell find the geolocated-cell with the highest value 
df = st_as_sf(x)
max_per_band = df %>% 
  pivot_longer(
    cols = -geometry
  ) %>% 
  group_by(name) %>% 
  filter(value==max(value))

st_centroid(max_per_band)

            geometry name           value
 *        <POINT [m]> <chr>          <dbl>
 1   (288990 9119179) L7_ETMs.tif.V5   255
 2   (288990 9119179) L7_ETMs.tif.V6   255
 3 (297511.5 9118239) L7_ETMs.tif.V5   255
 4 (297511.5 9118239) L7_ETMs.tif.V6   255
 5   (296457 9117925) L7_ETMs.tif.V6   255
 6 (294376.5 9117127) L7_ETMs.tif.V5   255
 7 (294376.5 9117127) L7_ETMs.tif.V6   255
 8   (294348 9117099) L7_ETMs.tif.V1   255
 9   (294348 9117099) L7_ETMs.tif.V2   255
10   (294348 9117099) L7_ETMs.tif.V3   255
4
  • 1
    terra has the convenient where.max function which gets you 90% of the way there (returns cell numbers, but you can get coords when needed). You can convert stars to terra with terra::rast or work natively in terra?
    – Spacedman
    Commented Oct 31, 2023 at 17:28
  • sorry for the late question... Can I get the coords directly with that function or do I need to use the cell-index?
    – Lenn
    Commented Dec 4, 2023 at 11:33
  • I hesitated to make an answer with this because you are using stars in the Q and not terra. If terra is an option for you, how about making a new question and then we'll have solutions in stars (if the answer below is an answer, in which case it needs marking as accepted) and in terra....
    – Spacedman
    Commented Dec 4, 2023 at 12:05
  • ah yes of course! Sorry...
    – Lenn
    Commented Dec 4, 2023 at 12:24

1 Answer 1

2

I would take advantage of the fact that the representation of a stars object in R is an array so we can get the indices of cells where the band values are equal to the maximum value of the respective band. The code below then constructs an sf object by querying the coordinates of the cell centroids at the respective indices.

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
library(sf)

tif = system.file("tif/L7_ETMs.tif", package = "stars")
(x = read_stars(tif))
#> stars object with 3 dimensions and 1 attribute
#> attribute(s):
#>              Min. 1st Qu. Median     Mean 3rd Qu. Max.
#> L7_ETMs.tif     1      54     69 68.91242      86  255
#> dimension(s):
#>      from  to  offset delta                     refsys point x/y
#> x       1 349  288776  28.5 SIRGAS 2000 / UTM zone 25S FALSE [x]
#> y       1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y]
#> band    1   6      NA    NA                         NA    NA

x_dims <- st_get_dimension_values(x, which = 1, where = "center")
y_dims <- st_get_dimension_values(x, which = 2, where = "center")

band_max_indices <- apply(x[[1]], 3, function(x) which(x == max(x) , arr.ind = TRUE))

locs <- lapply(1:length(band_max_indices), function(i){
  band <- band_max_indices[[i]]
  m <- matrix(c(
    x_dims[band[ ,1]],
    y_dims[band[ ,2]]
  ), ncol = 2)
  
  p <- st_as_sf(as.data.frame(m), 
                coords = c(1,2), 
                crs = st_crs(x))
  
  p$band <- i
  p$x_index <- band[ ,2]
  p$y_index <- band[ ,1]
  p$value <- st_extract(x, p)[[1]][,i]
  p
})

(locs <- do.call(rbind, locs))
#> Simple feature collection with 61 features and 4 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 288990 ymin: 9110914 xmax: 297511.5 ymax: 9119179
#> Projected CRS: SIRGAS 2000 / UTM zone 25S
#> First 10 features:
#>    band x_index y_index value                 geometry
#> 1     1     129     196   255   POINT (294348 9117099)
#> 2     1     129     197   255 POINT (294376.5 9117099)
#> 3     1     318     182   255   POINT (293949 9111712)
#> 4     1     318     183   255 POINT (293977.5 9111712)
#> 5     1     319     180   255   POINT (293892 9111684)
#> 6     1     319     181   255 POINT (293920.5 9111684)
#> 7     1     319     182   255   POINT (293949 9111684)
#> 8     1     320     180   255   POINT (293892 9111655)
#> 9     1     320     181   255 POINT (293920.5 9111655)
#> 10    1     320     183   255 POINT (293977.5 9111655)

Created on 2023-11-01 with reprex v2.0.2

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.