Temporal rolling mean in stars

Question

I need to find temporal rolling means (or other stats - max, min, etc.) in stars spatial obects. For stars raster cubes, this would be directly analogous to the calc(rasterBrick, function(x) movingFun(x, 3, mean)) available in the raster package. I could just use raster instead of stars if all I had were raster cubes. However, I'm really liking the ability of stars to handle vector data as well, and I have some multipolygon data with time slices in stars objects that I would also like to get rolling statistics on. For example, I might try to calculate the 3-day rolling mean for the values in each polygon.

I can use something like aggregate(polygonStars, by = 'weeks', FUN = mean, na.rm = TRUE) to get the mean over discrete time units, but aggregate aggregates rather than rolls. I've tried using st_apply, but haven't had success with feeding it a rolling function like RcppRoll::roll_mean (I'm unclear what exactly it's doing when I feed it a rolling function since it wants a single output, but it ends up replacing a dimension instead of the attribute values).

I've tried just shoving the objects into the rolling functions from RcppRoll and zoo, but they expect vectors, and I haven't found a way to tell them to operate along the 3rd dimension. If that's possible, I suppose I could take the stars objects apart into vectors, do the rolling calculations, and then put them back together, but that would be cumbersome, error prone, and slow.

Ideally, I'd remain in stars, but if there's not the possibility of getting rolling temporal means from stars objects, I'll switch to raster for the rasters, and seek a suggestion for the polygons with temporal layers.

To illustrate using the example from https://r-spatial.github.io/stars/, How do I get the rolling means for each pixel in the raster cube prec and each polygon in the vector cube a, as indicated by the last lines of this code (which don't work because I'm using functions from raster to show what I'm looking for)?

# Modified from the data raster and vector time series example at
# https://r-spatial.github.io/stars/ to illustrate

library(stars)
library(tidyverse)

# Get the rasters: these look similar to what my data has
prec_file = system.file("nc/test_stageiv_xyt.nc", package = "stars")
(prec = read_ncdf(prec_file, curvilinear = c("lon", "lat"), ignore_bounds = TRUE))

# Get north carolina
sf::read_sf(system.file("gpkg/nc.gpkg", package = "sf"), "nc.gpkg") %>%
  st_transform(st_crs(prec)) -> nc # transform from NAD27 to WGS84
nc_outline = st_union(st_geometry(nc))

# Do the aggregate to make polygons with temporal layers
a = aggregate(prec, by = nc, FUN = max)

# Rolling means over time 
# These don't work because the functions are from `raster`, but this is the sort
# of thing I want to do.
rastRoll <- calc(prec, function(x) movingFun(x, 3, mean))
vectRoll <- calc(a, function(x) movingFun(x, 3, mean))

Answer 1

With stars of raster geometry, this is straightforward using st_apply. Any function, including raster::movingFun, can be applied on per-pixel basis, as follows:

rastRoll = st_apply(prec, c("x", "y"), function(x) raster::movingFun(x, 3, mean))
plot(rastRoll)

With stars of vector geometry, I couldn't find a way to use st_apply. However, the following workaround applies raster::movingFun directly on the matrix of stars layer values:

vectRoll = a
vectRoll[[1]] = t(apply(vectRoll[[1]], 1, function(x) raster::movingFun(x, 3, mean)))
plot(vectRoll, col = hcl.colors(12))