Add irrigation projection workflow

modules/data.remote/inst/ccmmf/documentation/projections/06-irrigation_predictions.md

@@ -0,0 +1,124 @@
+# Irrigation Projections
+
+## Overview
+
+This section is split into two scripts. `caladapt_climate.R` turns downloaded Cal-Adapt NetCDFs into a
+county daily climate lookup, and `irrigation_projection.R` runs the projection on that series using
+the projected crop and date products.
+
+The water balance is the same as the inventory, but everything going into it changes.
+Weather comes from downscaled climate projections rather than observations, planting and harvest dates
+come from the date projections, and peak greenness that the inventory observes has to be
+reconstructed.
+
+For simplicity, the prediction outputs are specific to GCM CESM2/ssp370, and the inventory resolves
+precipitation and ET to the parcel level, while the projection averages 45 km Cal-Adapt grid cells to
+one daily series per county. Parcel-level variation survives in crop type, soil water holding capacity,
+rooting depth, and phenology, but not in the weather.
+
+## Caladapt_climate.R
+
+Overall, this script reads two types of files, each generated by an earlier workflow and whose locations need to be specified in the `config.yml`:
+
+a. Cal-Adapt WRF NetCDFs, one per grid cell per year, named `<GCM>.<SSP>.<year>.nc`. Cells sit in
+directories named by their grid position (`x47_y52`), and the current download covers California at
+45 km, which comes to 198 cells.
+b. California county boundaries, pulled from `tigris`
+
+Assigning cells to counties: each grid cell directory holds a latitude and longitude for its centre,
+which are turned into points and intersected with county polygons. At 45 km a small county can contain
+no cell centre at all, so any county without one is assigned its nearest cell by projected distance.
+Every county therefore ends up with at least one. The lookup is written out alongside the climate file.
+
+Calculating reference ET: the hourly NetCDF variables — air temperature, wind speed, specific humidity,
+pressure, shortwave and longwave radiation, precipitation — are aggregated to daily, then daily ET0 is
+computed with the FAO-56 Penman-Monteith method. This is the same algorithm CIMIS uses, so the
+projected ET0 is on the same footing as the inventory's observed values.
+
+County aggregation: daily ET0 and precipitation are averaged across the cells assigned to each county,
+giving one series per county per day. Large counties average several cells; small ones are represented
+by a single cell, sometimes one whose centre falls outside the county.
+
+**Important:** the year range has to extend one year past the final projection year. A crop planted
+late in 2045 is harvested in 2046, and the water balance runs to the harvest date.
+
+Processing is parallel across files, which is the slow part of this script. The output is a single CSV
+of county daily climate.
+
+## Irrigation_projection.R
+
+Overall, this script reads six types of files to create the irrigation predictions:
+
+a. Projected crop identity files, one per scenario and year
+b. Projected planting and harvest files
+c. `assigned_year=<year>_gapfilled.parquet` from the matched directory, for historical peak positions
+d. The county daily climate file from `caladapt_climate.R`
+e. SSURGO weights and soil tables
+f. Crop water parameters from PEcAn's data.land module
+
+Assembling crop cycles: projected crops join to their planting dates on crop code. Harvest carries no
+crop code, so cycles are matched chronologically — both tables are sorted by their own date column and
+paired by position within each parcel-year. A count mismatch between the two is a hard stop.
+
+Reconstructing peak greenness: the water balance needs a peak date to shape the canopy curve, and
+nothing observes one in the future. Instead the inventory phenology record gives where the peak fell as
+a fraction of the season — the gap from onset to peak over the gap from onset to end — and that
+fraction is applied to the projected season length. Matching previous scripts, this has five fallback levels:
+1. county + crop code
+2. county + crop class
+3. crop code
+4. crop class
+5. statewide
+
+Crop water parameters: each crop maps to a BISM crop name, which supplies the crop coefficient curve,
+the management allowable depletion fraction, and the rooting depth. Where LandIQ codes have no exact
+BISM match a small proxy table assigns the nearest reasonable crop. Parcels with no mapping, no rooting
+depth, or no depletion fraction are dropped, and the counts by reason are printed.
+
+Soil water holding capacity: effective capacity is computed per parcel and rooting depth from SSURGO
+horizons, weighted across components and map units. Results are cached per county, so the first run
+pays the cost and later runs reuse it. Rice is exempt, since the flooded model does not use a bucket
+capacity.
+
+Running the balance: counties run in parallel with 6 workers, with the largest ones first. Within a county, parcels sharing a crop,
+planting date, peak date, and harvest date are grouped so the ET series is computed once for the group, slightly improving run time.
+Each county writes its own file on completion, so an interrupted run resumes by skipping counties that
+already have output.
+
+**Important: simplifications**
+
+1. Daily weather is uniform within a county, and comes from 45 km cells. This is the main trade off for reasonable run/computational times.
+
+2. The canopy curve is piecewise linear — 0.15 at planting, rising to 1 at peak, falling back toward
+0.15 at harvest.
+
+3. Peak greenness is derived rather than observed, from where the peak historically fell within the
+season for that crop and county.
+
+4. Unclassified fallow and idle land receive no irrigation.
+
+5. Parcels with no usable crop water parameters or no SSURGO capacity are dropped rather than given a
+fallback, so this product covers fewer parcels than the others.
+
+6. County files under `_county_tmp` are reused on a rerun, so you can pick up where you left off on irrigation prediction without having to restart. 
+
+Any change to crop cycles, dates, or soil inputs requires deleting them, or the run will mix old and new results.
+
+7. The script stops rather than carrying on if a crop has no planting or harvest date, if a harvest does
+not follow its planting, if a peak date is impossible, if a climate day is missing, or if an irrigation
+amount in the output is not positive.
+
+## Outputs
+
+`caladapt_climate.R` writes to your indvidually defined `work_root`: 
+
+`caladapt_county_daily_climate.csv` — County, date, GCM, SSP, ET0_mm, precip_mm, and daily temperatures
+`caladapt_county_site_lookup.csv` — which grid cell supplies each county
+
+`irrigation_projection.R` writes to scenario subdirectories under the irrigation output root, matching
+the inventory structure/column names:
+
+`<scenario>/irrigation_statewide_<2024-2045>.parquet` — event_type, parcel_id, date, amount_mm
+
+Per-county intermediate files remain under `<scenario>/_county_tmp/` after a run and are what makes
+resuming possible. They are not a deliverable.

modules/data.remote/inst/get_climate_data.R

@@ -0,0 +1,250 @@
+pacman::p_load(ncdf4, sf, tidyr, dplyr, fs, purrr, furrr, future, lubridate, tigris, FAO56, parallelly)
+
+# ---- setup ----
+config = config::get(config = "scc", file = "config.yml")
+output_dir = config$work_root
+base_dir = config$wrf_base_dir
+
+# Default config uses a relative extracted directory under work_root
+if (!grepl("^/", base_dir)) {
+  base_dir = file.path(config$work_root, base_dir)
+}
+
+years_keep = seq.int(config$start_year + 1L, config$end_year + 1L)
+
+gcm_keep = config$climate_gcm
+ssp_keep = config$climate_ssp
+
+dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
+
+# ---- helpers ----
+specific_humidity_to_ea = function(q, pressure_pa) {
+  P_kpa = pressure_pa / 1000
+  (q * P_kpa) / (0.622 + 0.378 * q)
+}
+
+calculate_net_radiation = function(sw_rad, lw_rad, temp_c, albedo = 0.23, emissivity = 0.98) {
+  sigma = 5.67e-8
+  Rns = (1 - albedo) * sw_rad
+  Rnl = lw_rad - emissivity * sigma * (temp_c + 273.15)^4
+  (Rns + Rnl) * 0.0036
+}
+
+parse_nc_time = function(time_steps, time_units, year_val) {
+  origin_string = sub("^.*since\\s+", "", time_units)
+  origin_time = lubridate::ymd_hms(origin_string, tz = "UTC", quiet = TRUE)
+
+  if (is.na(origin_time)) {
+    origin_time = lubridate::ymd_hms(sprintf("%d-01-01 00:00:00", year_val), tz = "UTC")
+  }
+
+  if (grepl("^seconds", time_units, ignore.case = TRUE)) {
+    origin_time + lubridate::seconds(time_steps)
+  } else if (grepl("^hours", time_units, ignore.case = TRUE)) {
+    origin_time + lubridate::hours(time_steps)
+  } else if (grepl("^days", time_units, ignore.case = TRUE)) {
+    origin_time + lubridate::days(time_steps)
+  } else {
+    stop("Unknown NetCDF time units: ", time_units)
+  }
+}
+
+calculate_daily_et0 = function(df_hourly, wind_height_m) {
+  df_daily = df_hourly |>
+    dplyr::mutate(
+      ea_kpa = specific_humidity_to_ea(.data$spec_hum, .data$pressure),
+      Rn_MJ = calculate_net_radiation(.data$sw_rad, .data$lw_rad, .data$temp_c),
+      date = as.Date(.data$datetime)
+    ) |>
+    dplyr::group_by(.data$date) |>
+    dplyr::summarise(
+      T_min = min(.data$temp_c, na.rm = TRUE),
+      T_max = max(.data$temp_c, na.rm = TRUE),
+      mean_temp_c = mean(.data$temp_c, na.rm = TRUE),
+      R_n = sum(.data$Rn_MJ, na.rm = TRUE),
+      e_a = mean(.data$ea_kpa, na.rm = TRUE),
+      P_kpa = mean(.data$pressure / 1000, na.rm = TRUE),
+      wind = mean(.data$wind_speed, na.rm = TRUE),
+      precip_mm = sum(.data$precip_mm, na.rm = TRUE),
+      .groups = "drop"
+    )
+
+  df_daily |>
+    dplyr::rowwise() |>
+    dplyr::mutate(
+      ET0_mm = {
+        gamma = 0.000665 * .data$P_kpa
+        if (wind_height_m == 2) {
+          FAO56::ETo_FPM(
+            R_n = .data$R_n, G = 0, gamma = gamma, u_2 = .data$wind,
+            e_a = .data$e_a, T_min = .data$T_min, T_max = .data$T_max
+          )
+        } else {
+          FAO56::ETo_FPM(
+            R_n = .data$R_n, G = 0, gamma = gamma, u_z = .data$wind,
+            z = wind_height_m, e_a = .data$e_a, T_min = .data$T_min, T_max = .data$T_max
+          )
+        }
+      }
+    ) |>
+    dplyr::ungroup()
+}
+
+process_nc_file = function(nc_path, wind_height_m = config$wind_height_m) {
+  nc = tryCatch(
+    ncdf4::nc_open(nc_path),
+    error = function(e) {
+      warning("Could not open: ", nc_path)
+      NULL
+    }
+  )
+
+  if (is.null(nc)) return(NULL)
+  on.exit(ncdf4::nc_close(nc))
+
+  file_name = basename(nc_path)
+  site_hash = basename(dirname(nc_path))
+  parts = unlist(strsplit(gsub("\\.nc$", "", file_name), "\\."))
+
+  if (length(parts) < 3) {
+    warning("Unexpected filename: ", file_name)
+    return(NULL)
+  }
+
+  model_name = parts[1]
+  scenario_value = parts[2]
+  year_val = as.integer(parts[3])
+
+  time_steps = as.numeric(ncdf4::ncvar_get(nc, "time"))
+  temp_k = as.numeric(ncdf4::ncvar_get(nc, "air_temperature"))
+  wind_speed = as.numeric(ncdf4::ncvar_get(nc, "wind_speed"))
+  precip = as.numeric(ncdf4::ncvar_get(nc, "precipitation_flux"))
+  spec_hum = as.numeric(ncdf4::ncvar_get(nc, "specific_humidity"))
+  lw_rad = as.numeric(ncdf4::ncvar_get(nc, "surface_downwelling_longwave_flux_in_air"))
+  sw_rad = as.numeric(ncdf4::ncvar_get(nc, "surface_downwelling_shortwave_flux_in_air"))
+  pressure = as.numeric(ncdf4::ncvar_get(nc, "air_pressure"))
+  time_units = ncdf4::ncatt_get(nc, "time", "units")$value
+  datetime = parse_nc_time(time_steps, time_units, year_val)
+
+  df_hourly = tibble::tibble(
+    datetime = datetime, temp_c = temp_k - 273.15, wind_speed = wind_speed,
+    sw_rad = sw_rad, lw_rad = lw_rad, spec_hum = spec_hum,
+    pressure = pressure, precip_mm = precip * 3600
+  )
+
+  calculate_daily_et0(df_hourly, wind_height_m) |>
+    dplyr::mutate(
+      site_hash = .env$site_hash,
+      model = .env$model_name,
+      scenario = .env$scenario_value,
+      year = lubridate::year(.data$date),
+      day_of_year = lubridate::yday(.data$date)
+    ) |>
+    dplyr::select(
+      "date", "site_hash", "model", "scenario", "ET0_mm", "precip_mm",
+      "mean_temp_c", "T_min", "T_max", dplyr::everything()
+    )
+}
+
+# ---------- geo assignments ----------
+#convert the Cal-Adapt lat/lon locations into county assignments, attach the county to the climate records,
+#and create a lookup table for the design points
+
+site_dirs = dir_ls(base_dir, type = "directory")
+
+##site_index will store each site_hash and their lat/lon info
+site_index = map_dfr(site_dirs, function(dir_path) {
+  nc_files = dir_ls(dir_path, glob = "*.nc")
+  if (!length(nc_files)) return(NULL)
+
+  nc = nc_open(nc_files[1])
+  lat = as.numeric(ncvar_get(nc, "latitude"))
+  lon = as.numeric(ncvar_get(nc, "longitude"))
+  nc_close(nc)
+
+  tibble(site_hash = basename(dir_path), lat = lat, lon = lon)
+})
+
+##just for a check, should say 198
+cat(sprintf("Indexed %d unique Cal-Adapt design-point sites.\n", nrow(site_index)))
+
+##get CA county info
+ca_counties = counties(state = "CA", cb = TRUE, class = "sf") |> select(County = NAME)
+
+##sites_sf now holds each site_hash's lat/lon and point geometry
+sites_sf = st_as_sf(site_index, coords = c("lon", "lat"), crs = 4326, remove = FALSE)
+
+#double check both layers use same CRS
+ca_counties = st_transform(ca_counties, st_crs(sites_sf))
+
+##assign each site hash its county name
+sites_with_county = st_join(sites_sf, ca_counties, join = st_intersects, left = TRUE) |> st_drop_geometry()
+
+##can see the counties that are still not assigned to a site
+missing_counties = ca_counties |> filter(!County %in% sites_with_county$County)
+print(missing_counties$County)
+
+##assign whats leftover using nearest neighbor
+sites_projected = st_transform(sites_sf, 3310)
+missing_counties_projected = st_transform(missing_counties, 3310)
+missing_county_points = st_point_on_surface(missing_counties_projected)
+nearest_site_index = st_nearest_feature(missing_county_points, sites_projected)
+
+missing_county_lookup = tibble(County = missing_counties$County,
+  site_hash = sites_projected$site_hash[nearest_site_index])
+
+internal_site_lookup = sites_with_county |> filter(!is.na(County)) |> select(County, site_hash)
+county_site_lookup = bind_rows(internal_site_lookup, missing_county_lookup) |> distinct()
+
+##now have a lookup table with sites and their county info
+write.csv(county_site_lookup, file.path(output_dir, "caladapt_county_site_lookup.csv"), row.names = FALSE)
+
+#double check every county now has at least one climate site
+unmatched_counties = setdiff(ca_counties$County, county_site_lookup$County)
+print(unmatched_counties)
+if (length(unmatched_counties)) stop("Some counties still have no Cal-Adapt site assignment.")
+
+# ---------- process required NetCDF files ----------
+site_dirs = fs::dir_ls(base_dir, type = "directory")
+
+all_nc_files = purrr::map(site_dirs,
+  \(site_dir) {
+    file.path(site_dir, paste0(gcm_keep, ".", ssp_keep, ".", years_keep, ".nc"))
+  }
+) |>
+  unlist()
+
+missing_files = all_nc_files[!file.exists(all_nc_files)]
+
+if (length(missing_files)) {
+  stop(length(missing_files), " expected climate files are missing. First missing file: ",
+    missing_files[1])
+}
+
+cat("Processing", length(all_nc_files), "files across", length(site_dirs), "grid cells.\n")
+
+plan(multisession, workers = max(1, parallelly::availableCores() - 1))
+
+daily_climate_dataset = future_map_dfr(all_nc_files, process_nc_file,
+  wind_height_m = config$wind_height_m, .progress = TRUE, .options = furrr_options(seed = TRUE))
+
+# ---------- attach counties to daily climate records ----------
+daily_climate_county = daily_climate_dataset |> inner_join(county_site_lookup, by = "site_hash")
+
+county_daily_climate = daily_climate_county |>
+  group_by(County, date, model, scenario) |>
+  summarise(
+    ET0_mm = mean(ET0_mm, na.rm = TRUE),
+    precip_mm = mean(precip_mm, na.rm = TRUE),
+    mean_temp_c = mean(mean_temp_c, na.rm = TRUE),
+    min_temp_c = mean(T_min, na.rm = TRUE),
+    max_temp_c = mean(T_max, na.rm = TRUE),
+    n_climate_sites = n_distinct(site_hash),
+    .groups = "drop"
+  ) |>
+  mutate(Year = year(date), Week = isoweek(date)) |>
+  rename(GCM = model, SSP = scenario)
+
+write.csv(county_daily_climate, file.path(output_dir, "caladapt_county_daily_climate.csv"), row.names = FALSE)
+
+message("DONE: ", file.path(output_dir, "caladapt_county_daily_climate.csv"))