Skip to content

Commit 2387b51

Browse files
infotrophinfotroph
authored
Merge pull request #3641 from infotroph/rothc-coupler
RothC coupler
2 parents 09ea6c9 + 7ef8546 commit 2387b51

39 files changed

+1729
-1
lines changed

‎Makefile‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ NCPUS ?= 1
44
BASE := logger utils db settings visualization qaqc remote workflow
55

66
MODELS := basgra biocro clm45 dalec dvmdostem ed fates gday jules linkages \
7-
ldndc lpjguess maat maespa sibcasa sipnet stics template
7+
ldndc lpjguess maat maespa rothc sibcasa sipnet stics template
88

99
MODULES := allometry assim.batch assim.sequential benchmark \
1010
data.atmosphere data.land data.remote \

‎Makefile.depends‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@ $(call depends,models/linkages): | .install/base/db .install/base/logger .instal
2121
$(call depends,models/lpjguess): | .install/base/logger .install/base/remote .install/base/utils
2222
$(call depends,models/maat): | .install/base/logger .install/base/remote .install/base/settings .install/base/utils .install/modules/data.atmosphere
2323
$(call depends,models/maespa): | .install/base/logger .install/base/remote .install/base/utils .install/modules/data.atmosphere
24+
$(call depends,models/rothc): | .install/base/logger .install/base/utils
2425
$(call depends,models/sibcasa): | .install/base/logger
2526
$(call depends,models/sipnet): | .install/base/logger .install/base/remote .install/base/utils .install/modules/data.atmosphere .install/modules/data.land
2627
$(call depends,models/stics): | .install/base/logger .install/base/remote .install/base/settings .install/base/utils

