NetCDF time aggregation using R stars package

Question

I have a monthly temperature netCDF dataset with aprox 20 years and I'd like to calculate a monthly climatology (average of all Jan, Feb, Mar...).

I know I can do this with stackapply from package raster. But I'm trying to learn the package stars. I tried using the aggregate function but I'm not able to figure out how to set the by argument correctly to apply the mean function to each month.

This is what my dataset looks like:

stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
 Tavg_monthly_UT_Brazil_v2_19800101_20151231.nc 
 Min.   :18.11                                  
 1st Qu.:25.00                                  
 Median :26.36                                  
 Mean   :25.85                                  
 3rd Qu.:26.85                                  
 Max.   :29.83                                  
 NA's   :56236                                  
dimension(s):
     from  to offset delta  refsys point                    values    
x       1 168 -74.25  0.25      NA    NA                      NULL [x]
y       1 162   6.25 -0.25      NA    NA                      NULL [y]
time    1 432     NA    NA POSIXct    NA 1980-01-16,...,2015-12-16    

Is aggregate the correct function to use? I also tried group_by but that does not work for stars objects.

Answer 1

I was able to calculate the monthly climatology. But I'm not sure it's the best approach.

monthly_mean <- function(month) {
  nc_time_series %>%
    filter(lubridate::month(time) == month) %>%
    st_apply(c(1,2), mean, na.rm = TRUE)
}

a <- lapply(1:12, monthly_mean) # calculate for each month

climate <- do.call(c, a) %>% # join list ellements
  merge() %>% # convert attribute to dimension
  setNames("Mean T") %>% # fix attribute name
  st_set_dimensions(3, values = 1:12) %>% # fix dimension size
  st_set_dimensions(names = c('x', 'y', 'month')) # set dimension name

Answer 2

I have added an example on how to aggregate along the time dimension in time steps. See this link to the stars-github-repo. For completeness I've copied it here:

# aggregate time dimension in format Date
tif = system.file("tif/L7_ETMs.tif", package = "stars")
t1 = as.Date("2018-07-31")
x = read_stars(c(tif, tif, tif, tif), along = list(time = c(t1, t1+1, t1+2, t1+3)), 
               RasterIO = list(nXOff = c(1), 
                               nYOff = c(1), 
                               nXSize = 50, 
                               nYSize = 50, 
                               bands = c(1:6)))
st_get_dimension_values(x, "time")
x_agg_time = aggregate(x, by = t1 + c(0, 2, 4), FUN = max)

# aggregate time dimension in format Date - interval
by_t = "2 days"
x_agg_time2 = aggregate(x, by = by_t, FUN = max)
st_get_dimension_values(x_agg_time2, "time")
x_agg_time - x_agg_time2

# aggregate time dimension in format POSIXct
x = st_set_dimensions(x, 4, values = as.POSIXct(c("2018-07-31",
                                                  "2018-08-01",
                                                  "2018-08-02",
                                                  "2018-08-03")),
                      names = "time")
by_t = as.POSIXct(c("2018-07-31", "2018-08-02"))
x_agg_posix = aggregate(x, by = by_t, FUN = max)
st_get_dimension_values(x_agg_posix, "time")
x_agg_time - x_agg_posix

Hope you can work off of this.

Added after comment:

# aggregate non continuous time dimension
x = read_stars(c(tif, tif, tif, tif), 
               along = list(time = c(t1, t1+10, t1+365, t1+700)), 
               RasterIO = list(nXOff = c(1), 
                               nYOff = c(1), 
                               nXSize = 50, 
                               nYSize = 50, 
                               bands = c(1:6)))
st_get_dimension_values(x, "time")
x_agg_time = aggregate(x, by = "1 year", FUN = max)
st_get_dimension_values(x_agg_time, "time")