-
Notifications
You must be signed in to change notification settings - Fork 1
/
check_the_numbers.R
66 lines (57 loc) · 1.89 KB
/
check_the_numbers.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
library(tidync)
library(ncdf4)
library(tidyverse)
library(lubridate)
library(units)
install_unit("plant", def="unitless")
e_file <- "/data/output/pecan_runs/MANDIFORE_big_run/MANDIFORE-SEUS-352/out/SA-median/analysis-E-2005-02-00-000000-g01.h5"
# -E- file ----------------------------------------------------------------
# look up metadata
e_nc <- nc_open(e_file)
ncatt_get(e_nc, "MMEAN_NPPDAILY_CO")
ncatt_get(e_nc, "MMEAN_NPPDAILY_PY")
ncatt_get(e_nc, "NPLANT")
ncatt_get(e_nc, "PFT")
ncatt_get(e_nc, "PACO_N")
ncatt_get(e_nc, "SIPA_N")
ncatt_get(e_nc, "PYSI_N")
# Extract cohort-level variables from E file
e_cohort <-
tidync(e_file) |>
activate("D0") |>
hyper_tibble(select_var = c(
"MMEAN_NPPDAILY_CO", #"Monthly mean - Net primary productivity - total [kgC/m2/yr]" but actual units are kgC/plant/yr https://github.com/EDmodel/ED2/issues/342
"PFT", #ED2 PFT number
"NPLANT" # "Plant density [plant/m2]"
)) |>
mutate(
MMEAN_NPPDAILY_CO = set_units(MMEAN_NPPDAILY_CO, "kg/plant/yr"),
NPLANT = set_units(NPLANT, "plant/m^2")
)
# Extract polygon-level variables from E file
e_polygon <-
tidync(e_file) |>
activate("D4")|>
hyper_tibble(select_var = c(
"AREA",
"PACO_ID", # "First index for patch
"PACO_N" # Cohort count in each patch
))
# rolling-join (not actually necessary in this single-patch example)
e_df <-
left_join(e_cohort, e_polygon, join_by(closest(phony_dim_0 >= PACO_ID)))
e_df
# Polygon-level total
target <-
tidync(e_file) |>
activate("D1") |>
hyper_tibble(select_var = c(
"MMEAN_NPPDAILY_PY"#, # "Monthly mean - Net primary productivity - total [kgC/m2/yr]"
)) |>
mutate(
MMEAN_NPPDAILY_PY = set_units(MMEAN_NPPDAILY_PY, "kg/m^2/yr")
) |> pull(MMEAN_NPPDAILY_PY)
# does it all add up?
tot <- e_df |>
summarize(npp_total = sum(MMEAN_NPPDAILY_CO*NPLANT*AREA)) |> pull(npp_total)
waldo::compare(tot, target, tolerance = 0.01)