Skip to content

LPJ-GUESS restart support: new write_restart, read_restart tweaks & CRU driver hooks #3533

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mdietze merged 3 commits into from
May 31, 2025
Merged

LPJ-GUESS restart support: new write_restart, read_restart tweaks & CRU driver hooks #3533

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a43309c
LPJ-GUESS: add functional write_restart + update read_restart & confi…
yinghaoSunn May 30, 2025
97c11f7
Delete dplyr %>% in write_restart.LPJGUESS
yinghaoSunn May 30, 2025
5b99c00
Fix some issues.
yinghaoSunn May 30, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions models/lpjguess/NAMESPACE
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,14 @@ export(met2model.LPJGUESS)
export(model2netcdf.LPJGUESS)
export(pecan2lpjguess)
export(readStateBinary)
export(read_binary_LPJGUESS)
export(read_restart.LPJGUESS)
export(split_inputs.LPJGUESS)
export(update_state_LPJGUESS)
export(write.config.LPJGUESS)
export(write.insfile.LPJGUESS)
export(write_binary_LPJGUESS)
export(write_restart.LPJGUESS)
importFrom(PEcAn.utils,days_in_year)
importFrom(Rcpp,sourceCpp)
importFrom(ncdf4,nc_close)
Expand Down
40 changes: 26 additions & 14 deletions models/lpjguess/R/read_restart.LPJGUESS.R
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,13 +1,25 @@

# developing
# outdir = "/fs/data2/output//PEcAn_1000010473/out"
# runid = 1002656839
# stop.time = "1960-12-31 23:59:59 UTC"
# load("/fs/data2/output/PEcAn_1000010473/SDAsettings_develop.Rdata")
# var.names = c("AGB.pft", "TotSoilCarb")
# load("/fs/data2/output/PEcAn_1000010473/SDAparams_develop.Rdata")