‎docker/depends/pecan_package_dependencies.csv‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -55,6 +55,7 @@
5555
"dplyr","*","models/biocro","Imports",FALSE
5656
"dplyr","*","models/ed","Imports",FALSE
5757
"dplyr","*","models/ldndc","Imports",FALSE
58+
"dplyr","*","models/rothc","Imports",FALSE
5859
"dplyr","*","models/sipnet","Imports",FALSE
5960
"dplyr","*","models/stics","Imports",FALSE
6061
"dplyr","*","modules/assim.sequential","Imports",FALSE
@@ -171,6 +172,7 @@
171172
"lubridate",">= 1.6.0","models/lpjguess","Imports",FALSE
172173
"lubridate",">= 1.6.0","models/maat","Imports",FALSE
173174
"lubridate",">= 1.6.0","models/maespa","Imports",FALSE
175+
"lubridate",">= 1.6.0","models/rothc","Imports",FALSE
174176
"lubridate",">= 1.6.0","models/sipnet","Imports",FALSE
175177
"lubridate",">= 1.6.0","modules/assim.batch","Imports",FALSE
176178
"lubridate",">= 1.6.0","modules/assim.sequential","Imports",FALSE
@@ -233,6 +235,7 @@
233235
"ncdf4","*","models/basgra","Imports",FALSE
234236
"ncdf4","*","models/dvmdostem","Imports",FALSE
235237
"ncdf4","*","models/ldndc","Imports",FALSE
238+
"ncdf4","*","models/rothc","Imports",FALSE
236239
"ncdf4","*","models/sibcasa","Imports",FALSE
237240
"ncdf4","*","models/stics","Imports",FALSE
238241
"ncdf4","*","modules/assim.sequential","Imports",FALSE
@@ -340,6 +343,7 @@
340343
"PEcAn.logger","*","models/lpjguess","Imports",TRUE
341344
"PEcAn.logger","*","models/maat","Imports",TRUE
342345
"PEcAn.logger","*","models/maespa","Imports",TRUE
346+
"PEcAn.logger","*","models/rothc","Imports",TRUE
343347
"PEcAn.logger","*","models/sibcasa","Imports",TRUE
344348
"PEcAn.logger","*","models/sipnet","Imports",TRUE
345349
"PEcAn.logger","*","models/stics","Imports",TRUE
@@ -433,6 +437,7 @@
433437
"PEcAn.utils",">= 1.4.8","models/ldndc","Imports",TRUE
434438
"PEcAn.utils",">= 1.4.8","models/stics","Imports",TRUE
435439
"PEcAn.utils",">= 1.4.8","models/template","Imports",TRUE
440+
"PEcAn.utils",">= 1.8.0","models/rothc","Imports",TRUE
436441
"PEcAn.visualization","*","modules/assim.sequential","Suggests",TRUE
437442
"PEcAn.visualization","*","modules/data.land","Imports",TRUE
438443
"PEcAn.visualization","*","modules/priors","Suggests",TRUE
@@ -497,6 +502,7 @@
497502
"rlang","*","models/biocro","Imports",FALSE
498503
"rlang","*","models/ed","Imports",FALSE
499504
"rlang","*","models/ldndc","Imports",FALSE
505+
"rlang","*","models/rothc","Imports",FALSE
500506
"rlang","*","models/sipnet","Imports",FALSE
501507
"rlang","*","modules/assim.sequential","Imports",FALSE
502508
"rlang","*","modules/benchmark","Imports",FALSE
@@ -539,6 +545,7 @@
539545
"roxygen2","== 7.3.2","models/lpjguess","Roxygen",FALSE
540546
"roxygen2","== 7.3.2","models/maat","Roxygen",FALSE
541547
"roxygen2","== 7.3.2","models/maespa","Roxygen",FALSE
548+
"roxygen2","== 7.3.2","models/rothc","Roxygen",FALSE
542549
"roxygen2","== 7.3.2","models/sibcasa","Roxygen",FALSE
543550
"roxygen2","== 7.3.2","models/sipnet","Roxygen",FALSE
544551
"roxygen2","== 7.3.2","models/stics","Roxygen",FALSE
@@ -574,6 +581,7 @@
574581
"sp","*","modules/data.land","Imports",FALSE
575582
"sp","*","modules/data.remote","Imports",FALSE
576583
"stats","*","base/qaqc","Imports",FALSE
584+
"stats","*","models/rothc","Imports",FALSE
577585
"stats","*","models/sipnet","Imports",FALSE
578586
"stats","*","modules/allometry","Imports",FALSE
579587
"stats","*","modules/assim.batch","Imports",FALSE
@@ -632,6 +640,7 @@
632640
"testthat",">= 2.0.0","base/utils","Suggests",FALSE
633641
"testthat",">= 2.0.0","models/biocro","Suggests",FALSE
634642
"testthat",">= 2.0.0","modules/benchmark","Suggests",FALSE
643+
"testthat",">= 3.0.0","models/rothc","Suggests",FALSE
635644
"testthat",">= 3.0.0","models/sibcasa","Suggests",FALSE
636645
"testthat",">= 3.0.4","base/qaqc","Suggests",FALSE
637646
"testthat",">= 3.1.0","modules/data.land","Suggests",FALSE
@@ -667,6 +676,7 @@
667676
"utils","*","models/ed","Imports",FALSE
668677
"utils","*","models/linkages","Imports",FALSE
669678
"utils","*","models/lpjguess","Imports",FALSE
679+
"utils","*","models/rothc","Imports",FALSE
670680
"utils","*","modules/allometry","Imports",FALSE
671681
"utils","*","modules/assim.batch","Imports",FALSE
672682
"utils","*","modules/assim.sequential","Suggests",FALSE
@@ -685,6 +695,7 @@
685695
"withr","*","base/workflow","Suggests",FALSE
686696
"withr","*","models/basgra","Suggests",FALSE
687697
"withr","*","models/ed","Suggests",FALSE
698+
"withr","*","models/rothc","Suggests",FALSE
688699
"withr","*","models/sibcasa","Suggests",FALSE
689700
"withr","*","models/sipnet","Suggests",FALSE
690701
"withr","*","modules/allometry","Suggests",FALSE

