1

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.

2

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
|improve this answer|||||
0

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")
|improve this answer|||||
  • Cool, I was not aware I could do that. But how about aggregating non-continuous intervals. For example, Jan-1980 + Jan-1981 + Jan-1982. Can it be done? – Daniel Dec 4 '19 at 19:45
  • I updated the answer. I think this is what you are trying to do. – przell Dec 5 '19 at 8:52
  • Thanks. But I'm not sure if I did not understand your answer or if I did not explain well what I was trying to do. Say I have 10 years of monthly data. I want to end up with only 12 images, the first being the average off all January images (all years). The second is the mean off all February and so on. – Daniel Dec 6 '19 at 18:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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