#' Read Restart for LPJGUESS
#'
#' @param outdir output directory
#' @param runid run ID
#' @param stop.time year that is being read
#' @param settings PEcAn settings object
#' @param var.names var.names to be extracted
#' @param params passed on to return value
#'
#' @return X_tmp vector of forecasts
#' @export
#' @examples
#' \dontrun{
#' rx <- read_restart.LPJGUESS(
#' outdir = "/projectnb/…/LPJ_output",
#' runid = "123456",
#' stop.time = as.POSIXct("2001-12-31 23:59:59", tz = "UTC"),
#' settings = settings,
#' var.names = c("AGB.pft"),
#' params = params)
#' }
#' @author Istem Fer, Yinghao Sun
read_restart.LPJGUESS <- function(outdir, runid, stop.time, settings, var.names, params){

# which LPJ-GUESS version, the structure of state file depends a lot on version
Expand All @@ -24,7 +36,6 @@ read_restart.LPJGUESS <- function(outdir, runid, stop.time, settings, var.names,
# read binary state file, takes a couple of minutes
Gridcell_container <- read_binary_LPJGUESS(outdir = file.path(outdir, runid),
version = lpjguess_ver)

forecast <- list()

# additional varnames for LPJ-GUESS?
Expand All @@ -33,8 +44,8 @@ read_restart.LPJGUESS <- function(outdir, runid, stop.time, settings, var.names,

if (var_name == "AGB.pft") {

cmass_sap_perpft <- calculateGridcellVariablePerPFT(model.state = Gridcell_container, variable = "cmass_sap")
cmass_heart_perpft <- calculateGridcellVariablePerPFT(model.state = Gridcell_container, variable = "cmass_heart")
cmass_sap_perpft <- calculateGridcellVariablePerPFT(model.state = Gridcell_container$state, variable = "cmass_sap")
cmass_heart_perpft <- calculateGridcellVariablePerPFT(model.state = Gridcell_container$state, variable = "cmass_heart")

cmass_wood <- cmass_sap_perpft + cmass_heart_perpft
cmass_wood <- PEcAn.utils::ud_convert(cmass_wood, "kg/m^2", "Mg/ha")
Expand All @@ -45,11 +56,12 @@ read_restart.LPJGUESS <- function(outdir, runid, stop.time, settings, var.names,
cmass_abvg_wood <- cmass_wood - cmass_blwg_wood

forecast[[length(forecast) + 1]] <- cmass_abvg_wood
names(forecast[[length(forecast)]]) <- paste0("AGB.pft.", unlist(Gridcell_container$meta_data$pft))
names(forecast[[length(forecast)]]) <- paste0("AGB.pft.", unlist(Gridcell_container$state$meta_data$pft))

}
}

# params$LPJGUESS_state include state, pos_list, siz_list
params$LPJGUESS_state <- Gridcell_container

PEcAn.logger::logger.info("Finished --", runid)
Expand Down
201 changes: 165 additions & 36 deletions models/lpjguess/R/read_state.R
View file
Open in desktop

Large diffs are not rendered by default.

233 changes: 119 additions & 114 deletions models/lpjguess/R/split_inputs.LPJGUESS.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,120 +14,125 @@
##' @importFrom PEcAn.utils days_in_year
##' @export
split_inputs.LPJGUESS <- function(settings, start.time, stop.time, inputs, overwrite = FALSE, outpath = NULL){

#### Lubridate start and end times
start.day <- lubridate::yday(start.time)
start.year <- lubridate::year(start.time)
end.day <- lubridate::yday(stop.time)
end.year <- lubridate::year(stop.time)

# Whole run period
run.start <- lubridate::year(settings$run$start.date)
run.end <- lubridate::year(settings$run$end.date)

#### Get met paths
met <- inputs
path <- dirname(met)
prefix <- substr(basename(met), 1, nchar(basename(met))-16) #assuming we'll always have "PREFIX.1920.2010.tmp"
if(is.null(outpath)){
outpath <- path
}
if(!dir.exists(outpath)) dir.create(outpath)

var.names <- c("tmp", "pre", "cld")
long.names <- c("air_temperature",
"precipitation_flux",
"surface_downwelling_shortwave_flux_in_air")

# !!! always full years with LPJ-GUESS !!!
files.in <- file.path(outpath, paste0(prefix, run.start, ".", run.end, ".", var.names, ".nc"))
files.out <- file.path(outpath, paste0(prefix, start.year, ".", end.year, ".", var.names, ".nc"))

if(file.exists(files.out[1]) & !overwrite){
return(files.out[1])
}

## open netcdf files
fnc.tmp <- ncdf4::nc_open(files.in[1])
fnc.pre <- ncdf4::nc_open(files.in[2])
fnc.cld <- ncdf4::nc_open(files.in[3])

## read climate data
nc.tmp <- ncdf4::ncvar_get(fnc.tmp, var.names[1])
nc.pre <- ncdf4::ncvar_get(fnc.pre, var.names[2])
nc.cld <- ncdf4::ncvar_get(fnc.cld, var.names[3])

# cut where
if(start.year == run.start){
years <- start.year:end.year
inds <- 1:sum(PEcAn.utils::days_in_year(years))
#### If using CRU input, return directly
if (grepl("\\.cru(\\.bin)?$", inputs, ignore.case = TRUE)) {
PEcAn.logger::logger.info(paste("Input is a CRU file:", inputs, "- returning path directly without splitting."))
return(inputs) # Without cropping, use the original file directly
}else{
### come back
}

# split
nc.tmp <- nc.tmp[1,1,inds]
nc.pre <- nc.pre[1,1,inds]
nc.cld <- nc.cld[1,1,inds]

var.list <- list(nc.tmp, nc.pre, nc.cld)

# not that these will be different than "K", "kg m-2 s-1", "W m-2"
var.units <- c(fnc.tmp$var$tmp$units,
fnc.pre$var$pre$units,
fnc.cld$var$cld$units)

# get other stuff to be written to ncdf

## retrieve lat/lon
lon <- ncdf4::ncvar_get(fnc.tmp, "lon")
lat <- ncdf4::ncvar_get(fnc.tmp, "lat")

# write back
## write climate data define dimensions

latdim <- ncdf4::ncdim_def(name = "lat", "degrees_north", as.double(lat))
londim <- ncdf4::ncdim_def(name = "lon", "degrees_east", as.double(lon))
timedim <- ncdf4::ncdim_def("time", paste0("days since ", start.year - 1, "-12-31", sep = ""), as.double(c(1:length(nc.tmp))))

fillvalue <- 9.96920996838687e+36

for (n in seq_along(var.names)) {
# define variable
var.def <- ncdf4::ncvar_def(name = var.names[n],
units = var.units[n],
dim = (list(londim, latdim, timedim)),
fillvalue, long.names[n],
verbose = FALSE,
prec = "float")

# create netCD file for LPJ-GUESS
ncfile <- ncdf4::nc_create(files.out[[n]], vars = var.def, force_v4 = TRUE)


# put variable, rep(...,each=4) is a hack to write the same data for all grids (which all are the
# same)
ncdf4::ncvar_put(ncfile, var.def, rep(var.list[[n]], each = 4))


# additional attributes for LPJ-GUESS
ncdf4::ncatt_put(nc = ncfile, varid = var.names[n], attname = "standard_name", long.names[n])

ncdf4::ncatt_put(nc = ncfile, varid = "lon", attname = "axis", "X")
ncdf4::ncatt_put(nc = ncfile, varid = "lon", attname = "standard_name", "longitude")

ncdf4::ncatt_put(nc = ncfile, varid = "lat", attname = "axis", "Y")
ncdf4::ncatt_put(nc = ncfile, varid = "lat", attname = "standard_name", "latitude")

ncdf4::ncatt_put(nc = ncfile, varid = "time", attname = "calendar", "gregorian")

ncdf4::nc_close(ncfile)
#### Lubridate start and end times
start.day <- lubridate::yday(start.time)
start.year <- lubridate::year(start.time)
end.day <- lubridate::yday(stop.time)
end.year <- lubridate::year(stop.time)

# Whole run period
run.start <- lubridate::year(settings$run$start.date)
run.end <- lubridate::year(settings$run$end.date)

#### Get met paths
met <- inputs
path <- dirname(met)
prefix <- substr(basename(met), 1, nchar(basename(met))-16) #assuming we'll always have "PREFIX.1920.2010.tmp"
if(is.null(outpath)){
outpath <- path
}
if(!dir.exists(outpath)) dir.create(outpath)

var.names <- c("tmp", "pre", "cld")
long.names <- c("air_temperature",
"precipitation_flux",
"surface_downwelling_shortwave_flux_in_air")

# !!! always full years with LPJ-GUESS !!!
files.in <- file.path(outpath, paste0(prefix, run.start, ".", run.end, ".", var.names, ".nc"))
files.out <- file.path(outpath, paste0(prefix, start.year, ".", end.year, ".", var.names, ".nc"))

if(file.exists(files.out[1]) & !overwrite){
return(files.out[1])
}

## open netcdf files
fnc.tmp <- ncdf4::nc_open(files.in[1])
fnc.pre <- ncdf4::nc_open(files.in[2])
fnc.cld <- ncdf4::nc_open(files.in[3])

## read climate data
nc.tmp <- ncdf4::ncvar_get(fnc.tmp, var.names[1])
nc.pre <- ncdf4::ncvar_get(fnc.pre, var.names[2])
nc.cld <- ncdf4::ncvar_get(fnc.cld, var.names[3])

# cut where
if(start.year == run.start){
years <- start.year:end.year
inds <- 1:sum(PEcAn.utils::days_in_year(years))
}else{
### come back
}

# split
nc.tmp <- nc.tmp[1,1,inds]
nc.pre <- nc.pre[1,1,inds]
nc.cld <- nc.cld[1,1,inds]

var.list <- list(nc.tmp, nc.pre, nc.cld)

# not that these will be different than "K", "kg m-2 s-1", "W m-2"
var.units <- c(fnc.tmp$var$tmp$units,
fnc.pre$var$pre$units,
fnc.cld$var$cld$units)

# get other stuff to be written to ncdf

## retrieve lat/lon
lon <- ncdf4::ncvar_get(fnc.tmp, "lon")
lat <- ncdf4::ncvar_get(fnc.tmp, "lat")

# write back
## write climate data define dimensions

latdim <- ncdf4::ncdim_def(name = "lat", "degrees_north", as.double(lat))
londim <- ncdf4::ncdim_def(name = "lon", "degrees_east", as.double(lon))
timedim <- ncdf4::ncdim_def("time", paste0("days since ", start.year - 1, "-12-31", sep = ""), as.double(c(1:length(nc.tmp))))

fillvalue <- 9.96920996838687e+36

for (n in seq_along(var.names)) {
# define variable
var.def <- ncdf4::ncvar_def(name = var.names[n],
units = var.units[n],
dim = (list(londim, latdim, timedim)),
fillvalue, long.names[n],
verbose = FALSE,
prec = "float")

# create netCD file for LPJ-GUESS
ncfile <- ncdf4::nc_create(files.out[[n]], vars = var.def, force_v4 = TRUE)


# put variable, rep(...,each=4) is a hack to write the same data for all grids (which all are the
# same)
ncdf4::ncvar_put(ncfile, var.def, rep(var.list[[n]], each = 4))


# additional attributes for LPJ-GUESS
ncdf4::ncatt_put(nc = ncfile, varid = var.names[n], attname = "standard_name", long.names[n])

ncdf4::ncatt_put(nc = ncfile, varid = "lon", attname = "axis", "X")
ncdf4::ncatt_put(nc = ncfile, varid = "lon", attname = "standard_name", "longitude")

ncdf4::ncatt_put(nc = ncfile, varid = "lat", attname = "axis", "Y")
ncdf4::ncatt_put(nc = ncfile, varid = "lat", attname = "standard_name", "latitude")

ncdf4::ncatt_put(nc = ncfile, varid = "time", attname = "calendar", "gregorian")

ncdf4::nc_close(ncfile)
}

# close nc
ncdf4::nc_close(fnc.tmp)
ncdf4::nc_close(fnc.pre)
ncdf4::nc_close(fnc.cld)

return(files.out[1])
}

# close nc
ncdf4::nc_close(fnc.tmp)
ncdf4::nc_close(fnc.pre)
ncdf4::nc_close(fnc.cld)

return(files.out[1])
} # split_inputs.LPJGUESS
34 changes: 20 additions & 14 deletions models/lpjguess/R/write.config.LPJGUESS.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@
##' @export
##' @author Istem Fer, Tony Gardella
write.config.LPJGUESS <- function(defaults, trait.values, settings, run.id, restart = NULL) {

# find out where to write run/ouput
rundir <- file.path(settings$host$rundir, run.id)
if (!file.exists(rundir)) {
Expand Down Expand Up @@ -194,16 +193,20 @@ write.insfile.LPJGUESS <- function(settings, trait.values, rundir, outdir, run.i
paramsins <- paramsins[-pftindx]
paramsins <- c(paramsins, unlist(write2pftblock))


# write clim file names

tmp.file <- settings$run$inputs$met$path
pre.file <- gsub(".tmp.nc", ".pre.nc", tmp.file)
cld.file <- gsub(".tmp.nc", ".cld.nc", tmp.file)

guessins <- gsub("@TEMP_FILE@", tmp.file, guessins)
guessins <- gsub("@PREC_FILE@", pre.file, guessins)
guessins <- gsub("@INSOL_FILE@", cld.file, guessins)
# # Past version: write clim file names (cf input)
# tmp.file <- settings$run$inputs$met$path
# pre.file <- gsub(".tmp.nc", ".pre.nc", tmp.file)
# cld.file <- gsub(".tmp.nc", ".cld.nc", tmp.file)
#
# guessins <- gsub("@TEMP_FILE@", tmp.file, guessins)
# guessins <- gsub("@PREC_FILE@", pre.file, guessins)
# guessins <- gsub("@INSOL_FILE@", cld.file, guessins)

# when using cru input, lpjguess will not use these clim files
cru.file <- settings$run$inputs$met$path
misc.file <- sub("\\.bin$", "misc.bin", cru.file)
guessins <- gsub("@MET_AND_SOIL_FILE@", cru.file, guessins)
guessins <- gsub("@MISC_FILE@", misc.file, guessins)

# create and write CO2 file
start.year <- lubridate::year(settings$run$start.date)
Expand Down Expand Up @@ -232,9 +235,12 @@ write.insfile.LPJGUESS <- function(settings, trait.values, rundir, outdir, run.i
utils::write.table(CO2, file = co2.file, row.names = FALSE, col.names = FALSE, sep = "\t", eol = "\n")
guessins <- gsub("@CO2_FILE@", co2.file, guessins)

# write soil file path
soil.file <- settings$run$inputs$soil$path
guessins <- gsub("@SOIL_FILE@", soil.file, guessins)
# # write soil file path
# # when using cru input, it's also climate file
# soil.file <- settings$run$inputs$soil$path
Copy link
Member

@mdietze mdietze May 30, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Similarly, curious why the code setting the soils data is commented out. Is this just temporary? Presumably we'll want PEcAn to have the ability to pass in site-specific soil information.

Copy link
Contributor Author

@yinghaoSunn yinghaoSunn May 30, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, it's temporary. We’ve switched to CRU input, and the .cru.bin file actually provides both met and soil data.

# misc.file <- sub("\\.bin$", "misc.bin", soil.file)
# guessins <- gsub("@SOIL_FILE@", soil.file, guessins)
# guessins <- gsub("@MISC_FILE@", misc.file, guessins)

settings$model$insfile <- file.path(settings$rundir, run.id, "guess.ins")

Expand Down
Loading
Loading