Behaviour of `aggregate.stars` when multiple overlaps occur in the same raster pixel

[agila5]

Dear all, I created this issue since I noticed that when I run aggregate.stars and one of the raster cells cover more than one polygon, then the aggregation process occurs only for one of the polygons. For example (see the NA at the end after I conver the output to sf object):

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE

m <- matrix(1, nrow = 1, ncol = 1)
dim(m) <- c(x = 1, y = 1)
s <- st_as_stars(m, dimensions = st_dimensions(x = 0:1, y = 0:1))
image(s)

p1 <- st_polygon(
  x = list(
    rbind(
      c(0.2, 0.2), 
      c(0.4, 0.2), 
      c(0.4, 0.4), 
      c(0.2, 0.4), 
      c(0.2, 0.2)
    )
  )
)

p2 <- p1 + c(0.4, 0.4)
polys <- st_sfc(p1, p2) 
plot(polys, add = TRUE)

aggregate(s, polys, FUN = mean, as_points = FALSE) |> st_as_sf()
#> Simple feature collection with 2 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 0.2 ymin: 0.2 xmax: 0.8 ymax: 0.8
#> CRS:           NA
#>    X                       geometry
#> 1  1 POLYGON ((0.2 0.2, 0.4 0.2,...
#> 2 NA POLYGON ((0.6 0.6, 0.8 0.6,...

^{Created on 2024-02-17 with reprex v2.0.2}

Is this a bug or expected behaviour? I notice that there are some comments regarding this on the .R files:

https://github.com/r-spatial/stars/blob/20f55f8162cfb4c4c308f600ddf7ee6d6ba9b7a0/R/aggregate.R#L133-L138

but I really don't get it...

[edzer]

To aggregate: a whole formed by combining several separate elements.

The question is what s is: is it "one", one pixel, one observation, or is it an area that has the value 1 everywhere? If it is one (and pixels are often thought of as being small, and singular, point-like), how can it be assigned to polys[1] when it already has been assigned to polys[2] (or vice versa) by combining?

If not expected, it is in any case documented (see the docs of ?aggregate.stars, argument by). You may get what you want by either

aggregate(s, polys, exact=TRUE)

or

st_interpolate_aw(s, polys, extensive=FALSE)

[edzer]

or

st_extract(s, st_centroid(polys)) |> st_as_sf()

[agila5]

Dear @edzer, thank you very much for your quick reply, your comments, and your suggestions!

I agree with you that aggregate() should only be used to aggregate/combine pixels into (larger) elements, and the examples I presented here are ambiguous wrt to the assignment of pixels to polygons. I asked this question since, for my current project, I need a flexible approach that permits upscaling/downscaling of pixels to larger/smaller spatial units and I hoped I could simply adopt aggregate(raster, polygons, ...). Unfortunately, this is not the case, but I will develop some code starting from aggregate()!

aggregate(s, polys, exact=TRUE)

btw, for some reason I currently don't understand, I get only NA when applying this approach. This is not relevant here but maybe I will open another issue as soon as I can easily showcase the problem.

[edzer]

here:

> aggregate(s, polys, exact=TRUE, FUN=mean) |> st_as_sf()
Simple feature collection with 2 features and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 0.2 ymin: 0.2 xmax: 0.8 ymax: 0.8
CRS:           NA
  X                       geometry
1 1 POLYGON ((0.2 0.2, 0.4 0.2,...
2 1 POLYGON ((0.6 0.6, 0.8 0.6,...

maybe reinstall exactextractr?

[agila5]

No, I mean I get NAs when applying that code to the "real" data (which are quite difficult to share...) Working on a reprex and, hopefully, will share here ASAP 😄

[agila5]

Just kidding, I simply forgot na.rm, sorry 🤦‍♂️🤦‍♂️🤦‍♂️