Weighted aggregation of stars object

I was wondering if it's possible to adjust the following code in such a way that the aggregate() function takes into account the overlapping area between each pixel defined in a stars object and the "cells" used for calculating the overlap.

For example:

# load packages
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

# define two grids starting from a general polygon
gpoly <- st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 1), c(0, 0)))))
grid1 <- st_sf(data.frame(x = 1:9), geometry = st_make_grid(gpoly, n = c(3, 3)))
grid2 <- st_make_grid(gpoly, n = c(2, 2))

# convert the first grid into stars format and plot them. 
st1 <- st_as_stars(grid1)
plot(st1, reset = FALSE)
plot(st_boundary(grid2), col = "red", add = TRUE, lwd = 3)

Now I would like to aggregate the values defined by st1 using grid2 (i.e. the red cells in the figure above). The following code performs such an aggregation, but each new value is computed by taking a simple average of the values of the pixels intersecting a given cell (since the merge is performed using join = st_intersects).

aggregate(
  st1, 
  by = grid2, 
  FUN = mean 
) 
#> stars object with 1 dimensions and 1 attribute
#> attribute(s):
#>    Min. 1st Qu. Median Mean 3rd Qu. Max.
#> x     3   4.125      6    6   7.875    9
#> dimension(s):
#>          from to point
#> geometry    1  4 FALSE
#>                                                                 values
#> geometry POLYGON ((0 0, 0.5 0, 0.5...,...,POLYGON ((0.5 0.5, 1 0.5,...

Is it possible to adjust the previous code such that the means are computed taking into account the area that overlaps pixels and cells? Like a weighted mean where the weights are the areas of overlap.

^{Created on 2023-10-25 with reprex v2.0.2}

Solution

You should use st_interpolate_aw, as in

st_interpolate_aw(st1, grid2, extensive = FALSE)