‎models/rothc/.Rbuildignore‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
Dockerfile
2+
model_info.json
3+
inst/workflow_example

‎models/rothc/DESCRIPTION‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
Package: PEcAn.RothC
2+
Type: Package
3+
Title: PEcAn Package for Integration of the RothC Model
4+
Version: 0.0.0.9000
5+
Authors@R: c(
6+
person("Chris", "Black", role = c("aut", "cre"), email = "chris@ckblack.org"),
7+
person("Rothamsted Research", role = c("cph"))
8+
)
9+
Description: This module provides functions to link the Rothamstead soil carbon model "RothC" to PEcAn.
10+
Uses RothC >= 2.0, available on GitHub <https://github.com/Rothamsted-Models/RothC_Code>
11+
as well as directly from the RothC team <https://www.rothamsted.ac.uk/rothamsted-carbon-model-rothc>.
12+
URL: https://pecanproject.github.io
13+
BugReports: https://github.com/PecanProject/pecan/issues
14+
Depends: R (>= 4.1)
15+
Imports:
16+
dplyr,
17+
lubridate (>= 1.6.0),
18+
ncdf4,
19+
PEcAn.logger,
20+
PEcAn.utils (>= 1.8.0),
21+
rlang,
22+
stats,
23+
utils
24+
Suggests:
25+
testthat (>= 3.0.0),
26+
withr
27+
SystemRequirements: RothC
28+
OS_type: unix
29+
License: BSD_3_clause + file LICENSE
30+
Copyright: Authors
31+
Encoding: UTF-8
32+
RoxygenNote: 7.3.2
33+
Config/testthat/edition: 3

‎models/rothc/Dockerfile‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
# this needs to be at the top, what version are we building
2+
ARG IMAGE_VERSION="latest"
3+
4+
5+
# ----------------------------------------------------------------------
6+
# BUILD PECAN FOR MODEL
7+
# ----------------------------------------------------------------------
8+
FROM pecan/models:${IMAGE_VERSION}
9+
10+
# ----------------------------------------------------------------------
11+
# INSTALL MODEL SPECIFIC PIECES
12+
# ----------------------------------------------------------------------
13+
14+
#RUN apt-get update \
15+
# && apt-get install -y --no-install-recommends \
16+
# python3.9 \
17+
# && rm -rf /var/lib/apt/lists/*
18+
19+
# ----------------------------------------------------------------------
20+
# SETUP FOR SPECIFIC MODEL
21+
# ----------------------------------------------------------------------
22+
23+
# Some variables that can be used to set control the docker build
24+
ARG MODEL_VERSION=2.1.0
25+
26+
RUN mkdir -p /tmp/rothc \
27+
&& curl -sSL https://github.com/Rothamsted-Models/RothC_Code/archive/refs/tags/v${MODEL_VERSION}.tar.gz \
28+
| tar xzf - -C /tmp/rothc \
29+
&& cd /tmp/rothc/RothC_Code-${MODEL_VERSION} \
30+
&& gfortran RothC.for Shell.for -o rothc_bin \
31+
&& mv rothc_bin /usr/local/bin/RothC_v${MODEL_VERSION} \
32+
&& rm -rf /tmp/rothc
33+
34+
# TODO hard-coding this path is probably wrong
35+
# Setup model_info file
36+
# @VERSION@ is replaced with model version in the model_info.json file
37+
# @BINARY@ is replaced with model binary in the model_info.json file
38+
COPY model_info.json /work/model.json
39+
RUN sed -i -e "s/@VERSION@/${MODEL_VERSION}/g" \
40+
-e "s#@BINARY@#/usr/local/bin/RothC_v${MODEL_VERSION}#g" /work/model.json

‎models/rothc/LICENSE‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
YEAR: 2024
2+
COPYRIGHT HOLDER: PEcAn Project
3+
ORGANIZATION: PEcAn Project, authors affiliations

