|
| 1 | +#!/usr/bin/env Rscript |
| 2 | + |
| 3 | +# --- Profiling Start --- |
| 4 | +# Rprof("profiling.out") |
| 5 | +# --- End Profiling Start --- |
| 6 | + |
| 7 | +# Minimal MVP: build mvp_events.json from ca_field_attributes.csv |
| 8 | +# - Input: data/ca_field_attributes.csv (columns: site_id, year, pft, ...) |
| 9 | +# - Output: data/mvp_events.json following data/pecan_events_schema_v0.1.0.json |
| 10 | +# - Events (minimal): |
| 11 | +# * planting: annual crops -> every site-year; woody perennials -> first observed year only |
| 12 | +# * harvest: all site-years |
| 13 | +# Each event includes only the schema-required fields per event_type. |
| 14 | + |
| 15 | +# --- Config --- |
| 16 | +data_dir <- "/projectnb2/dietzelab/ccmmf/data" |
| 17 | +field_attr_csv <- file.path(data_dir, "ca_field_attributes.csv") |
| 18 | +sample_output_json <- file.path(data_dir, "events/mvp_events.json") |
| 19 | +output_json <- file.path(data_dir, "events/events.json") |
| 20 | + |
| 21 | +# if TRUE, only generate for design points |
| 22 | +# TODO: generate full set for all sites to use in site selection and downscaling |
| 23 | +DESIGN_POINTS <- TRUE |
| 24 | + |
| 25 | +PRODUCTION <- FALSE # set TRUE for all sites, not needed if DESIGN_POINTS is TRUE |
| 26 | +if (PRODUCTION) { |
| 27 | + stop("This could be very slow; consider profiling and writing to db or arrow etc") |
| 28 | +} |
| 29 | +set.seed(123) |
| 30 | + |
| 31 | +ca_field_attributes <- vroom::vroom( |
| 32 | +field_attr_csv, |
| 33 | + show_col_types = FALSE |
| 34 | +) |
| 35 | + |
| 36 | +if (DESIGN_POINTS) { |
| 37 | + # design_points <- readr::read_csv("https://raw.githubusercontent.com/ccmmf/workflows/refs/heads/main/data/design_points.csv") |
| 38 | + # d <- update_design_point_site_ids(design_points, ca_field_attributes) |
| 39 | + # readr::write_csv(d, file.path(data_dir, "design_points.csv")) |
| 40 | + # readr::write_csv(d, "~/downscaling/data/design_points.csv") |
| 41 | + # design_points <- readr::read_csv(file.path(data_dir, "design_points.csv")) |
| 42 | + # use the one under version control |
| 43 | + design_points <- readr::read_csv("~/downscaling/data/design_points.csv") |
| 44 | + ca_field_attributes <- ca_field_attributes |> |
| 45 | + dplyr::filter(site_id %in% design_points$site_id) |
| 46 | +} else if (!PRODUCTION) { |
| 47 | + ca_field_attributes <- ca_field_attributes |> |
| 48 | + dplyr::slice_sample(n = 1000) |
| 49 | +} |
| 50 | + |
| 51 | +ca_fields <- ca_field_attributes |> |
| 52 | + dplyr::select(site_id, pft, crop) |> |
| 53 | + dplyr::distinct() |> |
| 54 | + tidyr::crossing(year = 2016:2024) |> |
| 55 | + dplyr::group_by(site_id) |> |
| 56 | + dplyr::mutate(first_year = min(year)) |> |
| 57 | + dplyr::ungroup() |
| 58 | + |
| 59 | +# Planting (annuals) |
| 60 | +planting_annual <- ca_fields |> |
| 61 | + dplyr::filter(pft == "annual crop") |> |
| 62 | + dplyr::transmute( |
| 63 | + event_type = "planting", |
| 64 | + date = paste0(year, "-03-15"), |
| 65 | + site_id = site_id, |
| 66 | + # required for planting |
| 67 | + leaf_c_kg_m2 = 0.05, |
| 68 | + crop = crop |
| 69 | + ) |
| 70 | + |
| 71 | +# Planting (woody): first year |
| 72 | +planting_woody <- ca_fields |> |
| 73 | + dplyr::filter(pft == "woody perennial crop") |> |
| 74 | + dplyr::filter(year == first_year) |> |
| 75 | + dplyr::transmute( |
| 76 | + event_type = "planting", |
| 77 | + date = paste0(year, "-03-15"), |
| 78 | + site_id = site_id, |
| 79 | + leaf_c_kg_m2 = 0.2, |
| 80 | + crop = crop |
| 81 | + ) |
| 82 | + |
| 83 | +# Fertilization |
| 84 | +fertilization <- ca_fields |> |
| 85 | + dplyr::transmute( |
| 86 | + event_type = "fertilization", |
| 87 | + date = paste0(year, "-02-11"), |
| 88 | + site_id = site_id, |
| 89 | + org_n_kg_m2 = 0.0, |
| 90 | + org_c_kg_m2 = 0.0, |
| 91 | + nh4_n_kg_m2 = 0.02, |
| 92 | + no3_n_kg_m2 = 0.03 |
| 93 | + ) |
| 94 | + |
| 95 | +# Organic Matter Addition |
| 96 | +organic_matter_addition <- ca_fields |> |
| 97 | + dplyr::transmute( |
| 98 | + event_type = "fertilization", |
| 99 | + date = paste0(year, "-03-11"), |
| 100 | + site_id = site_id, |
| 101 | + org_n_kg_m2 = 0.05, |
| 102 | + org_c_kg_m2 = 0.5, |
| 103 | + nh4_n_kg_m2 = 0.0, |
| 104 | + no3_n_kg_m2 = 0.0 |
| 105 | + ) |
| 106 | + |
| 107 | +# Harvest |
| 108 | +harvest <- ca_fields |> |
| 109 | + dplyr::transmute( |
| 110 | + event_type = "harvest", |
| 111 | + date = paste0(year, "-10-15"), |
| 112 | + site_id = site_id, |
| 113 | + frac_above_removed_0to1 = 0.10, |
| 114 | + frac_below_removed_0to1 = 0.0, |
| 115 | + frac_above_to_litter_0to1 = 0.0, |
| 116 | + frac_below_to_litter_0to1 = 0.0, |
| 117 | + crop = crop |
| 118 | + ) |
| 119 | + |
| 120 | +# Pruning (woody) |
| 121 | +pruning <- ca_fields |> |
| 122 | + dplyr::filter(pft == "woody perennial crop") |> |
| 123 | + dplyr::mutate(offset = year - first_year) |> |
| 124 | + dplyr::filter(offset %% 4 == 1) |> |
| 125 | + dplyr::transmute( |
| 126 | + event_type = "harvest", |
| 127 | + date = paste0(year, "-12-15"), |
| 128 | + site_id = site_id, |
| 129 | + frac_above_removed_0to1 = 0.30, |
| 130 | + frac_below_removed_0to1 = 0.0, |
| 131 | + frac_above_to_litter_0to1 = 0.0, |
| 132 | + frac_below_to_litter_0to1 = 0.0, |
| 133 | + crop = crop |
| 134 | + ) |
| 135 | + |
| 136 | +# Tillage |
| 137 | +tillage <- ca_fields |> |
| 138 | + dplyr::filter(pft == "annual crop") |> |
| 139 | + tidyr::crossing(till_suffix = c("-03-01", "-11-01")) |> |
| 140 | + dplyr::transmute( |
| 141 | + event_type = "tillage", |
| 142 | + date = paste0(year, till_suffix), |
| 143 | + site_id = site_id, |
| 144 | + tillage_eff_0to1 = 0.10 |
| 145 | + ) |
| 146 | + |
| 147 | +# Irrigation (both pfts): 3 per month for all months |
| 148 | +# TODO: Should annual crops skip irrigation during fallow season? |
| 149 | +months <- sprintf("%02d", 1:12) |
| 150 | +days <- c("05", "15", "25") |
| 151 | + |
| 152 | +irrigation <- ca_fields |> |
| 153 | + tidyr::crossing(month = months, day = days) |> |
| 154 | + dplyr::transmute( |
| 155 | + event_type = "irrigation", |
| 156 | + date = paste0(year, "-", month, "-", day), |
| 157 | + site_id = site_id, |
| 158 | + amount_mm = 40, |
| 159 | + method = "soil" |
| 160 | + ) |
| 161 | + |
| 162 | +# Combine and order by site/date |
| 163 | +events_all <- dplyr::bind_rows( |
| 164 | + planting_annual, planting_woody, |
| 165 | + harvest, pruning, |
| 166 | + tillage, irrigation, |
| 167 | + fertilization, organic_matter_addition |
| 168 | +) |> |
| 169 | + dplyr::arrange(site_id, date) |
| 170 | + |
| 171 | +# --- Build site objects per schema ------------------------------------------ |
| 172 | +# Helper: drop NULL/NA fields from a named list |
| 173 | +compact_list <- function(x) { |
| 174 | + Filter(function(v) !(is.null(v) || (length(v) == 1 && is.atomic(v) && is.na(v))), x) |
| 175 | +} |
| 176 | + |
| 177 | +sites <- unique(events_all$site_id) |
| 178 | + |
| 179 | +site_objs <- purrr::map(sites, function(sid) { |
| 180 | + evs_df <- events_all |> |
| 181 | + dplyr::filter(site_id == sid) |> |
| 182 | + dplyr::arrange(date) |
| 183 | + |
| 184 | + # Only include required fields for each event type |
| 185 | + evs_list <- purrr::pmap( |
| 186 | + evs_df, |
| 187 | + function(event_type, date, site_id, leaf_c_kg_m2 = NA_real_, frac_above_removed_0to1 = NA_real_, |
| 188 | + frac_below_removed_0to1 = NA_real_, frac_above_to_litter_0to1 = NA_real_, |
| 189 | + frac_below_to_litter_0to1 = NA_real_, amount_mm = NA_real_, method = NA_character_, |
| 190 | + tillage_eff_0to1 = NA_real_, org_c_kg_m2 = NA_real_, org_n_kg_m2 = NA_real_, |
| 191 | + nh4_n_kg_m2 = NA_real_, no3_n_kg_m2 = NA_real_, |
| 192 | + crop = NA_character_, ...) { |
| 193 | + base <- list(event_type = event_type, date = date) |
| 194 | + |
| 195 | + # Add required fields per event type |
| 196 | + if (event_type == "planting" && !is.na(leaf_c_kg_m2)) { |
| 197 | + base$leaf_c_kg_m2 <- leaf_c_kg_m2 |
| 198 | + if (!is.na(crop)) base$crop <- crop |
| 199 | + } |
| 200 | + if (event_type == "harvest" && !is.na(frac_above_removed_0to1)) { |
| 201 | + base$frac_above_removed_0to1 <- frac_above_removed_0to1 |
| 202 | + if (!is.na(frac_below_removed_0to1)) base$frac_below_removed_0to1 <- frac_below_removed_0to1 |
| 203 | + if (!is.na(frac_above_to_litter_0to1)) base$frac_above_to_litter_0to1 <- frac_above_to_litter_0to1 |
| 204 | + if (!is.na(frac_below_to_litter_0to1)) base$frac_below_to_litter_0to1 <- frac_below_to_litter_0to1 |
| 205 | + if (!is.na(crop)) base$crop <- crop |
| 206 | + } |
| 207 | + if (event_type == "irrigation" && !is.na(amount_mm) && !is.na(method)) { |
| 208 | + base$amount_mm <- amount_mm |
| 209 | + base$method <- method |
| 210 | + } |
| 211 | + if (event_type == "tillage" && !is.na(tillage_eff_0to1)) { |
| 212 | + base$tillage_eff_0to1 <- tillage_eff_0to1 |
| 213 | + } |
| 214 | + if (event_type == "fertilization" && !is.na(org_c_kg_m2)) { |
| 215 | + base$org_c_kg_m2 <- org_c_kg_m2 |
| 216 | + if (!is.na(org_n_kg_m2)) base$org_n_kg_m2 <- org_n_kg_m2 |
| 217 | + } |
| 218 | + |
| 219 | + compact_list(base) |
| 220 | + } |
| 221 | + ) |
| 222 | + list( |
| 223 | + pecan_events_version = "0.1.0", |
| 224 | + site_id = sid, |
| 225 | + events = evs_list |
| 226 | + ) |
| 227 | +}) |
| 228 | + |
| 229 | +# TODO add PEcAn Schema info |
| 230 | + |
| 231 | +# Validate JSON given schema |
| 232 | +# schema <- "data/pecan_events_schema_v0.1.0.json" |
| 233 | +# validator <- jsonvalidate::json_validator(schema) |
| 234 | +# json_txt_temp <- jsonlite::toJSON(site_objs, auto_unbox = TRUE) |
| 235 | +# if (!validator(json_txt_temp)) { |
| 236 | +# stop("JSON does not match schema") |
| 237 | +# } |
| 238 | + |
| 239 | +# --- Write JSON -------------------------------------------------------------- |
| 240 | + |
| 241 | +# Complete |
| 242 | +jsonlite::write_json(site_objs, path = output_json, pretty = FALSE, auto_unbox = TRUE) |
| 243 | +# Single site example |
| 244 | +jsonlite::write_json(site_objs[1:3], path = gsub(".json", "_3sites.json", output_json), pretty = TRUE, auto_unbox = TRUE) |
| 245 | +# When dealing with full dataset, may need to write to more performant files |
| 246 | +# #Sample |
| 247 | +# jsonlite::write_json(site_objs[1:100], path = sample_output_json, pretty = TRUE, auto_unbox = TRUE) |
| 248 | + |
| 249 | +# # Complete - compressed |
| 250 | +output_json_gz <- paste0(output_json, ".gz") |
| 251 | +gz_con <- gzfile(output_json_gz, "w") |
| 252 | +jsonlite::write_json(site_objs, path = gz_con, pretty = FALSE, auto_unbox = TRUE) |
| 253 | +close(gz_con) |
| 254 | + |
| 255 | +# --- Profiling End --- |
| 256 | +# Rprof(NULL) |
| 257 | +# summaryRprof("profiling.out") |
| 258 | +# --- End Profiling End --- |
0 commit comments