NetCDF aggregation using {terra}

aggregate_netcdf.R

library(terra)
library(stringr)

#---- static code ----

# read in the original data
r <- rast("LPX-Bern_SH1_gpppft.nc")

# first subset to only the required years
# 1978 - 2020
loc <- which(time(r) >= "1978-01-01")
r <- subset(r, loc)

#--- the part below is dynamic ----

# subset the data further based on PFT
# we need some additional information for this

# create index to aggregate by
# based on the layer names

# first dump the raw names in a data frame
df <- data.frame(names = names(r))

# split into parts and clean up using
# string substitution (keep all parts for reference)
df[,c("product", "PFT","index")] <- stringr::str_split_fixed(df$names, "_", 3)
df$PFT <- as.numeric(gsub("PFT=", "", df$PFT))

# add a month and year index
df$month <- format(time(r), "%m")
df$year <- format(time(r), "%Y")

# these are the locations of the
# PFT selection
pft_loc <- which(df$PFT <= 8) # range from 1-8 but can be adjusted

#---- subset for a particular PFT range ----

# subset the original raster r with
# this index, keep the original and
# create a new object
s <- subset(r, pft_loc)

# also subset the meta-data
df <- df[pft_loc,]

#---- sum by year-month across PFTs ----

year_month_idx <- paste(df$year, df$month)
year_month_total <- tapp(
  s,
  year_month_idx,
  fun = sum,
  na.rm = TRUE
  )

#---- create monthly average across years ----

month_idx <- data.frame(names = names(year_month_total))
month_idx[,c("year", "month")] <- stringr::str_split_fixed(month_idx$names, "\\.", 2)

month_mean <- tapp(
  year_month_total,
  month_idx$month,
  fun = mean,
  na.rm = TRUE
)

#---- plot monthly mean values summed across the first 8 PFT ----

plot(month_mean)