‎models/rothc/NAMESPACE‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
# Generated by roxygen2: do not edit by hand
2+
3+
export(met2model.RothC)
4+
export(model2netcdf.RothC)
5+
export(read_restart.RothC)
6+
export(write.config.RothC)
7+
export(write_restart.RothC)
8+
importFrom(rlang,.data)
9+
importFrom(rlang,.env)

‎models/rothc/NEWS.md‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
# PEcAn.RothC 0.0.0.9000
2+
3+
Initial development version

‎models/rothc/R/met2model.RothC.R‎

Lines changed: 185 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,185 @@
1+
#' Extract monthly weather from CF file for input to RothC
2+
#'
3+
#' Input files need to be named `<in.path>/<in.prefix>.YYYY.nc`
4+
#'
5+
#' Output files are named `<outfolder>/<in.prefix>.YY-mm.YY-mm.dat`
6+
#' with one line per month and columns for temperature, rainfall,
7+
#' and evaporation.
8+
#'
9+
#' Note that the created file contains only weather data and not any of the
10+
#' soil or management data needed for RothC's single combined input file.
11+
#' See `write.config.RothC()` for assembly into a model-ready RothC_input.dat`.
12+
#'
13+
#' @param in.path path on disk where CF files live
14+
#' @param in.prefix prefix for each file
15+
#' @param outfolder location where model specific output is written.
16+
#' @param start_date,end_date When to start and end output.
17+
#' Specify as exact dates, but output will be padded to whole months.
18+
#' @param overwrite logical: replace output files if they already exist?
19+
#' @return data frame summarizing file metadata
20+
#' @export
21+
#' @author Chris Black
22+
met2model.RothC <- function(in.path,
23+
in.prefix,
24+
outfolder,
25+
start_date,
26+
end_date,
27+
overwrite = FALSE) {
28+
29+
PEcAn.logger::logger.info("START met2model.RothC")
30+
31+
start_date <- as.Date(start_date)
32+
end_date <- as.Date(end_date)
33+
start_year <- strftime(start_date, "%Y")
34+
end_year <- strftime(end_date, "%Y")
35+
year_regex <- paste(start_year:end_year, collapse = "|")
36+
37+
if (grepl("\\.nc$", in.prefix)) {
38+
# Assume it's the full filename rather than a prefix
39+
# NB also means we assume it contains the whole requested date range
40+
name_pattern <- in.prefix
41+
} else {
42+
name_pattern <- paste0(in.prefix, "\\.(", year_regex, ")\\.nc$")
43+
}
44+
45+
nc_files <- list.files(in.path, pattern = name_pattern, full.names = TRUE)
46+
47+
if (length(nc_files) == 0) {
48+
PEcAn.logger::logger.severe(
49+
"No files found matching ", in.prefix,
50+
"for years", start_year, ":", end_year,
51+
"; cannot process data."
52+
)
53+
}
54+
55+
# TODO complain (fail?) here if some but not all years found
56+
57+
out_filename <- paste(
58+
in.prefix,
59+
strftime(start_date, "%Y-%m"),
60+
strftime(end_date, "%Y-%m"),
61+
"dat",
62+
sep = "."
63+
)
64+
out_path <- file.path(outfolder, out_filename)
65+
results <- data.frame(file = out_path,
66+
host = Sys.getenv(
67+
"FQDN",
68+
unset = Sys.info()[["nodename"]]
69+
),
70+
mimetype = "text/tab-separated-values",
71+
formatname = "RothC.dat",
72+
startdate = start_date,
73+
enddate = end_date,
74+
dbfile.name = out_filename,
75+
stringsAsFactors = FALSE)
76+
PEcAn.logger::logger.info("internal results")
77+
PEcAn.logger::logger.info(results)
78+
79+
if (file.exists(out_path) && !overwrite) {
80+
PEcAn.logger::logger.debug(
81+
"File '", out_path, "' already exists, skipping to next file."
82+
)
83+
return(invisible(results))
84+
}
85+
86+
if (!file.exists(outfolder)) {
87+
dir.create(outfolder)
88+
}
89+
90+
met <- nc_files |>
91+
lapply(
92+
read_nc,
93+
varnames = c("air_temperature", "precipitation_flux", "specific_humidity")
94+
) |>
95+
do.call(what = "rbind")
96+
97+
# TODO probably need more care with partial months here:
98+
# we check if data extends to start/end, but not whether that includes enough
99+
# days to treat as a whole month.
100+
# e.g. if start_date = YYYY-01-31 and data starts YYYY-01-30 ->
101+
# current code will aggregate those two days as if they were all of January.
102+
# Consider failing if >n days missing in any output month?
103+
first_month <- lubridate::floor_date(start_date, unit = "month")
104+
last_month <- lubridate::ceiling_date(end_date, unit = "month")
105+
met <- met[(met$timestamp >= first_month) & (met$timestamp < last_month), ]
106+
if (as.Date(min(met$timestamp)) > start_date
107+
|| as.Date(max(met$timestamp)) < end_date) {
108+
PEcAn.logger::logger.severe(
109+
"input (",
110+
paste(range(met$timestamp), collapse = " to "),
111+
") does not cover requested time window (",
112+
start_date, "to", end_date, ")"
113+
)
114+
}
115+
116+
timestep <- met$timestamp |>
117+
diff(units = "secs") |>
118+
mean() |>
119+
as.numeric()
120+
121+
met$year <- lubridate::year(met$timestamp)
122+
met$month <- lubridate::month(met$timestamp)
123+
met$Tmp_C <- met$air_temperature |>
124+
PEcAn.utils::ud_convert("K", "degC")
125+
met$Rain_mm <- met$precipitation_flux * timestep # kg/m2/sec -> mm total
126+
met$Evap_mm <- 0# TODO... how to convert Qair to pan evaporation?
127+
128+
met_monthly <- merge(
129+
stats::aggregate(met, Tmp_C ~ year + month, mean),
130+
stats::aggregate(met, cbind(Rain_mm, Evap_mm) ~ year + month, sum),
131+
sort = FALSE # would treat months as strings; sort as numbers below instead
132+
)
133+
met_monthly <- met_monthly[order(met_monthly$year, met_monthly$month), ]
134+
135+
utils::write.table(
136+
# as.data.frame to write integer columns as ints not floats
137+
x = format(as.data.frame(met_monthly), digits = 4),
138+
file = out_path,
139+
quote = FALSE,
140+
sep = "\t",
141+
row.names = FALSE,
142+
col.names = TRUE
143+
)
144+
145+
results
146+
}
147+
148+
149+
150+
151+
152+
# slurp named vars from one PEcAn nc into a dataframe with timestamp
153+
#
154+
# TODO could read other dimensions if present too, but consider if worth it --
155+
# maybe this function is better for files where only the time dimension varies
156+
#
157+
# if vars = NULL, read all of them
158+
read_nc <- function(ncfile, varnames = NULL) {
159+
160+
nc <- ncdf4::nc_open(ncfile)
161+
on.exit(ncdf4::nc_close(nc), add = TRUE)
162+
163+
timestamps <- PEcAn.utils::cf2datetime(
164+
nc$dim$time$vals,
165+
nc$dim$time$units
166+
)
167+
168+
if (is.null(varnames)) {
169+
varnames <- names(nc$var)
170+
}
171+
172+
var_values <- lapply(
173+
varnames,
174+
ncdf4::ncvar_get,
175+
nc = nc
176+
)
177+
178+
# todo handle this case (multi-loc files?)
179+
stopifnot(all(sapply(var_values, length) == length(timestamps)))
180+
181+
var_values |>
182+
stats::setNames(varnames) |>
183+
as.data.frame() |>
184+
transform(timestamp = timestamps)
185+
}

0 commit comments

Comments
 (0)