Created
November 29, 2022 15:19
-
-
Save khufkens/6f828e051a3942d01261aa8886034aa0 to your computer and use it in GitHub Desktop.
NetCDF aggregation using {terra}
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment