From c8b965adaf81f9beca5c1cd9dab5cbffd060947d Mon Sep 17 00:00:00 2001 From: Ananya Kulkarni Date: Mon, 6 Jul 2026 16:09:34 -0400 Subject: [PATCH 1/3] updates to scripts --- modules/data.remote/inst/predict_and_store.R | 1186 ++++++++++++++---- modules/data.remote/inst/scenarios.R | 1126 ++++++++--------- modules/data.remote/inst/transition_matrix.R | 195 +-- 3 files changed, 1576 insertions(+), 931 deletions(-) diff --git a/modules/data.remote/inst/predict_and_store.R b/modules/data.remote/inst/predict_and_store.R index 5bb6f6750b8..13853e5b6c9 100644 --- a/modules/data.remote/inst/predict_and_store.R +++ b/modules/data.remote/inst/predict_and_store.R @@ -1,299 +1,1013 @@ -##loop through all the design points, predicts up until 2050 (or another specified end_year), and stores them in a -#landIQ style dataset +#goes through all counties and predicts up until end_year using the optimized matrices & stores county-separated +#LandIQ-style datasets inside bau and nbs target subfolders setwd("/projectnb/dietzelab/ananyak") -library(data.table) -library(sf) -library(ggplot2) -library(tigris) -##-------setup------- -##read necessary files and check files are formatted properly -year_states = fread("year_states.csv") -tmat_df = fread("transition_matrix.csv") +library(data.table) -row_id_col = colnames(tmat_df)[1] -states = colnames(tmat_df)[-1] +##-------scenario setup------- +run_scenario = "BAU" #"BAU" or "NBS_Targets" -tmat_final = as.matrix(tmat_df[, -1, with = FALSE]) -rownames(tmat_final) = tmat_df[[row_id_col]] -colnames(tmat_final) = states -storage.mode(tmat_final) = "double" +prediction_root_dir = "county_landiq_predictions" +scenario_prediction_dir = file.path(prediction_root_dir, run_scenario) -stopifnot(all(rownames(tmat_final) == colnames(tmat_final))) -setDT(year_states) -tmat_year = tmat_final -start_info = year_states[, .SD[which.max(year)], by = parcel_id] +dir.create(scenario_prediction_dir, recursive = TRUE, showWarnings = FALSE) -#clean classes and set in same order as transition matrix -start_info[, dominant_crop := trimws(as.character(dominant_crop))] -states = rownames(tmat_final) +start_year = 2023L +end_year = 2045L -##specify how long into the future the user wants to predict to -end_year = 2050 +crop_matrix_dir = file.path("county_optimized_matrices", run_scenario) +crop_matrix_pattern = paste0("_crop_matrix_", run_scenario, "\\.csv$") +target_dir = file.path("county_optimized_targets", run_scenario) -##--------design point predictions--------- -design_points = fread('/projectnb/dietzelab/ccmmf/management/design_points_landiq_2018-2023.csv') +use_till_targets = dir.exists(target_dir) && + length(list.files(target_dir, pattern = paste0("crop_by_till_targets_.*_", run_scenario, "_", end_year, "\\.csv$") + )) > 0 -design_points[, parcel_id := as.character(parcel_id)] -start_info[, parcel_id := as.character(parcel_id)] -year_states[, parcel_id := as.character(parcel_id)] +##-------helper functions------- +safe_county_name = function(x) {gsub("[^A-Za-z0-9_]+", "_", x)} -#filter to only design point parcels -start_info = start_info[parcel_id %in% design_points$parcel_id] -year_states_design = year_states[parcel_id %in% design_points$parcel_id] +read_tmat = function(path) { + + tmat_df = fread(path) + + row_id_col = colnames(tmat_df)[1] + states = as.character(colnames(tmat_df)[-1]) + row_states = as.character(tmat_df[[row_id_col]]) + + tmat_final = as.matrix(tmat_df[, -1, with = FALSE]) + rownames(tmat_final) = row_states + colnames(tmat_final) = states + storage.mode(tmat_final) = "double" + + stopifnot(all(rownames(tmat_final) == colnames(tmat_final))) + + return(tmat_final)} -##function that pulls transition matrix probabilities -all_preds = start_info[, { +repair_transition_matrix = function(A, matrix_name = "matrix") { + + A[is.na(A)] = 0 + + #clamp tiny numerical optimizer noise + if (any(A < 0, na.rm = TRUE)) { + warning( + matrix_name, + " has negative probabilities. Minimum = ", + min(A, na.rm = TRUE), + ". Clamping negatives to 0.") + A[A < 0] = 0} - current_state = dominant_crop - years = seq(year + 1, end_year) + if (any(A > 1, na.rm = TRUE)) { + warning( + matrix_name, + " has probabilities > 1. Maximum = ", + max(A, na.rm = TRUE), + ". Clamping values above 1 to 1.") + A[A > 1] = 1} - preds = character(length(years)) - probs = numeric(length(years)) + row_sums = rowSums(A) + zero_rows = names(row_sums)[is.na(row_sums) | row_sums == 0] - for (i in seq_along(years)) { - p = tmat_year[current_state, ] + if (length(zero_rows) > 0) { + warning( + matrix_name, " has zero-sum rows. Setting those rows to self-loop: ", + paste(zero_rows, collapse = ", ")) - if (sum(p) == 0) { - preds[i] = NA_character_ - probs[i] = NA_real_ - next - } + for (s in zero_rows) { + A[s, ] = 0 + A[s, s] = 1}} + + row_sums = rowSums(A) + A = sweep(A, 1, row_sums, "/") + + return(A)} + +load_crop_matrices = function(crop_matrix_dir, crop_matrix_pattern, run_scenario) { + + matrix_files = list.files(crop_matrix_dir, pattern = crop_matrix_pattern, full.names = TRUE) + + if (length(matrix_files) == 0) { + stop("No optimized crop matrix files found in: ", crop_matrix_dir)} + + transition_mats = list() + + for (f in matrix_files) { - idx = which.max(p) - next_state = states[idx] + matrix_name = sub(paste0("_crop_matrix_", run_scenario, "\\.csv$"), "", basename(f)) - preds[i] = next_state - probs[i] = p[idx] + A = read_tmat(f) + A = repair_transition_matrix(A, paste0("optimized crop matrix ", matrix_name)) - current_state = next_state + transition_mats[[matrix_name]] = A} + + return(transition_mats)} + +load_crop_targets = function(target_dir, run_scenario, end_year) { + + pattern = paste0("crop_targets_.*_", run_scenario, "_", end_year, "\\.csv$") + files = list.files(target_dir, pattern = pattern, full.names = TRUE) + + if (length(files) == 0) { + stop("No crop target files found in: ", target_dir) } - .( - year = years, - pred_class = preds, - pred_prob = probs - ) - -}, by = parcel_id] - -#store the actual classes from 2018-2023 -actual_hist = year_states_design[, .(parcel_id, year, pred_class = NA_character_, - actual_class = dominant_crop)] - -#future predictions -preds_future = all_preds[, .( - parcel_id, - year, - pred_class, - actual_class = NA_character_)] - -##-----plotting------ -#combine for comparisons -plot_data = rbind(actual_hist, preds_future, fill = TRUE) - -sample_pids = sample(unique(plot_data$parcel_id), 1) - -plot_subset = plot_data[parcel_id %in% sample_pids] - -plot_subset[, pred_class := factor(pred_class, levels = states)] -plot_subset[, actual_class := factor(actual_class, levels = states)] - -ggplot(plot_subset, aes(x = year)) + - - geom_point( - data = plot_subset[!is.na(pred_class)], - aes(y = pred_class, color = pred_class), - size = 3 - ) + - - geom_point( - data = plot_subset[!is.na(actual_class)], - aes(y = actual_class), - shape = 1, - size = 3, - color = "black" - ) + - - facet_wrap(~parcel_id, ncol = 2) + - - scale_x_continuous( - limits = c(2018, end_year), - breaks = seq(2018, end_year, by = 2) - ) + - - ggtitle(sprintf("Predicted after Actual crop classes for parcel %s", sample_pids))+ - theme_minimal() - -##-------store predictions in landIQ style------- -##refers back to sublcass data and assigns the predictions a given sublcass based on fequencies -set.seed(123) - -#cleaning the original design point and prediction columns to be safe -design_points[, CLASS := as.character(CLASS)] -design_points[, SUBCLASS := as.character(SUBCLASS)] -preds_future[, parcel_id := as.character(parcel_id)] -preds_future[, pred_class := as.character(pred_class)] - -#getting the last observed landiq row for each parcel for current class/subclass before predictions start -last_obs = design_points[ - order(year, season), - .SD[.N], - by = parcel_id -][ - , .( - parcel_id, - last_CLASS = CLASS, - last_SUBCLASS = SUBCLASS - ) -] + out = rbindlist(lapply(files, fread), fill = TRUE) + + out[, county_safe := as.character(county_safe)] + out[, crop_state := as.character(crop_state)] + + if ("target_acres_used_for_opt" %in% names(out)) { + out[, target_acres := as.numeric(target_acres_used_for_opt)] + } else { + out[, target_acres := as.numeric(target_acres_raw)] + } + + out = out[ + !is.na(county_safe) & + !is.na(crop_state), + .(target_acres = sum(target_acres, na.rm = TRUE)), + by = .(county_safe, crop_state) + ] + + return(out[]) +} -#subclass probability table: if a class transition happens, draw subclass based on observed frequencies -subclass_probs = design_points[ - !is.na(CLASS) & !is.na(SUBCLASS), - .N, - by = .(CLASS, SUBCLASS) -] +make_matrix_powers = function(tmat, n_years) { + + states = rownames(tmat) + powers = vector("list", n_years) + + A_power = diag(length(states)) + rownames(A_power) = states + colnames(A_power) = states + + for (i in seq_len(n_years)) { + A_power = A_power %*% tmat + rownames(A_power) = states + colnames(A_power) = states + powers[[i]] = A_power + } + + return(powers) +} -subclass_probs[, prob := N / sum(N), by = CLASS] +predict_county_to_targets = function(start_info, tmat, targets_cty, + start_year, end_year, + state_col = "crop_class") { + + dt = copy(start_info) + states = rownames(tmat) + n_years = end_year - start_year + years = seq(start_year + 1L, end_year) + + dt = dt[get(state_col) %in% states] + dt[, start_CLASS := as.character(get(state_col))] + dt[, ACRES := as.numeric(ACRES)] + dt = dt[!is.na(ACRES) & ACRES > 0] + + if (nrow(dt) == 0) return(data.table()) + + powers = make_matrix_powers(tmat, n_years) + A_final = powers[[n_years]] + + # full target vector; states not in scenario target file become zero + target_vec = setNames(rep(0, length(states)), states) + if (nrow(targets_cty) > 0) { + matched = intersect(targets_cty$crop_state, states) + target_vec[matched] = targets_cty[match(matched, crop_state), target_acres] + } + + # force target total onto this county's modeled parcel-acre total + county_total = sum(dt$ACRES, na.rm = TRUE) + if (sum(target_vec, na.rm = TRUE) > 0) { + target_vec = target_vec / sum(target_vec, na.rm = TRUE) * county_total + } else { + # fallback: if no target file exists, preserve current acreage + cur = dt[, .(acres = sum(ACRES, na.rm = TRUE)), by = start_CLASS] + target_vec[cur$start_CLASS] = cur$acres + } + + # start by keeping all parcels in their observed class + dt[, final_CLASS := start_CLASS] + dt[, locked := FALSE] + + get_current = function(x) { + cur = x[, .(acres = sum(ACRES, na.rm = TRUE)), by = final_CLASS] + out = setNames(rep(0, length(states)), states) + out[cur$final_CLASS] = cur$acres + return(out) + } + + # greedily move parcels from surplus states into deficit states + for (to_state in states[order(-target_vec)]) { + + current_vec = get_current(dt) + need = target_vec[to_state] - current_vec[to_state] + + if (!is.finite(need) || need <= 0) next + + surplus_states = names(current_vec)[current_vec > target_vec + 1e-6] + surplus_states = setdiff(surplus_states, to_state) + + candidates = dt[ + !locked & + final_CLASS %in% surplus_states + ] + + if (nrow(candidates) == 0) next + + candidates[, prob_to := A_final[start_CLASS, to_state]] + candidates[, prob_current := A_final[start_CLASS, final_CLASS]] + candidates[, score := prob_to - prob_current] + + setorder(candidates, -score, -prob_to) + + candidates[, cum_acres := cumsum(ACRES)] + take_ids = candidates[cum_acres <= need | shift(cum_acres, fill = 0) < need, parcel_id] + + if (length(take_ids) > 0) { + dt[parcel_id %in% take_ids, `:=`( + final_CLASS = to_state, + locked = TRUE + )] + } + } + + # build annual time series: keep start class until conversion year + pred_list = vector("list", nrow(dt)) + + for (i in seq_len(nrow(dt))) { + + parcel = dt$parcel_id[i] + start_state = dt$start_CLASS[i] + final_state = dt$final_CLASS[i] + + pred_state = rep(start_state, n_years) + pred_prob = numeric(n_years) + + if (final_state == start_state) { + + for (yy in seq_len(n_years)) { + pred_prob[yy] = powers[[yy]][start_state, start_state] + } + + } else { + + conversion_idx = which(sapply(seq_len(n_years), function(k) { + powers[[k]][start_state, final_state] >= powers[[k]][start_state, start_state] + })) + + if (length(conversion_idx) == 0) { + conversion_idx = n_years + } else { + conversion_idx = conversion_idx[1] + } + + pred_state[conversion_idx:n_years] = final_state + + for (yy in seq_len(n_years)) { + pred_prob[yy] = powers[[yy]][start_state, pred_state[yy]] + } + } + + pred_list[[i]] = data.table( + parcel_id = parcel, + year = years, + CLASS = pred_state, + prob_crop_class = as.numeric(pred_prob) + ) + } + + out = rbindlist(pred_list, fill = TRUE) + return(out[]) +} -draw_subclass = function(class_name) { +predict_grouped_markov_to_targets = function(year_states, transition_mats, + crop_targets_2045, + group_col, + start_year, end_year, + state_col = "crop_class") { + + dt = copy(year_states) + all_preds = list() - choices = subclass_probs[CLASS == class_name] + groups = intersect(unique(na.omit(dt[[group_col]])), names(transition_mats)) - if (nrow(choices) == 0) { - return(NA_character_) + for (g in groups) { + + message("Predicting crop class with target allocation for county: ", g) + + dt_g = dt[get(group_col) == g] + tmat_g = transition_mats[[g]] + + start_info = dt_g[ + year <= start_year, + .SD[which.max(year)], + by = parcel_id + ] + + targets_g = crop_targets_2045[county_safe == g] + + preds_g = predict_county_to_targets( + start_info = start_info, + tmat = tmat_g, + targets_cty = targets_g, + start_year = start_year, + end_year = end_year, + state_col = state_col + ) + + if (nrow(preds_g) == 0) next + + preds_g[, (group_col) := g] + all_preds[[g]] = preds_g } - sample( - choices$SUBCLASS, - size = 1, - prob = choices$prob - ) + return(rbindlist(all_preds, fill = TRUE)) } -#adding previous observed class/subclass to predictions -preds_subclass = merge( - preds_future[, .(parcel_id, year, CLASS = pred_class)], - last_obs, - by = "parcel_id", - all.x = TRUE -) -setorder(preds_subclass, parcel_id, year) -preds_subclass[, SUBCLASS := NA_character_] +##--------load & clean all_data--------- +all_data = fread("all_data.csv") + +if ("V1" %in% names(all_data)) {all_data[, V1 := NULL]} + +if ("state" %in% names(all_data) && !"till_state" %in% names(all_data)) { + setnames(all_data, "state", "till_state")} + +if ("non_dom_prob" %in% names(all_data) && !"till_non_dom_prob" %in% names(all_data)) { + setnames(all_data, "non_dom_prob", "till_non_dom_prob")} + +setDT(all_data) + +required_cols = c("parcel_id", "year", "county", "crop_class", "ACRES") +missing_cols = setdiff(required_cols, names(all_data)) + +if (length(missing_cols) > 0) { + stop("all_data is missing required columns: ", paste(missing_cols, collapse = ", "))} + +all_data[, `:=`( + parcel_id = as.character(parcel_id), year = as.integer(year), county = as.character(county), + crop_class = trimws(as.character(crop_class)), ACRES = as.numeric(ACRES))] + +if (!"till_state" %in% names(all_data)) {all_data[, till_state := NA_character_]} else { + all_data[, till_state := trimws(as.character(till_state))]} + +if (!"county_geoid" %in% names(all_data)) {all_data[, county_geoid := NA_character_]} + +if (!"season" %in% names(all_data)) {all_data[, season := NA_integer_]} + +all_data[, county_safe := safe_county_name(county)] + +##--------load lookup table--------- + +lookup = fread("/projectnb/dietzelab/ccmmf/management/LandIQ_cropCode_lookup_table.csv") + +lookup[, `:=`( + CLASS = as.character(CLASS), SUBCLASS = as.character(SUBCLASS))] + +lookup_subclass = unique(lookup[, .( + CLASS, SUBCLASS, CLASS_desc, SUBCLASS_desc, PFT)], + by = c("CLASS", "SUBCLASS")) + +##--------load historical LandIQ source with subclass--------- + +with_subclass = fread("crops_full_counties.csv") + +if ("V1" %in% names(with_subclass)) {with_subclass[, V1 := NULL]} + +#standardize class, sublcass, season, and county name columns +if (!"CLASS" %in% names(with_subclass) && "crop_class" %in% names(with_subclass)) { + with_subclass[, CLASS := as.character(crop_class)]} + +if (!"SUBCLASS" %in% names(with_subclass) && "subclass" %in% names(with_subclass)) { + setnames(with_subclass, "subclass", "SUBCLASS")} + +if (!"county_safe" %in% names(with_subclass) && "county" %in% names(with_subclass)) { + with_subclass[, county_safe := safe_county_name(county)]} + +if (!"season" %in% names(with_subclass)) {with_subclass[, season := 0L]} + +required_subclass_cols = c("parcel_id", "year", "county_safe", "CLASS", "SUBCLASS") +missing_subclass_cols = setdiff(required_subclass_cols, names(with_subclass)) + +if (length(missing_subclass_cols) > 0) { + stop("crops_full_counties.csv is missing required columns: ", + paste(missing_subclass_cols, collapse = ", "), "\nAvailable columns: ", + paste(names(with_subclass), collapse = ", "))} + +with_subclass[, `:=`( + parcel_id = as.character(parcel_id), year = as.integer(year), + county_safe = as.character(county_safe), CLASS = trimws(as.character(CLASS)), + SUBCLASS = trimws(as.character(SUBCLASS)), season = as.integer(season))] -#assigning subclass year by year -for (p in unique(preds_subclass$parcel_id)) { +with_subclass[CLASS == "" | CLASS == "NA", CLASS := NA_character_] +with_subclass[SUBCLASS == "" | SUBCLASS == "NA", SUBCLASS := NA_character_] +with_subclass[is.na(season), season := 0L] + +message("Loaded subclass source with ", nrow(with_subclass[!is.na(SUBCLASS)]), + " rows with non-NA SUBCLASS.") + +if (!"CLASS" %in% names(all_data)) {all_data[, CLASS := as.character(crop_class)]} + +assign_predicted_subclass = function(future_landiq, subclass_obs, lookup_subclass, + crop_col = "CLASS", group_col = "county_safe", + start_year = 2023L) { + + dt = copy(future_landiq) + obs = copy(subclass_obs) + + dt[, orig_order := .I] + + obs[, `:=`( + parcel_id = as.character(parcel_id), year = as.integer(year), + county_safe = as.character(county_safe), CLASS = as.character(CLASS), + SUBCLASS = as.character(SUBCLASS))] + + if (!"season" %in% names(obs)) { + obs[, season := 0L] + } else { + obs[, season := as.integer(season)] + obs[is.na(season), season := 0L] + } + + #only use historical data up to the prediction start year + obs = obs[year <= start_year] + + dt[, parcel_id := as.character(parcel_id)] + dt[, CLASS := as.character(get(crop_col))] + dt[, (group_col) := as.character(get(group_col))] + + #remove old lookup/subclass fields if function is accidentally rerun + old_lookup_cols = intersect(c("SUBCLASS", "CLASS_desc", "SUBCLASS_desc", "PFT"), + names(dt)) + + if (length(old_lookup_cols) > 0) {dt[, (old_lookup_cols) := NULL]} + + #global subclass probabilities + global_probs = obs[ + !is.na(CLASS) & !is.na(SUBCLASS), + .N, by = .(CLASS, SUBCLASS)] + + if (nrow(global_probs) > 0) {global_probs[, prob := N / sum(N), by = CLASS]} + + #county-specific subclass probabilities + group_probs = obs[ + !is.na(CLASS) & !is.na(SUBCLASS), + .N, by = .(county_safe, CLASS, SUBCLASS)] + + if (nrow(group_probs) > 0) {group_probs[, prob := N / sum(N), by = .(county_safe, CLASS)]} + + #lookup fallback if a predicted class has no observed subclass distribution + lookup_probs = unique(lookup_subclass[ + !is.na(CLASS) & !is.na(SUBCLASS), + .(CLASS, SUBCLASS)]) + + if (nrow(lookup_probs) > 0) {lookup_probs[, N := 1L] + lookup_probs[, prob := 1 / .N, by = CLASS]} + + # Last observed class/subclass per parcel + county from subclass source + last_obs_source = obs[ + !is.na(CLASS) & !is.na(SUBCLASS)] + + last_obs = last_obs_source[ + order(year, season), + .SD[.N], by = .(parcel_id, county_safe)] + + if (!"SUBCLASS" %in% names(last_obs)) {last_obs[, SUBCLASS := NA_character_]} + + if (!"CLASS" %in% names(last_obs)) {last_obs[, CLASS := NA_character_]} + + last_obs = last_obs[ + , + .( + parcel_id, county_safe, last_CLASS = CLASS, last_SUBCLASS = SUBCLASS)] + + dt = merge(dt, last_obs, by = c("parcel_id", "county_safe"), all.x = TRUE) + + setorder(dt, parcel_id, year) - idxs = which(preds_subclass$parcel_id == p) + #find runs of same predicted CLASS within each parcel + dt[, prev_CLASS := shift(CLASS), by = .(parcel_id, county_safe)] + dt[is.na(prev_CLASS), prev_CLASS := last_CLASS] - prev_class = preds_subclass$last_CLASS[idxs[1]] - prev_subclass = preds_subclass$last_SUBCLASS[idxs[1]] + dt[, new_run := fifelse(is.na(CLASS), FALSE, is.na(prev_CLASS) | CLASS != prev_CLASS)] - for (i in idxs) { + dt[, run_id := cumsum(new_run), by = .(parcel_id, county_safe)] + + dt[, SUBCLASS := NA_character_] + + #if predicted class matches last observed class, carry forward real historical subclass + dt[run_id == 0 & + !is.na(CLASS) & + !is.na(last_CLASS) & + CLASS == last_CLASS & + !is.na(last_SUBCLASS), + SUBCLASS := last_SUBCLASS] + + #for remaining class-runs, draw subclass from county/class distribution + run_table = unique(dt[ + !is.na(CLASS) & is.na(SUBCLASS), + .(parcel_id, county_safe, run_id, CLASS)]) + + if (nrow(run_table) > 0) { + + run_table[, drawn_SUBCLASS := NA_character_] - current_class = preds_subclass$CLASS[i] + draw_groups = unique(run_table[, .(county_safe, CLASS)]) - if (is.na(current_class)) { - preds_subclass$SUBCLASS[i] = NA_character_ + for (ii in seq_len(nrow(draw_groups))) { + + cty = draw_groups$county_safe[ii] + cls = draw_groups$CLASS[ii] + + idx = which(run_table$county_safe == cty & run_table$CLASS == cls) + + choices = group_probs[county_safe == cty & CLASS == cls] + + if (nrow(choices) == 0) { + choices = global_probs[CLASS == cls] + } - } else if (!is.na(prev_class) && current_class == prev_class) { + if (nrow(choices) == 0) { + choices = lookup_probs[CLASS == cls] + } - #if a crop class does not transition, keep the subclass its currently in - preds_subclass$SUBCLASS[i] = prev_subclass + if (!("prob" %in% names(choices))) { + choices[, prob := NA_real_] + } + choices = choices[ + !is.na(SUBCLASS) & + !is.na(prob) & + is.finite(prob) & + prob > 0 + ] + + if (nrow(choices) > 0) { + choices[, prob := prob / sum(prob)] + + run_table$drawn_SUBCLASS[idx] = sample( + choices$SUBCLASS, + size = length(idx), + replace = TRUE, + prob = choices$prob + ) + } + } + + dt = merge(dt, run_table[, .(parcel_id, county_safe, run_id, drawn_SUBCLASS)], + by = c("parcel_id", "county_safe", "run_id"), all.x = TRUE) + + dt[is.na(SUBCLASS), SUBCLASS := drawn_SUBCLASS] + dt[, drawn_SUBCLASS := NULL]} + + helper_cols = intersect(c("last_CLASS", "last_SUBCLASS", "prev_CLASS", "new_run", "run_id"), + names(dt)) + + dt[, (helper_cols) := NULL] + + if (!"SUBCLASS" %in% names(dt)) {dt[, SUBCLASS := NA_character_]} + + if (!"CLASS" %in% names(dt)) {dt[, CLASS := NA_character_]} + + #add lookup descriptions/PFT + dt = merge(dt, lookup_subclass, by = c("CLASS", "SUBCLASS"), all.x = TRUE) + + setorder(dt, orig_order) + dt[, orig_order := NULL] + + return(dt[])} + +##--------crop prediction functions--------- +predict_markov = function(start_info, tmat, end_year, id_col = "parcel_id", time_col = "year", state_col = "crop_class") { + + dt = copy(start_info) + + if (!(id_col %in% names(dt))) {stop("id_col not found in start_info: ", id_col)} + if (!(time_col %in% names(dt))) {stop("time_col not found in start_info: ", time_col)} + if (!(state_col %in% names(dt))) {stop("state_col not found in start_info: ", state_col)} + + if (id_col != "parcel_id") {setnames(dt, id_col, "parcel_id")} + if (time_col != "year") {setnames(dt, time_col, "year")} + + states = rownames(tmat) + + dt = dt[get(state_col) %in% states] + + if (nrow(dt) == 0) {return(data.table())} + + #clean matrix + tmat[is.na(tmat)] = 0 + tmat[tmat < 0] = 0 + + row_sums = rowSums(tmat) + for (s in states[row_sums <= 0 | is.na(row_sums)]) { + tmat[s, ] = 0 + tmat[s, s] = 1} + + tmat = sweep(tmat, 1, rowSums(tmat), "/") + + preds = dt[, { + + start_state = get(state_col) + years = seq(year + 1L, end_year) + n_years = length(years) + + if (is.na(start_state) || !(start_state %in% states) || n_years == 0) { + .(year = integer(), CLASS = character(), prob_crop_class = numeric()) } else { - #if crop class transitions, draw subclass from new class distribution - preds_subclass$SUBCLASS[i] = draw_subclass(current_class) + # Track cumulative transition probabilities from start year to each future year + A_power = diag(length(states)) + rownames(A_power) = states + colnames(A_power) = states + + yearly_probs = vector("list", n_years) + + for (i in seq_len(n_years)) { + A_power = A_power %*% tmat + rownames(A_power) = states + colnames(A_power) = states + + p_i = as.numeric(A_power[start_state, ]) + names(p_i) = states + p_i[is.na(p_i)] = 0 + + if (sum(p_i) > 0) { + p_i = p_i / sum(p_i) + } + + yearly_probs[[i]] = p_i + } + + final_probs = yearly_probs[[n_years]] + final_state = names(final_probs)[which.max(final_probs)] + final_prob = final_probs[final_state] + + # Default: parcel stays constant + pred_state = rep(start_state, n_years) + pred_prob = sapply(yearly_probs, function(p) p[start_state]) + + # If final most-likely state is different, allow ONE conversion only + if (!is.na(final_state) && final_state != start_state) { + + # convert only when final state becomes at least as likely as start state + conversion_idx = which(sapply(yearly_probs, function(p) { + p[final_state] >= p[start_state] + })) + + if (length(conversion_idx) == 0) { + conversion_idx = n_years + } else { + conversion_idx = conversion_idx[1] + } + + pred_state[conversion_idx:n_years] = final_state + pred_prob[conversion_idx:n_years] = sapply( + yearly_probs[conversion_idx:n_years], + function(p) p[final_state] + ) + } + + .( + year = years, + CLASS = pred_state, + prob_crop_class = as.numeric(pred_prob) + ) } - #update previous class/subclass for next predicted year - prev_class = current_class - prev_subclass = preds_subclass$SUBCLASS[i] - } + }, by = parcel_id] + + return(preds) } -#adding class/subclass descriptions from lookup table -lookup = fread('/projectnb/dietzelab/ccmmf/management/LandIQ_cropCode_lookup_table.csv') +predict_grouped_markov = function(year_states, transition_mats, group_col, start_year, end_year, id_col = "parcel_id", + time_col = "year", state_col = "crop_class") { + + dt = copy(year_states) + + if (!(group_col %in% names(dt))) {stop("group_col not found in year_states: ", group_col)} + + if (id_col != "parcel_id") {setnames(dt, id_col, "parcel_id")} + + if (time_col != "year") {setnames(dt, time_col, "year")} + + all_preds = list() + + groups = intersect(unique(na.omit(dt[[group_col]])), names(transition_mats)) + + if (length(groups) == 0) {stop("No overlapping groups between all_data and transition matrices.")} + + for (g in groups) { + + message("Predicting crop class for county: ", g) + + dt_g = dt[get(group_col) == g] + tmat_g = transition_mats[[g]] + + # latest observed state for each parcel up to start_year + start_info = dt_g[year <= start_year, + .SD[which.max(year)], by = parcel_id] + + preds_g = predict_markov(start_info = start_info, tmat = tmat_g, end_year = end_year, state_col = state_col) + + if (nrow(preds_g) == 0) { + next} + + preds_g[, (group_col) := g] + + all_preds[[g]] = preds_g} + + return(rbindlist(all_preds, fill = TRUE))} + +##--------till target functions--------- -lookup[, CLASS := as.character(CLASS)] -lookup[, SUBCLASS := as.character(SUBCLASS)] +load_crop_by_till_targets = function(target_dir, run_scenario, end_year) { + + pattern = paste0("crop_by_till_targets_.*_", run_scenario, "_", end_year, "\\.csv$") + + files = list.files(target_dir, pattern = pattern, full.names = TRUE) + + if (length(files) == 0) {stop("No crop-by-till target files found in: ", target_dir)} + + out = rbindlist(lapply(files, fread), fill = TRUE) + + out[, county_safe := as.character(county_safe)] + out[, crop_state := as.character(crop_state)] + out[, till_state := as.character(till_state)] + out[, target_acres_raw := as.numeric(target_acres_raw)] + + return(out)} -lookup_subclass = unique(lookup[, .( - CLASS, - SUBCLASS, - CLASS_desc, - SUBCLASS_desc, - PFT -)]) - -preds_subclass = merge( - preds_subclass, - lookup_subclass, - by = c("CLASS", "SUBCLASS"), - all.x = TRUE +build_annual_till_targets = function(all_data, till_targets_2045, start_year, end_year) { + + latest_obs = all_data[year <= start_year & !is.na(crop_class) & !is.na(till_state), + .SD[which.max(year)], by = parcel_id] + + baseline_till = latest_obs[ + , + .(baseline_acres = sum(ACRES, na.rm = TRUE)), by = .(county_safe, crop_state = crop_class, till_state)] + + if (nrow(baseline_till) > 0) {baseline_till[, baseline_share := baseline_acres / sum(baseline_acres), + by = .(county_safe, crop_state)] + } else { + baseline_till[, baseline_share := numeric()]} + + target_till = copy(till_targets_2045) + + target_till[, target_share := target_acres_raw / sum(target_acres_raw), by = .(county_safe, crop_state)] + + all_counties = unique(target_till$county_safe) + all_crops = unique(target_till$crop_state) + all_till_states = unique(c(baseline_till$till_state, target_till$till_state)) + + annual_till_targets = CJ(county_safe = all_counties, crop_state = all_crops, till_state = all_till_states, + year = seq(start_year + 1L, end_year)) + + annual_till_targets = merge(annual_till_targets, baseline_till[, .(county_safe, crop_state, till_state, baseline_share)], + by = c("county_safe", "crop_state", "till_state"), all.x = TRUE) + + annual_till_targets = merge(annual_till_targets, target_till[, .(county_safe, crop_state, till_state, target_share)], + by = c("county_safe", "crop_state", "till_state"), all.x = TRUE) + + annual_till_targets[is.na(baseline_share), baseline_share := 0] + annual_till_targets[is.na(target_share), target_share := 0] + + annual_till_targets[, ramp := (year - start_year) / (end_year - start_year)] + + annual_till_targets[, till_share := + (1 - ramp) * baseline_share + ramp * target_share] + + annual_till_targets[ + , + share_sum := sum(till_share, na.rm = TRUE), by = .(county_safe, crop_state, year)] + + annual_till_targets[share_sum > 0, till_share := till_share / share_sum] + + annual_till_targets[share_sum <= 0, till_share := 1 / .N, by = .(county_safe, crop_state, year)] + + annual_till_targets[, share_sum := NULL] + + return(annual_till_targets)} + +assign_till_by_targets = function(future_landiq, annual_till_targets) { + + dt = copy(future_landiq) + targets_dt = copy(annual_till_targets) + + dt[, CLASS := as.character(CLASS)] + dt[, ACRES := as.numeric(ACRES)] + dt[is.na(ACRES) | !is.finite(ACRES) | ACRES < 0, ACRES := 0] + + targets_dt[, county_safe := as.character(county_safe)] + targets_dt[, crop_state := as.character(crop_state)] + targets_dt[, till_state := as.character(till_state)] + targets_dt[, year := as.integer(year)] + targets_dt[, till_share := as.numeric(till_share)] + + targets_dt = targets_dt[ + !is.na(till_state) & + !is.na(till_share) & + is.finite(till_share) & + till_share > 0] + + setkey(targets_dt, county_safe, year, crop_state) + + dt[, row_id := .I] + + out = dt[, { + + cty = county_safe[1] + yr = year[1] + cls = CLASS[1] + + targets = targets_dt[.(cty, yr, cls), nomatch = 0] + temp = copy(.SD) + + if (nrow(targets) == 0 || sum(targets$till_share, na.rm = TRUE) <= 0) { + + temp[, till_state := NA_character_] + temp[, prob_till_state := NA_real_] + temp + + } else { + + targets = copy(targets) + targets[, till_share := till_share / sum(till_share)] + setorder(targets, till_state) + + group_total_acres = sum(temp$ACRES, na.rm = TRUE) + + if (!is.finite(group_total_acres) || group_total_acres <= 0) { + + assigned = sample(targets$till_state, size = nrow(temp), replace = TRUE, + prob = targets$till_share) + + prob_lookup = setNames(targets$till_share, targets$till_state) + + temp[, till_state := assigned] + temp[, prob_till_state := prob_lookup[till_state]] + temp + + } else { + + #randomize row order so the same parcels are not always assigned first + temp[, rand_order := runif(.N)] + setorder(temp, rand_order) + + targets[, target_acres := till_share * group_total_acres] + targets[, upper_acres := cumsum(target_acres)] + + #use cumulative parcel acreage midpoint for assignment + temp[, cum_mid_acres := cumsum(ACRES) - ACRES / 2] + + idx = findInterval(temp$cum_mid_acres, targets$upper_acres) + 1L + idx[idx < 1L] = 1L + idx[idx > nrow(targets)] = nrow(targets) + + temp[, till_state := targets$till_state[idx]] + + prob_lookup = setNames(targets$till_share, targets$till_state) + temp[, prob_till_state := prob_lookup[till_state]] + + temp[, c("rand_order", "cum_mid_acres") := NULL] + temp + } + } + + }, by = .(county_safe, year, CLASS)] + + setorder(out, row_id) + out[, row_id := NULL] + + return(out[])} +##--------run predictions--------- +crop_mats = load_crop_matrices(crop_matrix_dir = crop_matrix_dir, crop_matrix_pattern = crop_matrix_pattern, + run_scenario = run_scenario) + +missing_mats = setdiff(unique(all_data$county_safe), names(crop_mats)) + +if (length(missing_mats) > 0) { + message("Counties in all_data without crop matrices: ", length(missing_mats))} + +crop_targets_2045 = load_crop_targets( + target_dir = target_dir, + run_scenario = run_scenario, + end_year = end_year ) -#adding parcel information from original design points -parcel_meta = design_points[ - order(year, season), - .SD[.N], - by = parcel_id +future_crop = predict_grouped_markov_to_targets( + year_states = all_data, + transition_mats = crop_mats, + crop_targets_2045 = crop_targets_2045, + group_col = "county_safe", + start_year = start_year, + end_year = end_year, + state_col = "crop_class" +) + +if (nrow(future_crop) == 0) {stop("No future crop predictions were generated.")} + +##--------add parcel metadata--------- +parcel_meta = all_data[year <= start_year, + .SD[which.max(year)], by = parcel_id ][ - , .(parcel_id, site_id, lon, lat) + , + .( + parcel_id, county, county_geoid, county_safe, ACRES, till_state_2023 = till_state)] + +future_landiq = merge(future_crop, parcel_meta, by = c("parcel_id", "county_safe"), all.x = TRUE) + +##--------assign till states--------- +if (use_till_targets) { + + message("Assigning till states using ", run_scenario, " crop-by-till targets.") + + till_targets_2045 = load_crop_by_till_targets(target_dir = target_dir, run_scenario = run_scenario, end_year = end_year) + + annual_till_targets = build_annual_till_targets(all_data = all_data, till_targets_2045 = till_targets_2045, + start_year = start_year, end_year = end_year) + + future_landiq = assign_till_by_targets(future_landiq = future_landiq, annual_till_targets = annual_till_targets) + +} else { + + message(run_scenario, " has no crop-by-till target files. Carrying forward 2023 observed till_state.") + + future_landiq[, till_state := till_state_2023] + future_landiq[, prob_till_state := 1]} + +##--------final formatting--------- +future_landiq[, season := NA_integer_] +future_landiq[, source := "predicted"] +future_landiq[, scenario := run_scenario] + +future_landiq = future_landiq[year >= start_year + 1L & year <= end_year] + +#add subclass + lookup fields +future_landiq = assign_predicted_subclass( future_landiq = future_landiq, + subclass_obs = with_subclass, + lookup_subclass = lookup_subclass, + crop_col = "CLASS", group_col = "county_safe", + start_year = start_year) + +#clean optional helper columns +if ("till_state_2023" %in% names(future_landiq)) {future_landiq[, till_state_2023 := NULL]} + +#make sure optional columns exist before selecting +optional_cols = c("SUBCLASS", "CLASS_desc", "SUBCLASS_desc", "PFT") + +for (cc in optional_cols) {if (!(cc %in% names(future_landiq))) { + future_landiq[, (cc) := NA_character_]}} + +final_cols = c("parcel_id", "county", "county_geoid", "county_safe", "year", "season", "CLASS", "SUBCLASS", "CLASS_desc", + "SUBCLASS_desc", "PFT", "till_state", "prob_crop_class", "prob_till_state", "ACRES", "source", "scenario") + +future_landiq = future_landiq[, ..final_cols] + +setorder(future_landiq, county_safe, parcel_id, year) + + +pred_2045_check = future_landiq[ + year == end_year, + .(predicted_acres = sum(ACRES, na.rm = TRUE)), + by = .(county_safe, CLASS) ] -preds_landiq = merge( - preds_subclass, - parcel_meta, - by = "parcel_id", - all.x = TRUE -) +target_2045_check = copy(crop_targets_2045) +setnames(target_2045_check, "crop_state", "CLASS") -#leaving season blank for future years -preds_landiq[, season := NA_integer_] - -#keep the columns that are in the design points file -preds_landiq = preds_landiq[, .( - site_id, - parcel_id, - lon, - lat, - year, - season, - CLASS, - SUBCLASS, - CLASS_desc, - SUBCLASS_desc, - PFT -)] - -#saving observed and predicted crop classes for the design points as separate files - -design_points_observed = copy(design_points[year >= 2018 & year <= 2023]) -design_points_predicted = copy(preds_landiq[year >= 2024 & year <= end_year]) - -design_points_observed[, source := "observed"] -design_points_predicted[, source := "predicted"] - -setorder(design_points_observed, parcel_id, year, season) -setorder(design_points_predicted, parcel_id, year, season) - -fwrite( - design_points_observed, - "design_points_observed.csv" +check_2045 = merge( + pred_2045_check, + target_2045_check, + by = c("county_safe", "CLASS"), + all = TRUE ) -fwrite( - design_points_predicted, - sprintf("predicted_2024_%s.csv", end_year) +check_2045[is.na(predicted_acres), predicted_acres := 0] +check_2045[is.na(target_acres), target_acres := 0] +check_2045[, diff_acres := predicted_acres - target_acres] +check_2045[, abs_diff_acres := abs(diff_acres)] + +check_path = file.path( + scenario_prediction_dir, + paste0("prediction_target_check_", run_scenario, "_", end_year, ".csv") ) +fwrite(check_2045, check_path) +message("Wrote 2045 target check to: ", check_path) + +##--------write county-separated outputs--------- +for (cty in unique(na.omit(future_landiq$county_safe))) { + + out_path = file.path(scenario_prediction_dir, paste0(cty, "_predicted_2024_", end_year, ".csv")) + + fwrite(future_landiq[county_safe == cty], out_path)} + +message("Finished writing predictions to: ", scenario_prediction_dir) \ No newline at end of file diff --git a/modules/data.remote/inst/scenarios.R b/modules/data.remote/inst/scenarios.R index 4cf8d137593..47b48513bf6 100644 --- a/modules/data.remote/inst/scenarios.R +++ b/modules/data.remote/inst/scenarios.R @@ -1,592 +1,596 @@ -#start from original transition matrix (tmat_final) and build initial state vector for all classes to find the new optimized -#transition matrix based on a user specified goal and given constraints. +setwd("/projectnb/dietzelab/ananyak") library(data.table) +library(readxl) library(nloptr) -library(expm) -library(dplyr) -library(ggplot2) -library(reshape2) -library(networkD3) - -setwd("/projectnb/dietzelab/ananyak") - -##------------setup----------------- -#load transition matrix and reformat - #keep row labels/states before dropping first column - #drop V1/class column to keep numeric matrix - #check matrix is square and rows sum to 1 -tmat_df = fread("transition_matrix.csv") - -row_id_col = colnames(tmat_df)[1] -states = colnames(tmat_df)[-1] - -A_orig = as.matrix(tmat_df[, ..states]) -rownames(A_orig) = tmat_df[[row_id_col]] -colnames(A_orig) = states -storage.mode(A_orig) = "double" - -n = nrow(A_orig) - -stopifnot(nrow(A_orig) == ncol(A_orig)) -stopifnot(all(rownames(A_orig) == colnames(A_orig))) -stopifnot(all(abs(rowSums(A_orig) - 1) < 1e-8)) - -#have to reread original data and add back acres column for initial state vector -path_management = "/projectnb/dietzelab/ccmmf/management" -path_landiq_v4 = "/projectnb/dietzelab/ccmmf/LandIQ-harmonized-v4.1" - -lookup = fread(file.path(path_management, "LandIQ_cropCode_lookup_table.csv")) -ag_classes = unique(lookup[is_agricultural == TRUE, as.character(CLASS)]) - -year_min = 2018L -year_max = 2023L - -crops_full = as.data.table( - arrow::open_dataset(file.path(path_landiq_v4, "crops_all_years.parq")) |> - filter(year >= year_min, year <= year_max, CLASS %in% ag_classes) |> - select(parcel_id, year, season, CLASS, SUBCLASS, ACRES) |> - collect() -) - -crops_full[, `:=`( - parcel_id = as.character(parcel_id), - year = as.integer(year), - season = as.integer(season), - CLASS = as.character(CLASS), - SUBCLASS = as.character(SUBCLASS), - ACRES = as.integer(ACRES) -)] - -setDT(crops_full) - -##filter to 2023 only (last observed year = prediction starting year) -crops_full_2023 = crops_full[year == 2023] - -#sum total acres by class, then convert to % -total_land = sum(crops_full_2023$ACRES, na.rm = TRUE) -land_by_class = aggregate(ACRES ~ CLASS, data = crops_full_2023, FUN = sum, na.rm = TRUE) -land_by_class$class_land_percs = land_by_class$ACRES / total_land - -#X0 has to be in the same order as the transition matrix states -X0_named = setNames(land_by_class$class_land_percs, land_by_class$CLASS) - -#reorder to exactly match transition matrix states -X0 = X0_named[states] - -#classes in tmat but missing in 2023 get 0 -X0[is.na(X0)] = 0 - -stopifnot(all(names(X0) == states)) -stopifnot(abs(sum(X0) - 1) < 1e-8) - -##drop v1 column from tmat for exact ordering -tmat_df = tmat_df %>% select(-V1) - -##-----------scenario goal inputs----------- - -#example: x% increase in crop class __ after n years/steps -target_crop = "V" - -target_val = X0[target_crop] * 1.30 -steps = 15 - -##---------objective and constraints------------- - #1.Minimize the "distance" from the original matrix - #2.x is a vector of the matrix elements, length n^2 - -obj_fun = function(x) { - A_new = matrix(x, nrow = n, byrow = TRUE) - sum((A_new - A_orig)^2) -} - -constr_fun = function(x) { - A_new = matrix(x, nrow = n, byrow = TRUE) +library(expm) + +##-------setup------- +scenario_name = "NBS_Targets" #bau or nbs targets +start_year = 2023L +end_year = 2045L +steps = end_year - start_year + +lambda_target = 1e6 +maxeval_optimizer = 50000 +maxtime_optimizer = 300 +scale_crop_targets_to_x0 = TRUE +nominal_zero_acres = 0.01 + +run_all_counties = TRUE #change to false if focusing on specific county/counties + +##to test specific counties if needed +#counties_manual = c("Fresno") + +scenario_out_dir = file.path("county_optimized_matrices", scenario_name) +target_out_dir = file.path("county_optimized_targets", scenario_name) + +dir.create(scenario_out_dir, recursive = TRUE, showWarnings = FALSE) +dir.create(target_out_dir, recursive = TRUE, showWarnings = FALSE) + +##-------helper functions------- +safe_county_name = function(x) {gsub("[^A-Za-z0-9_]+", "_", x)} + +normalize_crop_key = function(x) { + x = tolower(trimws(as.character(x))) + x = gsub("&", "and", x) + x = gsub("[[:punct:]]+", " ", x) + x = gsub("\\s+", " ", x) + trimws(x)} + +check_required_cols = function(dt, required_cols, dt_name) { + missing_cols = setdiff(required_cols, names(dt)) + if (length(missing_cols) > 0) { + stop(dt_name, " is missing required columns: ", paste(missing_cols, collapse = ", "), + "\nAvailable columns: ", paste(names(dt), collapse = ", "))}} + +read_tmat = function(path) { + tmat_df = fread(path) + row_id_col = colnames(tmat_df)[1] + states = as.character(colnames(tmat_df)[-1]) + row_states = as.character(tmat_df[[row_id_col]]) + tmat_final = as.matrix(tmat_df[, -1, with = FALSE]) + rownames(tmat_final) = row_states + colnames(tmat_final) = states + storage.mode(tmat_final) = "double" + stopifnot(all(rownames(tmat_final) == colnames(tmat_final))) + return(tmat_final)} + +repair_transition_matrix = function(A, matrix_name = "matrix") { + A[is.na(A)] = 0 - X_end = X0 %*% (A_new %^% steps) - colnames(X_end) = states + row_sums = rowSums(A) + zero_rows = names(row_sums)[is.na(row_sums) | row_sums == 0] - target_const = X_end[1, target_crop] - target_val - row_sums_const = rowSums(A_new) - 1 - - c(target_const, row_sums_const) -} - -##-----------run optimizer------------------ -#starting point = original matrix flattened into a vector -init_x = as.vector(t(A_orig)) - -#COBYLA for non-linear constraints without needing derivatives -res = nloptr(x0 = init_x, - eval_f = obj_fun, - eval_g_eq = constr_fun, - lb = rep(0, n^2), # Probabilities can't be negative - ub = rep(1, n^2), # Probabilities can't exceed 1 - opts = list( - algorithm = "NLOPT_LN_COBYLA", - xtol_rel = 1e-5, - ##how many max iterations - maxeval = 25000, - print_level = 1 - )) - -##-------outputs-------- - -A_final = matrix(res$solution, nrow = n, byrow = TRUE) -rownames(A_final) = states -colnames(A_final) = states - -X_end_orig = X0 %*% (A_orig %^% steps) -X_end_final = X0 %*% (A_final %^% steps) - -colnames(X_end_orig) = states -colnames(X_end_final) = states - -print("Optimizer status:") -print(res$status) -print(res$message) - -print("Optimized Matrix A:") -print(round(A_final, 4)) - -target_crop -X0[target_crop] -target_val -X_end_orig[1, target_crop] -X_end_final[1, target_crop] - -round(rowSums(A_final), 8) -max(abs(A_final - A_orig)) - -constr_fun(res$solution) - -##--------------visualizations------------------ - -## time series distributions -get_dist_over_time = function(X0, A, steps, states, scenario_name) { - out = lapply(0:steps, function(t) { - if (t == 0) { - Xt = X0 - } else { - Xt = as.numeric(X0 %*% (A %^% t)) - } + if (length(zero_rows) > 0) { + warning(matrix_name, " has zero-sum rows. Setting those rows to self-loop: ", + paste(zero_rows, collapse = ", ")) - data.table( - step = t, - year = 2023 + t, - CLASS = states, - prop_land = as.numeric(Xt), - scenario = scenario_name - ) - }) - - rbindlist(out) -} - -dist_orig = get_dist_over_time(X0, A_orig, steps, states, "Original transition matrix") -dist_final = get_dist_over_time(X0, A_final, steps, states, "Optimized transition matrix") - -dist_all = rbind(dist_orig, dist_final) - -#1.stacked area chart of projected land distribution by crop class (from optimized matrix) -##*changes are pretty small, so this graph is better for showing the final overall distribution, not comparisons -ggplot(dist_final, - aes(x = year, y = prop_land, fill = CLASS)) + - geom_area() + - labs( - title = "Optimized projected land distribution by crop class", - x = "Year", - y = "Fraction of total land", - fill = "Class" - ) + - theme_minimal() - -#2.differences in original vs optimized land share -final_compare = merge( - dist_orig[step == steps, .(CLASS, orig_prop = prop_land)], - dist_final[step == steps, .(CLASS, opt_prop = prop_land)], - by = "CLASS" -) - -final_compare[, change := opt_prop - orig_prop] - -ggplot(final_compare, - aes(x = reorder(CLASS, change), y = change)) + - geom_col() + - coord_flip() + - labs( - title = paste("Change in projected land share after", steps, "steps"), - x = "Class", - y = "Optimized - Original" - ) + - theme_minimal() - -#3.change in target class (isolated graph) -focus = dist_all[CLASS == target_crop] - -ggplot(focus, aes(x = year, y = prop_land, color = scenario)) + - geom_line(linewidth = 1.2) + - geom_point(size = 2) + - labs( - title = sprintf("Projected land share for class %s", target_crop), - x = "Year", - y = "Fraction of total land", - color = "Scenario" - ) + - theme_minimal() - -#-------------------sankey diagram-------------------------# -#where the 2023 land distribution ends up after the target number of steps -#under the optimized transition matrix - -A_use = A_final -scenario_name = "Optimized" - -#n-step transition matrix -A_steps = A_use %^% steps - -#flow table from 2023 class to final class -flows = as.data.table(as.table(A_steps)) -setnames(flows, c("source_class", "target_class", "transition_prob")) -flows[, source_land := as.numeric(X0[source_class])] -flows[, value := source_land * transition_prob] - -#**remove tiny flows so sankey is readable -flows = flows[value > 0.001] - -#*make separate node labels for start and end -flows[, source := paste0(source_class, " 2023")] -flows[, target := paste0(target_class, " ", 2023 + steps)] - -nodes = data.table(name = unique(c(flows$source, flows$target))) - -flows[, source_id := match(source, nodes$name) - 1] -flows[, target_id := match(target, nodes$name) - 1] - -links = flows[, .( - source = source_id, - target = target_id, - value = value -)] - -sankeyNetwork( - Links = links, - Nodes = nodes, - Source = "source", - Target = "target", - Value = "value", - NodeID = "name", - fontSize = 13, - nodeWidth = 25, - sinksRight = TRUE -) - -#####isolated sankey with only the target crop after time steps - -A_steps = A_final %^% steps - -flows = as.data.table(as.table(A_steps)) -setnames(flows, c("source_class", "target_class", "transition_prob")) - -flows[, source_land := as.numeric(X0[source_class])] -flows[, value := source_land * transition_prob] - -# only show land ending in target class -flows = flows[target_class == target_crop] - -flows[, source := paste0(source_class, " 2023")] -flows[, target := paste0(target_class, " ", 2023 + steps)] - -nodes = data.table(name = unique(c(flows$source, flows$target))) - -flows[, source_id := match(source, nodes$name) - 1] -flows[, target_id := match(target, nodes$name) - 1] - -links = flows[, .( - source = source_id, - target = target_id, - value = value -)] - -sankeyNetwork( - Links = links, - Nodes = nodes, - Source = "source", - Target = "target", - Value = "value", - NodeID = "name", - fontSize = 13, - nodeWidth = 25, - sinksRight = TRUE -) - -##-----predictions with the new optimized matrix----- - -tmat_year = A_final -states = rownames(tmat_year) - -year_states = fread("year_states.csv") -setDT(year_states) - -year_states[, parcel_id := as.character(parcel_id)] -year_states[, dominant_crop := trimws(as.character(dominant_crop))] - -start_info = year_states[, .SD[which.max(year)], by = parcel_id] - -design_points = fread('/projectnb/dietzelab/ccmmf/management/design_points_landiq_2018-2023.csv') -design_points[, parcel_id := as.character(parcel_id)] - -start_info = start_info[parcel_id %in% design_points$parcel_id] -year_states_design = year_states[parcel_id %in% design_points$parcel_id] + for (s in zero_rows) {A[s, ] = 0 + A[s, s] = 1}} + + row_sums = rowSums(A) + A = sweep(A, 1, row_sums, "/") + + return(A)} -end_year = 2050 +write_tmat = function(A, path) {out = as.data.table(A, keep.rownames = "state") +fwrite(out, path)} -all_preds_optimized = start_info[, { +build_x0_last_observed = function(all_data, county_name, start_year, states, state_col = "crop_class") { + dt = copy(all_data[county_safe == county_name & year <= start_year]) + if (nrow(dt) == 0) stop("No data found for county up to start_year: ", county_name, " / ", start_year) + + dt[, state_value := trimws(as.character(get(state_col)))] + latest_dt = dt[!is.na(state_value), .SD[which.max(year)], by = parcel_id] + latest_dt = latest_dt[state_value %in% states] + + x0_dt = latest_dt[, .(acres = sum(ACRES, na.rm = TRUE), n_parcels = uniqueN(parcel_id), + min_obs_year = min(year, na.rm = TRUE), max_obs_year = max(year, na.rm = TRUE)), + by = state_value] + + X0_vec = setNames(rep(0, length(states)), states) + matched_states = intersect(x0_dt$state_value, states) + X0_vec[matched_states] = x0_dt[match(matched_states, state_value), acres] + + X0_mat = matrix(X0_vec, nrow = 1) + colnames(X0_mat) = states + attr(X0_mat, "latest_observed_summary") = x0_dt + return(X0_mat)} + +scale_target_to_x0_total = function(target_vec, X0) { + target_total = sum(target_vec, na.rm = TRUE) + x0_total = sum(X0, na.rm = TRUE) + if (target_total <= 0 || x0_total <= 0) return(target_vec) + target_vec / target_total * x0_total} + +prep_target_for_opt = function(target_vec, X0, scale_to_x0_total = TRUE, nominal_zero_acres = 0.01) { + target_for_opt = target_vec + if (scale_to_x0_total) target_for_opt = scale_target_to_x0_total(target_for_opt, X0) + target_for_opt[!is.na(target_for_opt) & target_for_opt == 0] = nominal_zero_acres + if (scale_to_x0_total && sum(target_for_opt, na.rm = TRUE) > 0) { + target_for_opt = target_for_opt / sum(target_for_opt, na.rm = TRUE) * sum(X0, na.rm = TRUE)} + return(target_for_opt)} + +##-----optimizer----- +make_full_target_vec = function(raw_target_vec, states) { + out = setNames(rep(0, length(states)), states) + matched = intersect(names(raw_target_vec), states) + out[matched] = as.numeric(raw_target_vec[matched]) + return(out)} + +optimize_county_matrix = function(cty, A_orig, X0, target_vec, steps, + lambda_target = 1e6, target_vec_report = NULL, + maxeval = maxeval_optimizer, + maxtime = maxtime_optimizer) { + + states = rownames(A_orig) + n = length(states) + + #make sure target includes every crop state in the matrix + target_vec = make_full_target_vec(target_vec, states) + if (is.null(target_vec_report)) { + target_vec_report = target_vec + } else { + target_vec_report = make_full_target_vec(target_vec_report, states) + } - current_state = dominant_crop - years = seq(year + 1, end_year) + target_states = states - preds = character(length(years)) - probs = numeric(length(years)) + #pack only first n-1 columns of each row, last column is calculated as 1 - rowSum(first n-1). + pack_A = function(A) {as.vector(t(A[, 1:(n - 1), drop = FALSE]))} - for (i in seq_along(years)) { + unpack_x = function(x) { + A_part = matrix(x, nrow = n, ncol = n - 1, byrow = TRUE) + A_last = 1 - rowSums(A_part) + A_new = cbind(A_part, A_last) - p = tmat_year[current_state, ] + rownames(A_new) = states + colnames(A_new) = states - if (sum(p) == 0 || all(is.na(p))) { - preds[i] = NA_character_ - probs[i] = NA_real_ - next + return(A_new)} + + init_x = pack_A(A_orig) + + obj_fun = function(x) { + A_new = unpack_x(x) + + if (any(!is.finite(A_new)) || any(A_new < -1e-8) || any(A_new > 1 + 1e-8)) { + return(1e20) } - idx = which.max(p) - next_state = states[idx] + X_end = X0 %*% (A_new %^% steps) + colnames(X_end) = states - preds[i] = next_state - probs[i] = p[idx] + matrix_change_penalty = sum((A_new - A_orig)^2) - current_state = next_state - } + X_end_share = as.numeric(X_end[1, states]) / sum(X_end[1, states]) + target_share = as.numeric(target_vec[states]) / sum(target_vec[states]) + + target_error_penalty = sum((X_end_share - target_share)^2) + + matrix_change_penalty + lambda_target * target_error_penalty + } + #inequality constraint: sum(first n-1 row probs) <= 1, guarantees the last probability is >= 0. + constr_fun = function(x) {A_part = matrix(x, nrow = n, ncol = n - 1, byrow = TRUE) + rowSums(A_part) - 1} - .( - year = years, - pred_class = preds, - pred_prob = probs - ) + res = nloptr(x0 = init_x, eval_f = obj_fun, eval_g_ineq = constr_fun, + lb = rep(0, n * (n - 1)), ub = rep(1, n * (n - 1)), + opts = list(algorithm = "NLOPT_LN_COBYLA", xtol_rel = 1e-5, maxeval = maxeval, + maxtime = maxtime, print_level = 0)) -}, by = parcel_id] - -preds_future = all_preds_optimized[, .( - parcel_id, - year, - pred_class, - actual_class = NA_character_ -)] - -##-----make similar graphs using optimized scenarios----- -#(same script format as initial prediction code) - -actual_hist = year_states_design[, .( - parcel_id, - year, - pred_class = NA_character_, - actual_class = dominant_crop -)] - -plot_data_optimized = rbind(actual_hist, preds_future, fill = TRUE) - -sample_pids = sample(unique(plot_data_optimized$parcel_id), 1) - -plot_subset_optimized = plot_data_optimized[parcel_id %in% sample_pids] - -plot_subset_optimized[, pred_class := factor(pred_class, levels = states)] -plot_subset_optimized[, actual_class := factor(actual_class, levels = states)] - -ggplot(plot_subset_optimized, aes(x = year)) + - - geom_point( - data = plot_subset_optimized[!is.na(pred_class)], - aes(y = pred_class, color = pred_class), - size = 3 - ) + - - geom_point( - data = plot_subset_optimized[!is.na(actual_class)], - aes(y = actual_class), - shape = 1, - size = 3, - color = "black" - ) + - - facet_wrap(~parcel_id, ncol = 2) + - - scale_x_continuous( - limits = c(2018, end_year), - breaks = seq(2018, end_year, by = 2) - ) + - - labs( - title = sprintf( - "Optimized scenario predictions after actual crop classes for parcel %s", - sample_pids - ), - x = "Year", - y = "Crop class", - color = "Predicted class" - ) + - - theme_minimal() - -##-------store optimized predictions in landIQ style------- -#(still same script format as initial prediction code) - -design_points[, CLASS := as.character(CLASS)] -design_points[, SUBCLASS := as.character(SUBCLASS)] -preds_future[, parcel_id := as.character(parcel_id)] -preds_future[, pred_class := as.character(pred_class)] - -last_obs = design_points[ - order(year, season), - .SD[.N], - by = parcel_id -][ - , .( - parcel_id, - last_CLASS = CLASS, - last_SUBCLASS = SUBCLASS - ) -] - -subclass_probs = design_points[ - !is.na(CLASS) & !is.na(SUBCLASS), - .N, - by = .(CLASS, SUBCLASS) -] - -subclass_probs[, prob := N / sum(N), by = CLASS] - -draw_subclass = function(class_name) { + A_final = unpack_x(res$solution) - choices = subclass_probs[CLASS == class_name] + #clean optimizer numerical noise before writing matrix + A_final[is.na(A_final)] = 0 + A_final[!is.finite(A_final)] = 0 - if (nrow(choices) == 0) { - return(NA_character_) - } + if (any(A_final < 0, na.rm = TRUE)) { + warning(cty, " optimized matrix had negative values. Minimum = ", + min(A_final, na.rm = TRUE), ". Clamping negatives to 0.") + A_final[A_final < 0] = 0} - sample( - choices$SUBCLASS, - size = 1, - prob = choices$prob - ) -} - -preds_subclass_optimized = merge( - preds_future[, .(parcel_id, year, CLASS = pred_class)], - last_obs, - by = "parcel_id", - all.x = TRUE -) - -setorder(preds_subclass_optimized, parcel_id, year) - -preds_subclass_optimized[, SUBCLASS := NA_character_] - -for (p in unique(preds_subclass_optimized$parcel_id)) { + if (any(A_final > 1, na.rm = TRUE)) { + warning(cty, " optimized matrix had values > 1. Maximum = ", + max(A_final, na.rm = TRUE), ". Clamping values above 1 to 1.") + A_final[A_final > 1] = 1} - idxs = which(preds_subclass_optimized$parcel_id == p) + row_sums = rowSums(A_final) - prev_class = preds_subclass_optimized$last_CLASS[idxs[1]] - prev_subclass = preds_subclass_optimized$last_SUBCLASS[idxs[1]] + zero_rows = names(row_sums)[is.na(row_sums) | row_sums == 0] - for (i in idxs) { - - current_class = preds_subclass_optimized$CLASS[i] + if (length(zero_rows) > 0) { + warning(cty, " optimized matrix had zero rows after cleanup. Setting to self-loop: ", + paste(zero_rows, collapse = ", ")) - if (is.na(current_class)) { - preds_subclass_optimized$SUBCLASS[i] = NA_character_ - - } else if (!is.na(prev_class) && current_class == prev_class) { - - # if class does not change, keep previous subclass - preds_subclass_optimized$SUBCLASS[i] = prev_subclass - - } else { + for (s in zero_rows) {A_final[s, ] = 0 + A_final[s, s] = 1}} + + A_final = sweep(A_final, 1, rowSums(A_final), "/") + + rownames(A_final) = states + colnames(A_final) = states + + X_end_orig = X0 %*% (A_orig %^% steps) + X_end_final = X0 %*% (A_final %^% steps) + colnames(X_end_orig) = states + colnames(X_end_final) = states + + summary = data.table(county_safe = cty, target_state = target_states, + start_acres = as.numeric(X0[1, target_states]), + target_acres_raw = as.numeric(target_vec_report[target_states]), + target_acres_used_for_opt = as.numeric(target_vec[target_states]), + original_projected_acres = as.numeric(X_end_orig[1, target_states]), + optimized_projected_acres = as.numeric(X_end_final[1, target_states]), + raw_difference_after_optimization = + as.numeric(X_end_final[1, target_states]) - as.numeric(target_vec_report[target_states]), + opt_difference_after_optimization = + as.numeric(X_end_final[1, target_states]) - as.numeric(target_vec[target_states]), + optimizer_status = res$status, optimizer_message = res$message, + max_matrix_change = max(abs(A_final - A_orig)), + row_sum_error = max(abs(rowSums(A_final) - 1))) + + return(list(A_final = A_final, summary = summary, res = res))} +##-------scenario crop to transition-state mapping------- +scenario_crop_map_single = data.table( + Crop = c("All Other Berries", "Strawberries (Fresh Market)", "All Other Fruit Crops", + "All Other Nut Crops", "Almonds", "Pome Fruit", "Stone Fruit", "Citrus", + "Grapes Dried, Raisins", "Grapes, Table", "Grapes, Wine", "Fallow"), + crop_state = c("T", "T", "D", "D", "D", "D", "D", "C", "V", "V", "V", "X")) + +scenario_crop_map_single[, crop_key := normalize_crop_key(Crop)] + +scenario_crop_map_split = data.table( + Crop = c("All Other Field Crops (Incl. Pasture /Rangeland)", "Annual Cropland"), + split_group = c("field_pasture", "annual_cropland")) + +scenario_crop_map_split[, crop_key := normalize_crop_key(Crop)] + +get_split_states = function(split_group, crop_states) { + if (split_group == "field_pasture") return(intersect(c("F", "P"), crop_states)) + if (split_group == "annual_cropland") return(intersect(c("F", "G", "T", "R"), crop_states)) + return(character(0))} + +get_x0_split_weights = function(all_data, cty, start_year, split_states) { + dt = copy(all_data[county_safe == cty & year <= start_year]) + if (nrow(dt) == 0) { + weight = rep(1 / length(split_states), length(split_states)) + return(data.table(crop_state = split_states, split_weight = weight))} + + latest_dt = dt[!is.na(crop_class), .SD[which.max(year)], by = parcel_id] + latest_dt = latest_dt[crop_class %in% split_states] + + if (nrow(latest_dt) == 0) { + weight = rep(1 / length(split_states), length(split_states)) + return(data.table(crop_state = split_states, split_weight = weight))} + + out = latest_dt[, .(x0_acres = sum(ACRES, na.rm = TRUE)), by = crop_class] + out = merge(data.table(crop_state = split_states), out, by.x = "crop_state", by.y = "crop_class", all.x = TRUE) + out[is.na(x0_acres), x0_acres := 0] + + if (sum(out$x0_acres) == 0) out[, split_weight := 1 / .N] else out[, split_weight := x0_acres / sum(x0_acres)] + return(out[, .(crop_state, split_weight)])} + +expand_scenario_rows_to_crop_states = function(scenarios, all_data, cty, end_year, start_year, crop_states) { + scen_cty = copy(scenarios[county_safe == cty & Year == end_year]) + if (nrow(scen_cty) == 0) return(list(expanded = data.table(), unmatched = data.table())) + + scen_cty[, scenario_row_id := .I] + scen_cty[, crop_key := normalize_crop_key(Crop)] + + single = merge(scen_cty, scenario_crop_map_single[, .(crop_key, crop_state)], by = "crop_key", all.x = FALSE) + if (nrow(single) > 0) { + single[, split_group := "single"] + single[, split_weight := 1]} + + split_rows = merge(scen_cty, scenario_crop_map_split[, .(crop_key, split_group)], by = "crop_key", all.x = FALSE) + split_expanded_list = list() + + if (nrow(split_rows) > 0) { + for (sg in unique(split_rows$split_group)) { + rows_sg = split_rows[split_group == sg] + split_states = get_split_states(sg, crop_states) + if (length(split_states) == 0) { + warning("No valid split states found for split group: ", sg) + next + } - # if class changes, draw subclass from observed subclass distribution - preds_subclass_optimized$SUBCLASS[i] = draw_subclass(current_class) - } - - prev_class = current_class - prev_subclass = preds_subclass_optimized$SUBCLASS[i] + weights = get_x0_split_weights(all_data = all_data, cty = cty, start_year = start_year, split_states = split_states) + expanded = CJ(scenario_row_id = rows_sg$scenario_row_id, crop_state = weights$crop_state) + expanded = merge(expanded, rows_sg, by = "scenario_row_id", all.x = TRUE, allow.cartesian = TRUE) + expanded = merge(expanded, weights, by = "crop_state", all.x = TRUE) + split_expanded_list[[sg]] = expanded}} + + split_expanded = if (length(split_expanded_list) > 0) rbindlist(split_expanded_list, fill = TRUE) else data.table() + expanded = rbindlist(list(single, split_expanded), fill = TRUE) + + if (nrow(expanded) > 0) { + expanded = expanded[crop_state %in% crop_states] + expanded[, `:=`( + Acres_Total_mapped = Acres_Total * split_weight, + Tilled_acres_mapped = `Tilled acres` * split_weight, + Reduced_till_acres_mapped = `Reduced till acres (CPS 345)` * split_weight, + No_till_acres_mapped = `No till acres (CPS 329)` * split_weight)]} + + unmatched = scen_cty[ + !(crop_key %in% scenario_crop_map_single$crop_key) & + !(crop_key %in% scenario_crop_map_split$crop_key), + .(scenario_row_id, Crop, crop_key, Acres_Total)] + + if (nrow(unmatched) > 0) { + warning("Some scenario crops are not mapped: ", paste(unique(unmatched$Crop), collapse = ", "))} + + return(list(expanded = expanded, unmatched = unmatched))} + +build_scenario_crop_targets = function(scenarios, all_data, cty, end_year, start_year, crop_states) { + expanded_info = expand_scenario_rows_to_crop_states(scenarios = scenarios, all_data = all_data, cty = cty, + end_year = end_year, start_year = start_year, + crop_states = crop_states) + expanded = expanded_info$expanded + if (nrow(expanded) == 0) return(NULL) + + target_dt = expanded[, .(target_acres_raw = sum(Acres_Total_mapped, na.rm = TRUE), + scenario_crops = paste(sort(unique(Crop)), collapse = "; "), + n_scenario_rows = uniqueN(scenario_row_id)), by = crop_state] + target_dt = target_dt[crop_state %in% crop_states] + if (nrow(target_dt) == 0) return(NULL) + + target_vec = setNames(target_dt$target_acres_raw, target_dt$crop_state) + return(list(target_vec = target_vec, target_dt = target_dt, + unmatched = expanded_info$unmatched, expanded_rows = expanded))} + +make_crop_by_till_targets = function(expanded_rows) { + no_till_targets = expanded_rows[, .(target_acres_raw = sum(No_till_acres_mapped, na.rm = TRUE), + scenario_crops = paste(sort(unique(Crop)), collapse = "; ")), + by = .(county_safe, Year, crop_state)] + no_till_targets[, till_state := "no_till"] + + low_till_targets = expanded_rows[, .(target_acres_raw = sum(Reduced_till_acres_mapped, na.rm = TRUE), + scenario_crops = paste(sort(unique(Crop)), collapse = "; ")), + by = .(county_safe, Year, crop_state)] + low_till_targets[, till_state := "low_till"] + + high_till_targets = expanded_rows[, .(target_acres_raw = sum(Tilled_acres_mapped, na.rm = TRUE), + scenario_crops = paste(sort(unique(Crop)), collapse = "; ")), + by = .(county_safe, Year, crop_state)] + high_till_targets[, till_state := "high_till"] + + out = rbindlist(list(no_till_targets, low_till_targets, high_till_targets), fill = TRUE) + setcolorder(out, c("county_safe", "Year", "crop_state", "till_state", "target_acres_raw", "scenario_crops")) + return(out)} + +##-------validation helpers------- +check_matrix = function(A, matrix_name = "matrix") { + out = data.table(matrix_name = matrix_name, min_value = min(A, na.rm = TRUE), + max_value = max(A, na.rm = TRUE), + min_row_sum = min(rowSums(A), na.rm = TRUE), + max_row_sum = max(rowSums(A), na.rm = TRUE), + max_row_sum_error = max(abs(rowSums(A) - 1), na.rm = TRUE), + any_negative = any(A < -1e-10, na.rm = TRUE), + any_over_one = any(A > 1 + 1e-10, na.rm = TRUE)) + print(out) + if (out$any_negative) warning(matrix_name, " has negative probabilities.") + if (out$any_over_one) warning(matrix_name, " has probabilities > 1.") + if (out$max_row_sum_error > 1e-6) warning(matrix_name, " rows do not sum to 1.") + return(out)} + +check_mapping_totals = function(raw_rows, expanded_rows) { + raw_crop_total = sum(raw_rows$Acres_Total, na.rm = TRUE) + mapped_crop_total = sum(expanded_rows$Acres_Total_mapped, na.rm = TRUE) + out = data.table(raw_crop_total = raw_crop_total, mapped_crop_total = mapped_crop_total, + crop_total_diff = mapped_crop_total - raw_crop_total) + print(out) + if (abs(out$crop_total_diff) > 1e-4) warning("Mapped crop acres do not equal raw scenario crop acres.") + return(out)} + +##-------load & clean all_data------- +all_data = fread("all_data.csv") +if ("V1" %in% names(all_data)) all_data[, V1 := NULL] +setDT(all_data) + +required_all_data_cols = c("parcel_id", "year", "county", "crop_class", "ACRES") +check_required_cols(all_data, required_all_data_cols, "all_data") + +all_data[, `:=`(parcel_id = as.character(parcel_id), + year = as.integer(year), county = as.character(county), + crop_class = trimws(as.character(crop_class)), ACRES = as.numeric(ACRES))] +all_data[, county_safe := safe_county_name(county)] + +##-------load & clean scenario sheet------- +scenarios = as.data.table(read_excel("MAGiC scenarios_FINAL.xlsx", sheet = scenario_name)) +setnames(scenarios, names(scenarios), trimws(names(scenarios))) + +required_scenario_cols = c("Crop", "County", "Year", "Acres_Total", + "Tilled acres", "Reduced till acres (CPS 345)", "No till acres (CPS 329)") +check_required_cols(scenarios, required_scenario_cols, "scenario sheet") + +scenarios[, `:=`(Crop = trimws(as.character(Crop)), County = trimws(as.character(County)), + Year = as.integer(Year), Acres_Total = as.numeric(Acres_Total), + `Tilled acres` = as.numeric(`Tilled acres`), + `Reduced till acres (CPS 345)` = as.numeric(`Reduced till acres (CPS 345)`), + `No till acres (CPS 329)` = as.numeric(`No till acres (CPS 329)`))] + +scenarios[, county_safe := safe_county_name(County)] + +##-------run one county------- +run_county = function(focus_county) { + focus_county_safe = safe_county_name(focus_county) + message("\n==============================") + message("Running county: ", focus_county_safe) + message("==============================") + + if (!(focus_county_safe %in% scenarios$county_safe)) stop("County not found in scenario sheet: ", focus_county_safe) + if (!(focus_county_safe %in% all_data$county_safe)) stop("County not found in all_data: ", focus_county_safe) + + county_year_acres = all_data[county_safe == focus_county_safe, + .(n_rows = .N, n_parcels = uniqueN(parcel_id), + total_acres = sum(ACRES, na.rm = TRUE)), by = year][order(year)] + print(county_year_acres) + + county_scenario_rows = scenarios[county_safe == focus_county_safe & Year == end_year] + if (nrow(county_scenario_rows) == 0) stop("No scenario rows found for county/year: ", focus_county_safe, " / ", end_year) + + raw_target_path = file.path(target_out_dir, paste0("raw_scenario_rows_", focus_county_safe, "_", scenario_name, "_", end_year, ".csv")) + fwrite(county_scenario_rows, raw_target_path) + + crop_matrix_file = file.path("county_crop_matrices", paste0(focus_county_safe, "_crop_matrix.csv")) + if (!file.exists(crop_matrix_file)) stop("Crop matrix file not found: ", crop_matrix_file) + + A_crop_orig = repair_transition_matrix(read_tmat(crop_matrix_file), paste0("crop matrix ", focus_county_safe)) + crop_states = rownames(A_crop_orig) + + crop_target = build_scenario_crop_targets(scenarios = scenarios, all_data = all_data, cty = focus_county_safe, + end_year = end_year, start_year = start_year, crop_states = crop_states) + if (is.null(crop_target)) stop("No crop target vector could be built for county: ", focus_county_safe) + + crop_target_vec_raw = make_full_target_vec(crop_target$target_vec, crop_states) + X0_crop = build_x0_last_observed(all_data = all_data, county_name = focus_county_safe, + start_year = start_year, states = crop_states, state_col = "crop_class") + if (sum(X0_crop, na.rm = TRUE) <= 0) stop("X0 crop total is zero for county: ", focus_county_safe) + + crop_target_vec_opt = prep_target_for_opt(target_vec = crop_target_vec_raw, X0 = X0_crop, + scale_to_x0_total = scale_crop_targets_to_x0, + nominal_zero_acres = nominal_zero_acres) + + crop_targets_out = copy(crop_target$target_dt) + crop_targets_out[, `:=`(scenario_name = scenario_name, county_safe = focus_county_safe, + start_year = start_year, end_year = end_year, + target_state_col = "manual_scenario_crop_to_landiq_class_map", + target_acres_used_for_opt = as.numeric(crop_target_vec_opt[crop_state]))] + + crop_targets_path = file.path(target_out_dir, paste0("crop_targets_", focus_county_safe, "_", scenario_name, "_", end_year, ".csv")) + expanded_rows_path = file.path(target_out_dir, paste0("expanded_scenario_rows_", focus_county_safe, "_", scenario_name, "_", end_year, ".csv")) + fwrite(crop_targets_out, crop_targets_path) + fwrite(crop_target$expanded_rows, expanded_rows_path) + + crop_by_till_targets = make_crop_by_till_targets(crop_target$expanded_rows) + crop_by_till_targets_path = file.path(target_out_dir, paste0("crop_by_till_targets_", focus_county_safe, "_", scenario_name, "_", end_year, ".csv")) + fwrite(crop_by_till_targets, crop_by_till_targets_path) + + crop_by_till_check = crop_by_till_targets[, .(till_target_total = sum(target_acres_raw, na.rm = TRUE)), by = crop_state] + crop_by_till_check = merge(crop_by_till_check, crop_targets_out[, .(crop_state, crop_target_total = target_acres_raw)], + by = "crop_state", all = TRUE) + crop_by_till_check[, diff := till_target_total - crop_target_total] + + crop_by_till_check_path = file.path(target_out_dir, paste0("crop_by_till_check_", focus_county_safe, "_", scenario_name, "_", end_year, ".csv")) + fwrite(crop_by_till_check, crop_by_till_check_path) + + if (nrow(crop_target$unmatched) > 0) { + unmatched_path = file.path(target_out_dir, paste0("unmatched_scenario_crops_", focus_county_safe, "_", scenario_name, "_", end_year, ".csv")) + fwrite(crop_target$unmatched, unmatched_path) + warning("Wrote unmatched scenario crops to: ", unmatched_path) } + + scenario_total = scenarios[county_safe == focus_county_safe & Year == end_year, sum(Acres_Total, na.rm = TRUE)] + acreage_basis_check = data.table(county_safe = focus_county_safe, start_year = start_year, end_year = end_year, + x0_latest_observed_acres = sum(X0_crop), + scenario_target_acres = scenario_total, + scenario_minus_x0 = scenario_total - sum(X0_crop), + scenario_divided_by_x0 = scenario_total / sum(X0_crop), + target_acres_used_for_opt_total = sum(crop_target_vec_opt)) + print(acreage_basis_check) + + mapping_total_check = check_mapping_totals(raw_rows = county_scenario_rows, expanded_rows = crop_target$expanded_rows) + split_weight_check = crop_target$expanded_rows[, .(weight_sum = sum(split_weight, na.rm = TRUE)), by = scenario_row_id] + bad_weights = split_weight_check[abs(weight_sum - 1) > 1e-6] + if (nrow(bad_weights) > 0) warning("Some scenario split weights do not sum to 1 for county: ", focus_county_safe) + + orig_matrix_check = check_matrix(A_crop_orig, paste0("original crop matrix ", focus_county_safe)) + + message("Starting optimizer for: ", focus_county_safe) + opt_start_time = Sys.time() + + opt_crop = optimize_county_matrix(cty = focus_county_safe, A_orig = A_crop_orig, X0 = X0_crop, + target_vec = crop_target_vec_opt, steps = steps, + lambda_target = lambda_target, target_vec_report = crop_target_vec_raw, + maxeval = maxeval_optimizer, maxtime = maxtime_optimizer) + + message("Finished optimizer for: ", focus_county_safe, + " in ", round(as.numeric(difftime(Sys.time(), opt_start_time, units = "mins")), 2), " minutes") + + crop_out_path = file.path(scenario_out_dir, paste0(focus_county_safe, "_crop_matrix_", scenario_name, ".csv")) + write_tmat(opt_crop$A_final, crop_out_path) + + optimization_summary = opt_crop$summary + optimization_summary[, `:=`(scenario_name = scenario_name, matrix_type = "crop", focus_group = "crop_class", + focus_county = focus_county, start_year = start_year, end_year = end_year, + x0_rule = "latest_observed_crop_state_per_parcel_up_to_start_year", + target_note = ifelse(scale_crop_targets_to_x0, + "Crop targets scaled to X0 total for feasible row-stochastic optimization", + "Raw crop targets used directly"))] + + optimization_summary[, abs_error_opt := optimized_projected_acres - target_acres_used_for_opt] + optimization_summary[, pct_error_opt := abs_error_opt / pmax(abs(target_acres_used_for_opt), 1)] + optimization_summary[, abs_error_raw := optimized_projected_acres - target_acres_raw] + optimization_summary[, pct_error_raw := abs_error_raw / pmax(abs(target_acres_raw), 1)] + + total_opt_error_share = sum(abs( + optimization_summary$optimized_projected_acres / sum(optimization_summary$optimized_projected_acres) - + optimization_summary$target_acres_used_for_opt / sum(optimization_summary$target_acres_used_for_opt) + ), na.rm = TRUE) + + true_run_status = ifelse(opt_crop$res$status < 0, + "optimizer_failed", + ifelse(total_opt_error_share > 0.05, "poor_fit", "success")) + + summary_path = file.path(scenario_out_dir, paste0("optimization_summary_", focus_county_safe, "_", scenario_name, ".csv")) + fit_check_path = file.path(scenario_out_dir, paste0("optimization_fit_check_", focus_county_safe, "_", scenario_name, ".csv")) + fwrite(optimization_summary, summary_path) + fwrite(optimization_summary, fit_check_path) + + opt_matrix_check = check_matrix(opt_crop$A_final, paste0("optimized crop matrix ", focus_county_safe)) + + manifest = data.table(run_status = true_run_status, error_message = ifelse(true_run_status == "success", NA_character_, opt_crop$res$message), + scenario_name = scenario_name, focus_county = focus_county, + focus_county_safe = focus_county_safe, start_year = start_year, + end_year = end_year, steps = steps, + x0_total_acres = sum(X0_crop), scenario_target_acres = scenario_total, + target_acres_used_for_opt_total = sum(crop_target_vec_opt), + raw_scenario_rows_path = raw_target_path, + expanded_scenario_rows_path = expanded_rows_path, + crop_targets_path = crop_targets_path, + crop_by_till_targets_path = crop_by_till_targets_path, + crop_by_till_check_path = crop_by_till_check_path, + optimized_crop_matrix_path = crop_out_path, + optimization_summary_path = summary_path, + optimization_fit_check_path = fit_check_path, + x0_rule = "latest_observed_crop_state_per_parcel_up_to_start_year", + max_matrix_change = max(abs(opt_crop$A_final - A_crop_orig)), + row_sum_error = max(abs(rowSums(opt_crop$A_final) - 1)), + total_opt_error_share = total_opt_error_share) + + manifest_path = file.path(scenario_out_dir, paste0("run_manifest_", focus_county_safe, "_", scenario_name, ".csv")) + fwrite(manifest, manifest_path) + + message("Finished county: ", focus_county_safe) + return(manifest)} + +##-----all county loop----- +if (run_all_counties) { + counties_to_run = sort(intersect(unique(scenarios$county_safe), unique(all_data$county_safe))) +} else { + counties_to_run = safe_county_name(counties_manual) } -#add class/subclass descriptions -lookup = fread('/projectnb/dietzelab/ccmmf/management/LandIQ_cropCode_lookup_table.csv') - -lookup[, CLASS := as.character(CLASS)] -lookup[, SUBCLASS := as.character(SUBCLASS)] - -lookup_subclass = unique(lookup[, .( - CLASS, - SUBCLASS, - CLASS_desc, - SUBCLASS_desc, - PFT -)]) - -preds_subclass_optimized = merge( - preds_subclass_optimized, - lookup_subclass, - by = c("CLASS", "SUBCLASS"), - all.x = TRUE -) - -parcel_meta = design_points[ - order(year, season), - .SD[.N], - by = parcel_id -][ - , .(parcel_id, site_id, lon, lat) -] - -preds_landiq_optimized = merge( - preds_subclass_optimized, - parcel_meta, - by = "parcel_id", - all.x = TRUE -) - -preds_landiq_optimized[, season := NA_integer_] - -# keep LandIQ-style columns -preds_landiq_optimized = preds_landiq_optimized[, .( - site_id, - parcel_id, - lon, - lat, - year, - season, - CLASS, - SUBCLASS, - CLASS_desc, - SUBCLASS_desc, - PFT -)] - -design_points_predicted_optimized = copy( - preds_landiq_optimized[year >= 2024 & year <= end_year] -) - -design_points_predicted_optimized[, source := "predicted_optimized"] - -setorder(design_points_predicted_optimized, parcel_id, year, season) - -fwrite( - design_points_predicted_optimized, - sprintf( - "predicted_optimized_%s_2024_%s.csv", - target_crop, - end_year +all_manifests = rbindlist(lapply(counties_to_run, function(cty) { + tryCatch( + run_county(cty), + error = function(e) { + data.table( + run_status = "error", + error_message = conditionMessage(e), + scenario_name = scenario_name, + focus_county = cty, + focus_county_safe = safe_county_name(cty), + start_year = start_year, + end_year = end_year + ) + } ) -) \ No newline at end of file +}), fill = TRUE) + +all_manifest_path = file.path(scenario_out_dir, paste0("all_county_run_manifest_", scenario_name, ".csv")) + +fwrite(all_manifests, all_manifest_path) +print(all_manifests) diff --git a/modules/data.remote/inst/transition_matrix.R b/modules/data.remote/inst/transition_matrix.R index f07d3cf3daa..0a97f4de6b8 100644 --- a/modules/data.remote/inst/transition_matrix.R +++ b/modules/data.remote/inst/transition_matrix.R @@ -1,4 +1,4 @@ -##creates the transition matricies for all parcels in each county by creating crop class sequences +##creates the transition matrices for all parcels in each county by creating crop class sequences setwd("/projectnb/dietzelab/ananyak") library(arrow) @@ -15,73 +15,52 @@ ag_classes = unique(lookup[is_agricultural == TRUE, as.character(CLASS)]) year_min = 2018L year_max = 2023L -crops_full <- as.data.table( +crops_full = as.data.table( arrow::open_dataset(file.path(path_landiq_v4, "crops_all_years.parq")) |> filter(year >= year_min, year <= year_max, CLASS %in% ag_classes) |> - select(parcel_id, year, season, CLASS, SUBCLASS, centx, centy) |> - collect() -) + select(parcel_id, year, season, CLASS, SUBCLASS, centx, centy, ACRES) |> + collect()) crops_full[, `:=`( parcel_id = as.character(parcel_id), year = as.integer(year), season = as.integer(season), CLASS = as.character(CLASS), - SUBCLASS = as.character(SUBCLASS) -)] + ACRES = as.integer(ACRES), + SUBCLASS = as.character(SUBCLASS))] ##-----adding county to crops file----- library(sf) library(tigris) options(tigris_use_cache = TRUE) -#one row per parcel so repeated years/seasons do not duplicate spatial join +#one row per parcel so repeated years/seasons dont duplicate spatial join parcel_unique = crops_full[ - !is.na(centx) & !is.na(centy), - .SD[1], - by = parcel_id] + !is.na(centx) & !is.na(centy), .SD[1], by = parcel_id] -#convert centx/centy to sf; centx/centy in EPSG:3310 -parcel_sf = st_as_sf( - parcel_unique, - coords = c("centx", "centy"), - crs = 3310, - remove = FALSE) +#convert centx/centy to sf; and change from EPSG:3310 to 4326 +parcel_sf = st_as_sf(parcel_unique, coords = c("centx", "centy"), crs = 3310, remove = FALSE) #California county boundaries ca_counties = counties(state = "CA", cb = TRUE, class = "sf") |> st_transform(4326) ca_counties = ca_counties |> st_transform(3310) -#spatial join parcels to counties -parcel_county = st_join( - parcel_sf, - ca_counties[, c("NAME", "GEOID")]) +parcel_county = st_join(parcel_sf, ca_counties[, c("NAME", "GEOID")]) -#convert back to data.table and rename columns parcel_county_dt = as.data.table(st_drop_geometry(parcel_county)) setnames(parcel_county_dt, "NAME", "county") -#keep only parcel_id + county info parcel_county_lookup = parcel_county_dt[, .( - parcel_id, - county, - county_geoid = GEOID -)] - -# merge county back onto the full crops_full table -crops_full_county = merge( - crops_full, - parcel_county_lookup, - by = "parcel_id", - all.x = TRUE -) - -##-------order for sequence building - add season if needed------- + parcel_id, county, county_geoid = GEOID)] + +#merge county back onto the full crops_full table +crops_full_county = merge(crops_full, parcel_county_lookup, by = "parcel_id", all.x = TRUE) + +##-------order for sequence building------- setorder(crops_full_county, parcel_id, county, year, season) write.csv(crops_full_county, 'crops_full_counties.csv') - ##-------cleaning rules to lower the amount of messy/unrealistic sequences (for 'X' cases)------ fix_seq = function(seq) { parts = strsplit(seq, "-", fixed = TRUE)[[1]] @@ -95,9 +74,7 @@ fix_seq = function(seq) { for (i in 2:(n - 1)) { if (parts[i] == "X" && parts[i - 1] == parts[i + 1]) { parts[i] = parts[i - 1] - } - } - } + }}} #Rule 2: similar to rule 2 but with a 'short run' of X's i = 1 @@ -108,8 +85,7 @@ fix_seq = function(seq) { end = i - 1 run_len = end - start + 1 - if ( - run_len <= 2 && + if (run_len <= 2 && start > 1 && end < n && parts[start - 1] == parts[end + 1] @@ -118,8 +94,7 @@ fix_seq = function(seq) { } } else { i = i + 1 - } - } + }} #Rule 3: 'edge X's' where its all one class + one X -- replace with the majority if (n >= 2) { @@ -128,8 +103,7 @@ fix_seq = function(seq) { } if (parts[n] == "X") { parts[n] = parts[n - 1] - } - } + }} #Rule 4: remaining short X run with one valid neighbor side --> fill from that side i = 1 @@ -152,22 +126,16 @@ fix_seq = function(seq) { } } else { i = i + 1 - } - } - paste(parts, collapse = "-") -} - - + }} + paste(parts, collapse = "-")} ##-----sequence formatting by parcel----- crop_sequences = crops_full_county[ , .( - crop_sequence = paste(CLASS, collapse = "-"), - season_sequence = paste(season, collapse = "-") + crop_sequence = paste(CLASS, collapse = "-"), season_sequence = paste(season, collapse = "-") ), - by = .(county, county_geoid, parcel_id, year) -] + by = .(county, county_geoid, parcel_id, year, ACRES)] ##-------apply sequence-fixing rules------- crop_sequences[, crop_sequence := vapply(crop_sequence, fix_seq, character(1))] @@ -203,8 +171,7 @@ seq_lookup[, c("dominant_crop", "non_dom_prob") := { } } - .(dom, prob) -}] + .(dom, prob)}] crop_sequences = seq_lookup[crop_sequences, on = c("crop_sequence", "season_sequence")] @@ -212,28 +179,21 @@ crop_sequences = seq_lookup[crop_sequences, on = c("crop_sequence", "season_sequ ##parcel_id, year, state, non_dom_prob, optional grouping columns like county year_states = copy(crop_sequences)[,.( - county, - county_geoid, - parcel_id, - year, - state = dominant_crop, - non_dom_prob)] + county, county_geoid, parcel_id, year, state = dominant_crop, ACRES, non_dom_prob)] year_states[, state := trimws(as.character(state))] year_states[, parcel_id := as.character(parcel_id)] year_states[, year := as.integer(year)] -setorder(year_states, county, parcel_id, year) +year_states[, ACRES := as.integer(ACRES)] +setorder(year_states, county, parcel_id, year) +fwrite(year_states, "/projectnb/dietzelab/ananyak/crop_year_states_cleaned.csv") -##-----function for transition format------ -make_transitions = function( - year_states, - id_col = "parcel_id", - time_col = "year", - state_col = "state", - non_dom_col = "non_dom_prob", - min_weight = 0.05) { +##-----function for transition formatting------ +make_transitions = function(year_states, id_col = "parcel_id", time_col = "year", + state_col = "state", non_dom_col = "non_dom_prob", + min_weight = 0.05) { dt = copy(as.data.table(year_states)) @@ -261,23 +221,16 @@ make_transitions = function( !is.na(to) & next_time == time + 1] - transitions[, weight := pmax( - min_weight, - (1 - from_non_dom) * (1 - to_non_dom))] + transitions[, weight := pmax(min_weight, (1 - from_non_dom) * (1 - to_non_dom))] setnames(transitions, "id", id_col) setnames(transitions, "time", time_col) - return(transitions) -} + return(transitions)} ##-----function to make a transition matrix------ -make_transition_matrix = function( - dt, - states_all, - from_col = "from", - to_col = "to", - weight_col = "weight") { +make_transition_matrix = function(dt, states_all, from_col = "from", to_col = "to", + weight_col = "weight") { dt = copy(as.data.table(dt)) @@ -295,25 +248,18 @@ make_transition_matrix = function( by = .(from, to)] if (nrow(transitions_weighted) == 0) { - empty_mat = matrix( - 0, - nrow = length(states_all), - ncol = length(states_all), - dimnames = list(states_all, states_all)) + empty_mat = matrix(0, nrow = length(states_all), ncol = length(states_all), + dimnames = list(states_all, states_all)) return(empty_mat)} - tmat_counts = dcast( - transitions_weighted, - from ~ to, - value.var = "N", - fill = 0) + tmat_counts = dcast(transitions_weighted, from ~ to, value.var = "N", fill = 0) - ## add missing columns + #add missing columns missing_cols = setdiff(states_all, colnames(tmat_counts)) for (mc in missing_cols) { tmat_counts[[mc]] = 0} - ## add missing rows + #add missing rows missing_rows = setdiff(states_all, tmat_counts$from) if (length(missing_rows) > 0) { zero_rows = data.table(from = missing_rows) @@ -321,13 +267,13 @@ make_transition_matrix = function( zero_rows[[s]] = 0} tmat_counts = rbind(tmat_counts, zero_rows, fill = TRUE)} - ## order rows/cols + #order rows/cols tmat_counts[, ord := match(from, states_all)] setorder(tmat_counts, ord) tmat_counts[, ord := NULL] tmat_counts = tmat_counts[, c("from", states_all), with = FALSE] - ## convert to probability matrix + #convert to probability matrix rn = tmat_counts$from prob_mat = as.matrix(tmat_counts[, ..states_all]) storage.mode(prob_mat) = "double" @@ -350,51 +296,32 @@ make_transition_matrix = function( return(tmat_final)} ##-----function to make transition matrices by a category (crop, county, etc)------ -#avoids having to use split(... by=grouping) -make_grouped_transition_matrices = function( - transitions, - states_all, - group_cols) { - - transition_groups = split( - transitions, - by = group_cols, - keep.by = TRUE) - - transition_mats = lapply( - transition_groups, - make_transition_matrix, - states_all = states_all) +make_grouped_transition_matrices = function(transitions, states_all, group_cols) { + + transition_groups = split(transitions, by = group_cols, keep.by = TRUE) + + transition_mats = lapply(transition_groups, make_transition_matrix, states_all = states_all) return(transition_mats)} - -##-----using functions to make matrices by county------ -transitions_full = make_transitions( - year_states = year_states, - id_col = "parcel_id", - time_col = "year", - state_col = "state", - non_dom_col = "non_dom_prob") +##-----using functions to make crop transition matrices by county------ +transitions_full = make_transitions(year_states = year_states, id_col = "parcel_id", + time_col = "year", state_col = "state", + non_dom_col = "non_dom_prob") states_all = c("YP","D","X","T","G","F","P","C","I","V","R") -county_transition_mats = make_grouped_transition_matrices( - transitions = transitions_full, - states_all = states_all, - group_cols = "county") - +county_transition_mats = make_grouped_transition_matrices(transitions = transitions_full, + states_all = states_all, + group_cols = "county") ##-----save------ -dir.create("county_transition_matrices", showWarnings = FALSE) +dir.create("county_crop_matrices", showWarnings = FALSE) for (cty in names(county_transition_mats)) { safe_name = gsub("[^A-Za-z0-9_]+", "_", cty) - write.csv( - county_transition_mats[[cty]], - file = file.path( - "county_transition_matrices", - paste0(safe_name, "_transition_matrix.csv")), - row.names = TRUE)} \ No newline at end of file + write.csv(county_transition_mats[[cty]], + file = file.path("county_crop_matrices", paste0(safe_name, "_crop_matrix.csv")), + row.names = TRUE)} \ No newline at end of file From c53dcc6839f40d38cff6f2784f9f7ab5fef13c1b Mon Sep 17 00:00:00 2001 From: Ananya Kulkarni Date: Mon, 6 Jul 2026 16:40:30 -0400 Subject: [PATCH 2/3] script that assigns predicted planting and harvesting dates per county --- .../data.remote/inst/planting_harvest_dates.R | 310 ++++++++++++++++++ 1 file changed, 310 insertions(+) create mode 100644 modules/data.remote/inst/planting_harvest_dates.R diff --git a/modules/data.remote/inst/planting_harvest_dates.R b/modules/data.remote/inst/planting_harvest_dates.R new file mode 100644 index 00000000000..f5727d22f4f --- /dev/null +++ b/modules/data.remote/inst/planting_harvest_dates.R @@ -0,0 +1,310 @@ +setwd("/projectnb/dietzelab/ananyak") + +library(data.table) +library(arrow) +library(sf) +library(tigris) + +options(tigris_use_cache = TRUE) +options(arrow.unsafe_metadata = TRUE) + +phenology_dir = "/projectnb/dietzelab/ccmmf/management/phenology/matched_landiq_mslsp_v4.1" + +##-----loading files----- +#assigned parquets +phenology_files = list.files(phenology_dir, pattern = "^assigned_year=.*\\.parquet$", + full.names = TRUE) + +#combining them all +assigned_all = rbindlist(lapply(phenology_files, function(f) { + dt = as.data.table(read_parquet(f)) + #pull year from filename if year column is not already there + yr = as.integer(gsub(".*assigned_year=([0-9]{4})\\.parquet$", "\\1", f)) + dt[, source_year := yr] + + return(dt)}), + fill = TRUE) + +##-----inspect----- +names(assigned_all) + +View(assigned_all[, + .( + parcel_id, source_year, season, mslsp_cycle, landiq_CLASS, landiq_SUBCLASS, landiq_ADOY, mslsp_OGI, + mslsp_OGMn, mslsp_OGD, qc_adoy_vs_cycle, qc_cycle_status)]) + +##-----merge phenology file info with crops&county info----- +crops_and_counties = fread("/projectnb/dietzelab/ananyak/crops_full_counties.csv") + +if ("V1" %in% names(crops_and_counties)) {crops_and_counties[, V1 := NULL]} + +parcel_county_lookup = unique(crops_and_counties[!is.na(county), .(parcel_id, county, county_geoid)]) + +assigned_all[, parcel_id := as.character(parcel_id)] +parcel_county_lookup[, parcel_id := as.character(parcel_id)] + +assigned_county = merge(assigned_all, parcel_county_lookup, by = "parcel_id", all.x = TRUE) + +##-----convert Date columns to day-of-year----- +assigned_county[, planting_doy := as.integer(format(mslsp_OGI, "%j"))] +assigned_county[, harvest_doy := as.integer(format(mslsp_OGMn, "%j"))] +assigned_county[, peak_doy := as.integer(format(mslsp_Peak, "%j"))] + +#most likely already numeric but set just incase +assigned_county[, landiq_ADOY := as.numeric(landiq_ADOY)] + +##-----helper function that doesnt wrap if the dates being used go through new years ----- +#normal stats are skewed when dates cluster around New Year - ex) 365 and 5 are actually close but +#numerically far apart + +fix_wrap_doy = function(x, low_cutoff = 45, high_cutoff = 320) { + x = x[!is.na(x)] + + if (length(x) == 0) return(x) + + wraps = quantile(x, 0.05, na.rm = TRUE) <= low_cutoff && + quantile(x, 0.95, na.rm = TRUE) >= high_cutoff + + if (wraps) {x = fifelse(x <= low_cutoff, x + 365, x)} #treat january dates as after december dates + + return(x)} + +wrap_back_doy = function(x) {((round(x) - 1) %% 365) + 1} + +##-----then can build table with means and sd's----- +assigned_county[, crop_type := fifelse( + !is.na(landiq_SUBCLASS) & landiq_SUBCLASS != "", + paste0(landiq_CLASS, "_", landiq_SUBCLASS), + as.character(landiq_CLASS))] + +mvp = assigned_county[ + !is.na(county) & + !is.na(season) & + !is.na(crop_type) & + source_year >= 2018 & + source_year <= 2023, + { + planting_adj = fix_wrap_doy(planting_doy) + harvest_adj = fix_wrap_doy(harvest_doy) + peak_adj = fix_wrap_doy(peak_doy) + + .( + n_rows = .N, n_planting = sum(!is.na(planting_doy)), n_harvest = sum(!is.na(harvest_doy)), n_peak = sum(!is.na(peak_doy)), + + planting_wraparound = length(planting_adj) > 0 && max(planting_adj, na.rm = TRUE) > 365, + harvest_wraparound = length(harvest_adj) > 0 && max(harvest_adj, na.rm = TRUE) > 365, + + mean_planting_doy = wrap_back_doy(mean(planting_adj, na.rm = TRUE)), + sd_planting_doy = sd(planting_adj, na.rm = TRUE), + q05_planting_doy = wrap_back_doy(quantile(planting_adj, 0.05, na.rm = TRUE)), + q50_planting_doy = wrap_back_doy(quantile(planting_adj, 0.50, na.rm = TRUE)), + q95_planting_doy = wrap_back_doy(quantile(planting_adj, 0.95, na.rm = TRUE)), + + mean_harvest_doy = wrap_back_doy(mean(harvest_adj, na.rm = TRUE)), + sd_harvest_doy = sd(harvest_adj, na.rm = TRUE), + q05_harvest_doy = wrap_back_doy(quantile(harvest_adj, 0.05, na.rm = TRUE)), + q50_harvest_doy = wrap_back_doy(quantile(harvest_adj, 0.50, na.rm = TRUE)), + q95_harvest_doy = wrap_back_doy(quantile(harvest_adj, 0.95, na.rm = TRUE)), + + mean_peak_doy = wrap_back_doy(mean(peak_adj, na.rm = TRUE)), + sd_peak_doy = sd(peak_adj, na.rm = TRUE))}, + by = .(county, county_geoid, season, crop_type, landiq_CLASS, landiq_SUBCLASS)] + +setorder(mvp, county, season, crop_type) +#View(mvp) +#write.csv(mvp, 'mvp_draft1.csv') + +pheno = assigned_all[ + year >= 2018 & year <= 2023 & + assigned_by == "matched" & + !is.na(mslsp_OGI) & + !is.na(mslsp_OGMn) & + !is.na(parcel_id), + .( + parcel_id = as.character(parcel_id), year = as.integer(year), season, planting_date = as.IDate(mslsp_OGI), + harvest_date = as.IDate(mslsp_OGMn), peak_date = as.IDate(mslsp_Peak), landiq_CLASS,landiq_SUBCLASS)] + + +##-----assign all BAU and NBS predicted files planting and harvesting dates----- +#loop through the two subfolders in county_landiq_predictions and assign mean dates to each +#countys original file + +#helper: modal value +mode_value = function(x) { + x = x[!is.na(x)] + if (length(x) == 0) return(NA) + names(sort(table(x), decreasing = TRUE))[1]} + +#build crop_type in prediction files the same way as mvp +make_crop_type = function(dt) { + dt[, crop_type := fifelse( + !is.na(SUBCLASS) & SUBCLASS != "", + paste0(CLASS, "_", SUBCLASS), + as.character(CLASS))] + return(dt)} + +#convert DOY to actual date in prediction year +doy_to_date = function(year, doy) {as.IDate(as.Date(paste0(year, "-01-01")) + as.integer(round(doy)) - 1L)} + +normalize_geoid = function(x) { + x = as.character(x) + x = gsub("\\.0$", "", x) + x = fifelse(is.na(x) | x == "" | x == "NA", NA_character_, x) + x = fifelse(!is.na(x), sprintf("%05s", x), NA_character_) + return(x)} + +mvp[, county_geoid := normalize_geoid(county_geoid)] + +#county + crop_type lookup, collapsing across season +pheno_lookup_county_crop = mvp[ + !is.na(county) & + !is.na(crop_type) & + !is.na(mean_planting_doy) & + !is.na(mean_harvest_doy), + .( + assigned_season = mode_value(season), + + mean_planting_doy = round(weighted.mean(mean_planting_doy, w = n_planting, na.rm = TRUE)), + mean_harvest_doy = round(weighted.mean(mean_harvest_doy, w = n_harvest, na.rm = TRUE)), + mean_peak_doy = round(weighted.mean(mean_peak_doy, w = n_peak, na.rm = TRUE)), + + sd_planting_doy = weighted.mean(sd_planting_doy, w = n_planting, na.rm = TRUE), + sd_harvest_doy = weighted.mean(sd_harvest_doy, w = n_harvest, na.rm = TRUE), + sd_peak_doy = weighted.mean(sd_peak_doy, w = n_peak, na.rm = TRUE), + + n_pheno_rows = sum(n_rows, na.rm = TRUE), + n_planting = sum(n_planting, na.rm = TRUE), + n_harvest = sum(n_harvest, na.rm = TRUE), + n_peak = sum(n_peak, na.rm = TRUE), + + planting_wraparound = any(planting_wraparound, na.rm = TRUE), + harvest_wraparound = any(harvest_wraparound, na.rm = TRUE) + ), + by = .(county, county_geoid, crop_type, landiq_CLASS, landiq_SUBCLASS)] + +#broader fallback: crop_type only, across all counties +pheno_lookup_crop = mvp[ + !is.na(crop_type) & + !is.na(mean_planting_doy) & + !is.na(mean_harvest_doy), + .( + fallback_season = mode_value(season), + + fallback_planting_doy = round(weighted.mean(mean_planting_doy, w = n_planting, na.rm = TRUE)), + fallback_harvest_doy = round(weighted.mean(mean_harvest_doy, w = n_harvest, na.rm = TRUE)), + fallback_peak_doy = round(weighted.mean(mean_peak_doy, w = n_peak, na.rm = TRUE)), + + fallback_sd_planting_doy = weighted.mean(sd_planting_doy, w = n_planting, na.rm = TRUE), + fallback_sd_harvest_doy = weighted.mean(sd_harvest_doy, w = n_harvest, na.rm = TRUE), + fallback_sd_peak_doy = weighted.mean(sd_peak_doy, w = n_peak, na.rm = TRUE), + + fallback_n_pheno_rows = sum(n_rows, na.rm = TRUE), + fallback_planting_wraparound = any(planting_wraparound, na.rm = TRUE), + fallback_harvest_wraparound = any(harvest_wraparound, na.rm = TRUE) + ), + by = .(crop_type, landiq_CLASS, landiq_SUBCLASS)] + +scenario_dirs = c("BAU", "NBS_Targets") + +for (scen in scenario_dirs) { + + input_dir = file.path("county_landiq_predictions", scen) + output_dir = file.path("county_landiq_predictions_with_phenology", scen) + + dir.create(output_dir, recursive = TRUE, showWarnings = FALSE) + + pred_files = list.files(input_dir, pattern = paste0("_predicted_2024_2045", "\\.csv$"), full.names = TRUE) + + message("Processing ", scen, ": ", length(pred_files), " files") + + for (f in pred_files) { + + message(" Reading: ", basename(f)) + + pred = fread(f) + + pred[, `:=`( + parcel_id = as.character(parcel_id), county = as.character(county), county_geoid = normalize_geoid(county_geoid), + CLASS = as.character(CLASS), SUBCLASS = as.character(SUBCLASS), year = as.integer(year))] + + pred[SUBCLASS == "" | SUBCLASS == "NA", SUBCLASS := NA_character_] + pred = make_crop_type(pred) + + #merge county-specific phenology first + pred = merge(pred, pheno_lookup_county_crop, by.x = c("county", "county_geoid", "crop_type", "CLASS", "SUBCLASS"), + by.y = c("county", "county_geoid", "crop_type", "landiq_CLASS", "landiq_SUBCLASS"), all.x = TRUE) + + #merge crop-only fallback + pred = merge(pred, pheno_lookup_crop, by.x = c("crop_type", "CLASS", "SUBCLASS"), + by.y = c("crop_type", "landiq_CLASS", "landiq_SUBCLASS"), all.x = TRUE) + + #use county-specific first, fallback second + pred[is.na(mean_planting_doy), mean_planting_doy := fallback_planting_doy] + pred[is.na(mean_harvest_doy), mean_harvest_doy := fallback_harvest_doy] + pred[is.na(mean_peak_doy), mean_peak_doy := fallback_peak_doy] + + pred[is.na(sd_planting_doy), sd_planting_doy := fallback_sd_planting_doy] + pred[is.na(sd_harvest_doy), sd_harvest_doy := fallback_sd_harvest_doy] + pred[is.na(sd_peak_doy), sd_peak_doy := fallback_sd_peak_doy] + + pred[is.na(n_pheno_rows), n_pheno_rows := fallback_n_pheno_rows] + + pred[is.na(planting_wraparound), planting_wraparound := fallback_planting_wraparound] + pred[is.na(harvest_wraparound), harvest_wraparound := fallback_harvest_wraparound] + + #assign season if predicted season is currently NA + if ("season" %in% names(pred)) { + pred[is.na(season), season := assigned_season] + pred[is.na(season), season := fallback_season] + } else { + pred[, season := assigned_season] + pred[is.na(season), season := fallback_season]} + + #actual predicted dates + pred[, planting_date := doy_to_date(year, mean_planting_doy)] + pred[, peak_date := doy_to_date(year, mean_peak_doy)] + pred[, harvest_date := doy_to_date(year, mean_harvest_doy)] + + #if harvest wraps after New Year, push harvest into next calendar year + pred[ + !is.na(planting_date) & + !is.na(harvest_date) & + harvest_date < planting_date, + harvest_date := doy_to_date(year + 1L, mean_harvest_doy)] + + #same logic for peak if needed + pred[ + !is.na(planting_date) & + !is.na(peak_date) & + peak_date < planting_date & + !is.na(harvest_date) & + harvest_date > planting_date, + peak_date := doy_to_date(year + 1L, mean_peak_doy)] + + pred[, phenology_source := fifelse( + !is.na(assigned_season), + "county_crop_subclass_mean_2018_2023", + fifelse(!is.na(fallback_season), "crop_subclass_mean_2018_2023", NA_character_))] + + #remove fallback helper columns + drop_cols = c("fallback_season", "fallback_planting_doy", "fallback_harvest_doy", "fallback_peak_doy", + "fallback_sd_planting_doy", "fallback_sd_harvest_doy", "fallback_sd_peak_doy", "fallback_n_pheno_rows", + "fallback_planting_wraparound", "fallback_harvest_wraparound") + + drop_cols = intersect(drop_cols, names(pred)) + if (length(drop_cols) > 0) {pred[, (drop_cols) := NULL]} + + #optional: put key columns near front + front_cols = c("parcel_id", "county", "county_geoid", "county_safe", "year", "CLASS", "SUBCLASS", + "CLASS_desc", "SUBCLASS_desc", "PFT", "planting_date", "peak_date", "harvest_date", + "mean_planting_doy", "mean_peak_doy", "mean_harvest_doy", "till_state", "prob_crop_class", "prob_till_state", + "ACRES", "source", "scenario", "phenology_source") + + front_cols = intersect(front_cols, names(pred)) + other_cols = setdiff(names(pred), front_cols) + setcolorder(pred, c(front_cols, other_cols)) + + out_path = file.path(output_dir, basename(f)) + fwrite(pred, out_path)} + + message("Finished phenology assignment for: ", scen)} \ No newline at end of file From a91b89b82698d86fe2ee1f83a59b146dd9904e71 Mon Sep 17 00:00:00 2001 From: Ananya Kulkarni Date: Tue, 7 Jul 2026 19:25:57 -0400 Subject: [PATCH 3/3] updates --- modules/data.remote/inst/predict_and_store.R | 770 ++++++++----------- 1 file changed, 329 insertions(+), 441 deletions(-) diff --git a/modules/data.remote/inst/predict_and_store.R b/modules/data.remote/inst/predict_and_store.R index 13853e5b6c9..e978b2c415f 100644 --- a/modules/data.remote/inst/predict_and_store.R +++ b/modules/data.remote/inst/predict_and_store.R @@ -1,32 +1,27 @@ -#goes through all counties and predicts up until end_year using the optimized matrices & stores county-separated -#LandIQ-style datasets inside bau and nbs target subfolders +#Goes through all counties and predicts up until end_year using the optimized matrices. +#Stores county-separated LandIQ-style datasets inside BAU / NBS_Targets subfolders. setwd("/projectnb/dietzelab/ananyak") library(data.table) ##-------scenario setup------- -run_scenario = "BAU" #"BAU" or "NBS_Targets" +#Crop classes are predicted from county_optimized_matrices/_crop_matrix.csv. +#Tillage states are assigned from county_tillage_targets//crop_by_till_targets_*.csv. -prediction_root_dir = "county_landiq_predictions" -scenario_prediction_dir = file.path(prediction_root_dir, run_scenario) +run_scenarios = c("BAU Targets", "NBS Targets") -dir.create(scenario_prediction_dir, recursive = TRUE, showWarnings = FALSE) +prediction_root_dir = "county_landiq_predictions" +crop_matrix_dir = "county_optimized_matrices" +tillage_target_root = "county_tillage_targets" start_year = 2023L end_year = 2045L -crop_matrix_dir = file.path("county_optimized_matrices", run_scenario) -crop_matrix_pattern = paste0("_crop_matrix_", run_scenario, "\\.csv$") -target_dir = file.path("county_optimized_targets", run_scenario) - -use_till_targets = dir.exists(target_dir) && - length(list.files(target_dir, pattern = paste0("crop_by_till_targets_.*_", run_scenario, "_", end_year, "\\.csv$") - )) > 0 - ##-------helper functions------- safe_county_name = function(x) {gsub("[^A-Za-z0-9_]+", "_", x)} + read_tmat = function(path) { tmat_df = fread(path) @@ -44,37 +39,29 @@ read_tmat = function(path) { return(tmat_final)} + repair_transition_matrix = function(A, matrix_name = "matrix") { A[is.na(A)] = 0 - #clamp tiny numerical optimizer noise if (any(A < 0, na.rm = TRUE)) { - warning( - matrix_name, - " has negative probabilities. Minimum = ", - min(A, na.rm = TRUE), - ". Clamping negatives to 0.") + warning(matrix_name, " has negative probabilities. Minimum = ", + min(A, na.rm = TRUE), ". Clamping negatives to 0.") A[A < 0] = 0} if (any(A > 1, na.rm = TRUE)) { - warning( - matrix_name, - " has probabilities > 1. Maximum = ", - max(A, na.rm = TRUE), - ". Clamping values above 1 to 1.") + warning(matrix_name, " has probabilities > 1. Maximum = ", max(A, na.rm = TRUE), + ". Clamping values above 1 to 1.") A[A > 1] = 1} row_sums = rowSums(A) zero_rows = names(row_sums)[is.na(row_sums) | row_sums == 0] if (length(zero_rows) > 0) { - warning( - matrix_name, " has zero-sum rows. Setting those rows to self-loop: ", + warning(matrix_name, " has zero-sum rows. Setting those rows to self-loop: ", paste(zero_rows, collapse = ", ")) - for (s in zero_rows) { - A[s, ] = 0 + for (s in zero_rows) {A[s, ] = 0 A[s, s] = 1}} row_sums = rowSums(A) @@ -82,18 +69,17 @@ repair_transition_matrix = function(A, matrix_name = "matrix") { return(A)} -load_crop_matrices = function(crop_matrix_dir, crop_matrix_pattern, run_scenario) { + +load_crop_matrices = function(crop_matrix_dir) { - matrix_files = list.files(crop_matrix_dir, pattern = crop_matrix_pattern, full.names = TRUE) + matrix_files = list.files(crop_matrix_dir, pattern = "_crop_matrix\\.csv$", + full.names = TRUE) - if (length(matrix_files) == 0) { - stop("No optimized crop matrix files found in: ", crop_matrix_dir)} + if (length(matrix_files) == 0) {stop("No optimized crop matrix files found in: ", crop_matrix_dir)} transition_mats = list() - for (f in matrix_files) { - - matrix_name = sub(paste0("_crop_matrix_", run_scenario, "\\.csv$"), "", basename(f)) + for (f in matrix_files) {matrix_name = sub("_crop_matrix\\.csv$", "", basename(f)) A = read_tmat(f) A = repair_transition_matrix(A, paste0("optimized crop matrix ", matrix_name)) @@ -102,14 +88,14 @@ load_crop_matrices = function(crop_matrix_dir, crop_matrix_pattern, run_scenario return(transition_mats)} -load_crop_targets = function(target_dir, run_scenario, end_year) { + +load_crop_targets = function(target_dir, scenario_safe, end_year) { - pattern = paste0("crop_targets_.*_", run_scenario, "_", end_year, "\\.csv$") + pattern = paste0("crop_targets_.*_", scenario_safe, "_", end_year, "\\.csv$") files = list.files(target_dir, pattern = pattern, full.names = TRUE) - if (length(files) == 0) { - stop("No crop target files found in: ", target_dir) - } + if (length(files) == 0) {stop("No crop target files found in: ", target_dir, + " with pattern: ", pattern)} out = rbindlist(lapply(files, fread), fill = TRUE) @@ -119,18 +105,16 @@ load_crop_targets = function(target_dir, run_scenario, end_year) { if ("target_acres_used_for_opt" %in% names(out)) { out[, target_acres := as.numeric(target_acres_used_for_opt)] } else { - out[, target_acres := as.numeric(target_acres_raw)] - } + out[, target_acres := as.numeric(target_acres_raw)]} out = out[ - !is.na(county_safe) & - !is.na(crop_state), + !is.na(county_safe) & !is.na(crop_state), .(target_acres = sum(target_acres, na.rm = TRUE)), by = .(county_safe, crop_state) ] - return(out[]) -} + return(out[])} + make_matrix_powers = function(tmat, n_years) { @@ -145,12 +129,11 @@ make_matrix_powers = function(tmat, n_years) { A_power = A_power %*% tmat rownames(A_power) = states colnames(A_power) = states - powers[[i]] = A_power - } + powers[[i]] = A_power} - return(powers) -} + return(powers)} +##-------crop prediction functions------- predict_county_to_targets = function(start_info, tmat, targets_cty, start_year, end_year, state_col = "crop_class") { @@ -165,75 +148,111 @@ predict_county_to_targets = function(start_info, tmat, targets_cty, dt[, ACRES := as.numeric(ACRES)] dt = dt[!is.na(ACRES) & ACRES > 0] - if (nrow(dt) == 0) return(data.table()) + if (nrow(dt) == 0) {return(data.table())} powers = make_matrix_powers(tmat, n_years) A_final = powers[[n_years]] - # full target vector; states not in scenario target file become zero + county_total = sum(dt$ACRES, na.rm = TRUE) + + #Full scenario target vector, used only for diagnostics/checking. target_vec = setNames(rep(0, length(states)), states) - if (nrow(targets_cty) > 0) { - matched = intersect(targets_cty$crop_state, states) - target_vec[matched] = targets_cty[match(matched, crop_state), target_acres] - } - # force target total onto this county's modeled parcel-acre total - county_total = sum(dt$ACRES, na.rm = TRUE) - if (sum(target_vec, na.rm = TRUE) > 0) { - target_vec = target_vec / sum(target_vec, na.rm = TRUE) * county_total + if (nrow(targets_cty) > 0) {matched = intersect(targets_cty$crop_state, states) + target_vec[matched] = targets_cty[match(matched, crop_state), target_acres]} + + if (sum(target_vec, na.rm = TRUE) > 0) {target_vec = target_vec / sum(target_vec, na.rm = TRUE) * county_total} + + #Current county crop acreage vector + x0_dt = dt[, .(acres = sum(ACRES, na.rm = TRUE)), by = start_CLASS] + + X0_vec = setNames(rep(0, length(states)), states) + X0_vec[x0_dt$start_CLASS] = x0_dt$acres + + #Expected 2045 county acres from the optimized scenario matrix + matrix_target_vec = as.numeric(X0_vec %*% A_final) + names(matrix_target_vec) = states + matrix_target_vec[!is.finite(matrix_target_vec)] = 0 + + #Preserve county total + if (sum(matrix_target_vec, na.rm = TRUE) > 0) { + matrix_target_vec = matrix_target_vec / sum(matrix_target_vec, na.rm = TRUE) * county_total } else { - # fallback: if no target file exists, preserve current acreage - cur = dt[, .(acres = sum(ACRES, na.rm = TRUE)), by = start_CLASS] - target_vec[cur$start_CLASS] = cur$acres - } + matrix_target_vec = X0_vec} - # start by keeping all parcels in their observed class - dt[, final_CLASS := start_CLASS] - dt[, locked := FALSE] + #Assign parcels to final classes so county totals follow the optimized matrix projection. + dt[, final_CLASS := NA_character_] + dt[, prob_crop_class_final := NA_real_] - get_current = function(x) { - cur = x[, .(acres = sum(ACRES, na.rm = TRUE)), by = final_CLASS] - out = setNames(rep(0, length(states)), states) - out[cur$final_CLASS] = cur$acres - return(out) - } + remaining_target = copy(matrix_target_vec) - # greedily move parcels from surplus states into deficit states - for (to_state in states[order(-target_vec)]) { - - current_vec = get_current(dt) - need = target_vec[to_state] - current_vec[to_state] + for (to_state in states[order(-matrix_target_vec)]) { - if (!is.finite(need) || need <= 0) next + need = remaining_target[to_state] - surplus_states = names(current_vec)[current_vec > target_vec + 1e-6] - surplus_states = setdiff(surplus_states, to_state) + if (!is.finite(need) || need <= 0) { + next + } - candidates = dt[ - !locked & - final_CLASS %in% surplus_states - ] + candidates = copy(dt[is.na(final_CLASS)]) - if (nrow(candidates) == 0) next + if (nrow(candidates) == 0) { + break + } - candidates[, prob_to := A_final[start_CLASS, to_state]] - candidates[, prob_current := A_final[start_CLASS, final_CLASS]] - candidates[, score := prob_to - prob_current] + candidates[, prob_to := A_final[cbind(start_CLASS, rep(to_state, .N))]] + candidates[is.na(prob_to) | !is.finite(prob_to), prob_to := 0] - setorder(candidates, -score, -prob_to) + #Prefer parcels that the optimized matrix says are likely to become this class. + #Smaller acreage parcels help reduce overshoot. + setorder(candidates, -prob_to, ACRES) candidates[, cum_acres := cumsum(ACRES)] - take_ids = candidates[cum_acres <= need | shift(cum_acres, fill = 0) < need, parcel_id] + + take_ids = candidates[ + cum_acres <= need | shift(cum_acres, fill = 0) < need, + parcel_id] if (length(take_ids) > 0) { + dt[parcel_id %in% take_ids, `:=`( final_CLASS = to_state, - locked = TRUE + prob_crop_class_final = A_final[cbind(start_CLASS, rep(to_state, .N))] )] - } - } - - # build annual time series: keep start class until conversion year + + assigned_acres = dt[parcel_id %in% take_ids, sum(ACRES, na.rm = TRUE)] + remaining_target[to_state] = max(0, remaining_target[to_state] - assigned_acres)}} + + #Any unassigned parcels stay in their original class. + dt[is.na(final_CLASS), `:=`( + final_CLASS = start_CLASS, + prob_crop_class_final = A_final[cbind(start_CLASS, start_CLASS)])] + + #Diagnostic: compare optimized-matrix expected acres vs parcel-assigned acres. + assigned_vec_dt = dt[, .(assigned_acres = sum(ACRES, na.rm = TRUE)), by = final_CLASS] + + assigned_vec = setNames(rep(0, length(states)), states) + assigned_vec[assigned_vec_dt$final_CLASS] = assigned_vec_dt$assigned_acres + + matrix_assignment_check = data.table(CLASS = states, + start_acres = as.numeric(X0_vec[states]), + scenario_target_acres = as.numeric(target_vec[states]), + optimized_matrix_expected_acres = as.numeric(matrix_target_vec[states]), + parcel_assigned_acres = as.numeric(assigned_vec[states])) + + matrix_assignment_check[, matrix_assignment_diff_acres := + parcel_assigned_acres - optimized_matrix_expected_acres] + matrix_assignment_check[, matrix_assignment_abs_diff_acres := + abs(matrix_assignment_diff_acres)] + matrix_assignment_check[, scenario_target_diff_acres := + parcel_assigned_acres - scenario_target_acres] + matrix_assignment_check[, scenario_target_abs_diff_acres := + abs(scenario_target_diff_acres)] + + message("Optimized-matrix parcel assignment check:") + print(matrix_assignment_check[order(-matrix_assignment_abs_diff_acres)]) + + #Build annual time series: keep start class until conversion year. pred_list = vector("list", nrow(dt)) for (i in seq_len(nrow(dt))) { @@ -251,6 +270,8 @@ predict_county_to_targets = function(start_info, tmat, targets_cty, pred_prob[yy] = powers[[yy]][start_state, start_state] } + pred_prob[n_years] = dt$prob_crop_class_final[i] + } else { conversion_idx = which(sapply(seq_len(n_years), function(k) { @@ -268,6 +289,8 @@ predict_county_to_targets = function(start_info, tmat, targets_cty, for (yy in seq_len(n_years)) { pred_prob[yy] = powers[[yy]][start_state, pred_state[yy]] } + + pred_prob[n_years] = dt$prob_crop_class_final[i] } pred_list[[i]] = data.table( @@ -279,13 +302,14 @@ predict_county_to_targets = function(start_info, tmat, targets_cty, } out = rbindlist(pred_list, fill = TRUE) - return(out[]) -} + return(out[])} + predict_grouped_markov_to_targets = function(year_states, transition_mats, crop_targets_2045, group_col, - start_year, end_year, + start_year, + end_year, state_col = "crop_class") { dt = copy(year_states) @@ -293,9 +317,11 @@ predict_grouped_markov_to_targets = function(year_states, transition_mats, groups = intersect(unique(na.omit(dt[[group_col]])), names(transition_mats)) + if (length(groups) == 0) {stop("No overlapping groups between all_data and transition matrices.")} + for (g in groups) { - message("Predicting crop class with target allocation for county: ", g) + message("Predicting crop class with optimized scenario matrix for county: ", g) dt_g = dt[get(group_col) == g] tmat_g = transition_mats[[g]] @@ -303,31 +329,23 @@ predict_grouped_markov_to_targets = function(year_states, transition_mats, start_info = dt_g[ year <= start_year, .SD[which.max(year)], - by = parcel_id - ] + by = parcel_id] targets_g = crop_targets_2045[county_safe == g] - preds_g = predict_county_to_targets( - start_info = start_info, - tmat = tmat_g, - targets_cty = targets_g, - start_year = start_year, - end_year = end_year, - state_col = state_col - ) + preds_g = predict_county_to_targets(start_info = start_info, tmat = tmat_g, + targets_cty = targets_g, start_year = start_year, + end_year = end_year, state_col = state_col) - if (nrow(preds_g) == 0) next + if (nrow(preds_g) == 0) { + next + } preds_g[, (group_col) := g] all_preds[[g]] = preds_g } - return(rbindlist(all_preds, fill = TRUE)) -} - - - + return(rbindlist(all_preds, fill = TRUE))} ##--------load & clean all_data--------- all_data = fread("all_data.csv") @@ -348,11 +366,13 @@ if (length(missing_cols) > 0) { stop("all_data is missing required columns: ", paste(missing_cols, collapse = ", "))} all_data[, `:=`( - parcel_id = as.character(parcel_id), year = as.integer(year), county = as.character(county), + parcel_id = as.character(parcel_id), year = as.integer(year), county = as.character(county), crop_class = trimws(as.character(crop_class)), ACRES = as.numeric(ACRES))] -if (!"till_state" %in% names(all_data)) {all_data[, till_state := NA_character_]} else { - all_data[, till_state := trimws(as.character(till_state))]} +if (!"till_state" %in% names(all_data)) {all_data[, till_state := NA_character_] +} else { + all_data[, till_state := trimws(as.character(till_state))] +} if (!"county_geoid" %in% names(all_data)) {all_data[, county_geoid := NA_character_]} @@ -361,23 +381,20 @@ if (!"season" %in% names(all_data)) {all_data[, season := NA_integer_]} all_data[, county_safe := safe_county_name(county)] ##--------load lookup table--------- - lookup = fread("/projectnb/dietzelab/ccmmf/management/LandIQ_cropCode_lookup_table.csv") lookup[, `:=`( CLASS = as.character(CLASS), SUBCLASS = as.character(SUBCLASS))] -lookup_subclass = unique(lookup[, .( - CLASS, SUBCLASS, CLASS_desc, SUBCLASS_desc, PFT)], - by = c("CLASS", "SUBCLASS")) +lookup_subclass = unique(lookup[, .(CLASS, SUBCLASS, CLASS_desc, SUBCLASS_desc, PFT)], + by = c("CLASS", "SUBCLASS")) -##--------load historical LandIQ source with subclass--------- +##--------load historical LandIQ source with subclass--------- with_subclass = fread("crops_full_counties.csv") if ("V1" %in% names(with_subclass)) {with_subclass[, V1 := NULL]} -#standardize class, sublcass, season, and county name columns if (!"CLASS" %in% names(with_subclass) && "crop_class" %in% names(with_subclass)) { with_subclass[, CLASS := as.character(crop_class)]} @@ -394,11 +411,11 @@ missing_subclass_cols = setdiff(required_subclass_cols, names(with_subclass)) if (length(missing_subclass_cols) > 0) { stop("crops_full_counties.csv is missing required columns: ", - paste(missing_subclass_cols, collapse = ", "), "\nAvailable columns: ", - paste(names(with_subclass), collapse = ", "))} + paste(missing_subclass_cols, collapse = ", "), + "\nAvailable columns: ", paste(names(with_subclass), collapse = ", "))} with_subclass[, `:=`( - parcel_id = as.character(parcel_id), year = as.integer(year), + parcel_id = as.character(parcel_id), year = as.integer(year), county_safe = as.character(county_safe), CLASS = trimws(as.character(CLASS)), SUBCLASS = trimws(as.character(SUBCLASS)), season = as.integer(season))] @@ -411,6 +428,7 @@ message("Loaded subclass source with ", nrow(with_subclass[!is.na(SUBCLASS)]), if (!"CLASS" %in% names(all_data)) {all_data[, CLASS := as.character(crop_class)]} + assign_predicted_subclass = function(future_landiq, subclass_obs, lookup_subclass, crop_col = "CLASS", group_col = "county_safe", start_year = 2023L) { @@ -422,8 +440,7 @@ assign_predicted_subclass = function(future_landiq, subclass_obs, lookup_subclas obs[, `:=`( parcel_id = as.character(parcel_id), year = as.integer(year), - county_safe = as.character(county_safe), CLASS = as.character(CLASS), - SUBCLASS = as.character(SUBCLASS))] + county_safe = as.character(county_safe),CLASS = as.character(CLASS), SUBCLASS = as.character(SUBCLASS))] if (!"season" %in% names(obs)) { obs[, season := 0L] @@ -432,44 +449,35 @@ assign_predicted_subclass = function(future_landiq, subclass_obs, lookup_subclas obs[is.na(season), season := 0L] } - #only use historical data up to the prediction start year obs = obs[year <= start_year] dt[, parcel_id := as.character(parcel_id)] dt[, CLASS := as.character(get(crop_col))] dt[, (group_col) := as.character(get(group_col))] - #remove old lookup/subclass fields if function is accidentally rerun old_lookup_cols = intersect(c("SUBCLASS", "CLASS_desc", "SUBCLASS_desc", "PFT"), names(dt)) if (length(old_lookup_cols) > 0) {dt[, (old_lookup_cols) := NULL]} - #global subclass probabilities global_probs = obs[ !is.na(CLASS) & !is.na(SUBCLASS), .N, by = .(CLASS, SUBCLASS)] if (nrow(global_probs) > 0) {global_probs[, prob := N / sum(N), by = CLASS]} - #county-specific subclass probabilities group_probs = obs[ !is.na(CLASS) & !is.na(SUBCLASS), .N, by = .(county_safe, CLASS, SUBCLASS)] if (nrow(group_probs) > 0) {group_probs[, prob := N / sum(N), by = .(county_safe, CLASS)]} - #lookup fallback if a predicted class has no observed subclass distribution - lookup_probs = unique(lookup_subclass[ - !is.na(CLASS) & !is.na(SUBCLASS), - .(CLASS, SUBCLASS)]) + lookup_probs = unique(lookup_subclass[!is.na(CLASS) & !is.na(SUBCLASS), .(CLASS, SUBCLASS)]) if (nrow(lookup_probs) > 0) {lookup_probs[, N := 1L] lookup_probs[, prob := 1 / .N, by = CLASS]} - # Last observed class/subclass per parcel + county from subclass source - last_obs_source = obs[ - !is.na(CLASS) & !is.na(SUBCLASS)] + last_obs_source = obs[!is.na(CLASS) & !is.na(SUBCLASS)] last_obs = last_obs_source[ order(year, season), @@ -481,40 +489,34 @@ assign_predicted_subclass = function(future_landiq, subclass_obs, lookup_subclas last_obs = last_obs[ , - .( - parcel_id, county_safe, last_CLASS = CLASS, last_SUBCLASS = SUBCLASS)] + .(parcel_id, county_safe, last_CLASS = CLASS, last_SUBCLASS = SUBCLASS) + ] dt = merge(dt, last_obs, by = c("parcel_id", "county_safe"), all.x = TRUE) setorder(dt, parcel_id, year) - #find runs of same predicted CLASS within each parcel dt[, prev_CLASS := shift(CLASS), by = .(parcel_id, county_safe)] dt[is.na(prev_CLASS), prev_CLASS := last_CLASS] dt[, new_run := fifelse(is.na(CLASS), FALSE, is.na(prev_CLASS) | CLASS != prev_CLASS)] - dt[, run_id := cumsum(new_run), by = .(parcel_id, county_safe)] dt[, SUBCLASS := NA_character_] - #if predicted class matches last observed class, carry forward real historical subclass - dt[run_id == 0 & + dt[ + run_id == 0 & !is.na(CLASS) & !is.na(last_CLASS) & CLASS == last_CLASS & !is.na(last_SUBCLASS), SUBCLASS := last_SUBCLASS] - #for remaining class-runs, draw subclass from county/class distribution - run_table = unique(dt[ - !is.na(CLASS) & is.na(SUBCLASS), - .(parcel_id, county_safe, run_id, CLASS)]) + run_table = unique(dt[!is.na(CLASS) & is.na(SUBCLASS), .(parcel_id, county_safe, run_id, CLASS)]) if (nrow(run_table) > 0) { run_table[, drawn_SUBCLASS := NA_character_] - draw_groups = unique(run_table[, .(county_safe, CLASS)]) for (ii in seq_len(nrow(draw_groups))) { @@ -542,8 +544,7 @@ assign_predicted_subclass = function(future_landiq, subclass_obs, lookup_subclas !is.na(SUBCLASS) & !is.na(prob) & is.finite(prob) & - prob > 0 - ] + prob > 0] if (nrow(choices) > 0) { choices[, prob := prob / sum(prob)] @@ -561,176 +562,36 @@ assign_predicted_subclass = function(future_landiq, subclass_obs, lookup_subclas by = c("parcel_id", "county_safe", "run_id"), all.x = TRUE) dt[is.na(SUBCLASS), SUBCLASS := drawn_SUBCLASS] - dt[, drawn_SUBCLASS := NULL]} + dt[, drawn_SUBCLASS := NULL] + } helper_cols = intersect(c("last_CLASS", "last_SUBCLASS", "prev_CLASS", "new_run", "run_id"), - names(dt)) + names(dt)) dt[, (helper_cols) := NULL] - if (!"SUBCLASS" %in% names(dt)) {dt[, SUBCLASS := NA_character_]} + if (!"SUBCLASS" %in% names(dt)) { + dt[, SUBCLASS := NA_character_] + } - if (!"CLASS" %in% names(dt)) {dt[, CLASS := NA_character_]} + if (!"CLASS" %in% names(dt)) { + dt[, CLASS := NA_character_] + } - #add lookup descriptions/PFT dt = merge(dt, lookup_subclass, by = c("CLASS", "SUBCLASS"), all.x = TRUE) setorder(dt, orig_order) dt[, orig_order := NULL] return(dt[])} - -##--------crop prediction functions--------- -predict_markov = function(start_info, tmat, end_year, id_col = "parcel_id", time_col = "year", state_col = "crop_class") { - - dt = copy(start_info) - - if (!(id_col %in% names(dt))) {stop("id_col not found in start_info: ", id_col)} - if (!(time_col %in% names(dt))) {stop("time_col not found in start_info: ", time_col)} - if (!(state_col %in% names(dt))) {stop("state_col not found in start_info: ", state_col)} - - if (id_col != "parcel_id") {setnames(dt, id_col, "parcel_id")} - if (time_col != "year") {setnames(dt, time_col, "year")} - - states = rownames(tmat) - - dt = dt[get(state_col) %in% states] - - if (nrow(dt) == 0) {return(data.table())} - - #clean matrix - tmat[is.na(tmat)] = 0 - tmat[tmat < 0] = 0 - - row_sums = rowSums(tmat) - for (s in states[row_sums <= 0 | is.na(row_sums)]) { - tmat[s, ] = 0 - tmat[s, s] = 1} - - tmat = sweep(tmat, 1, rowSums(tmat), "/") - - preds = dt[, { - - start_state = get(state_col) - years = seq(year + 1L, end_year) - n_years = length(years) - - if (is.na(start_state) || !(start_state %in% states) || n_years == 0) { - .(year = integer(), CLASS = character(), prob_crop_class = numeric()) - } else { - - # Track cumulative transition probabilities from start year to each future year - A_power = diag(length(states)) - rownames(A_power) = states - colnames(A_power) = states - - yearly_probs = vector("list", n_years) - - for (i in seq_len(n_years)) { - A_power = A_power %*% tmat - rownames(A_power) = states - colnames(A_power) = states - - p_i = as.numeric(A_power[start_state, ]) - names(p_i) = states - p_i[is.na(p_i)] = 0 - - if (sum(p_i) > 0) { - p_i = p_i / sum(p_i) - } - - yearly_probs[[i]] = p_i - } - - final_probs = yearly_probs[[n_years]] - final_state = names(final_probs)[which.max(final_probs)] - final_prob = final_probs[final_state] - - # Default: parcel stays constant - pred_state = rep(start_state, n_years) - pred_prob = sapply(yearly_probs, function(p) p[start_state]) - - # If final most-likely state is different, allow ONE conversion only - if (!is.na(final_state) && final_state != start_state) { - - # convert only when final state becomes at least as likely as start state - conversion_idx = which(sapply(yearly_probs, function(p) { - p[final_state] >= p[start_state] - })) - - if (length(conversion_idx) == 0) { - conversion_idx = n_years - } else { - conversion_idx = conversion_idx[1] - } - - pred_state[conversion_idx:n_years] = final_state - pred_prob[conversion_idx:n_years] = sapply( - yearly_probs[conversion_idx:n_years], - function(p) p[final_state] - ) - } - - .( - year = years, - CLASS = pred_state, - prob_crop_class = as.numeric(pred_prob) - ) - } - - }, by = parcel_id] - - return(preds) -} - -predict_grouped_markov = function(year_states, transition_mats, group_col, start_year, end_year, id_col = "parcel_id", - time_col = "year", state_col = "crop_class") { - - dt = copy(year_states) - - if (!(group_col %in% names(dt))) {stop("group_col not found in year_states: ", group_col)} - - if (id_col != "parcel_id") {setnames(dt, id_col, "parcel_id")} - - if (time_col != "year") {setnames(dt, time_col, "year")} - - all_preds = list() - - groups = intersect(unique(na.omit(dt[[group_col]])), names(transition_mats)) - - if (length(groups) == 0) {stop("No overlapping groups between all_data and transition matrices.")} - - for (g in groups) { - - message("Predicting crop class for county: ", g) - - dt_g = dt[get(group_col) == g] - tmat_g = transition_mats[[g]] - - # latest observed state for each parcel up to start_year - start_info = dt_g[year <= start_year, - .SD[which.max(year)], by = parcel_id] - - preds_g = predict_markov(start_info = start_info, tmat = tmat_g, end_year = end_year, state_col = state_col) - - if (nrow(preds_g) == 0) { - next} - - preds_g[, (group_col) := g] - - all_preds[[g]] = preds_g} - - return(rbindlist(all_preds, fill = TRUE))} - ##--------till target functions--------- - -load_crop_by_till_targets = function(target_dir, run_scenario, end_year) { - - pattern = paste0("crop_by_till_targets_.*_", run_scenario, "_", end_year, "\\.csv$") +load_crop_by_till_targets = function(target_dir, scenario_safe, end_year) { + pattern = paste0("crop_by_till_targets_.*_", scenario_safe, "_", end_year, "\\.csv$") files = list.files(target_dir, pattern = pattern, full.names = TRUE) - if (length(files) == 0) {stop("No crop-by-till target files found in: ", target_dir)} + if (length(files) == 0) {stop("No crop-by-till target files found in: ", target_dir, + " with pattern: ", pattern)} out = rbindlist(lapply(files, fread), fill = TRUE) @@ -739,38 +600,44 @@ load_crop_by_till_targets = function(target_dir, run_scenario, end_year) { out[, till_state := as.character(till_state)] out[, target_acres_raw := as.numeric(target_acres_raw)] - return(out)} + return(out[])} build_annual_till_targets = function(all_data, till_targets_2045, start_year, end_year) { - latest_obs = all_data[year <= start_year & !is.na(crop_class) & !is.na(till_state), - .SD[which.max(year)], by = parcel_id] + latest_obs = all_data[ + year <= start_year & !is.na(crop_class) & !is.na(till_state), + .SD[which.max(year)], + by = parcel_id] baseline_till = latest_obs[ , .(baseline_acres = sum(ACRES, na.rm = TRUE)), by = .(county_safe, crop_state = crop_class, till_state)] - if (nrow(baseline_till) > 0) {baseline_till[, baseline_share := baseline_acres / sum(baseline_acres), - by = .(county_safe, crop_state)] + if (nrow(baseline_till) > 0) { + baseline_till[, baseline_share := baseline_acres / sum(baseline_acres), + by = .(county_safe, crop_state)] } else { - baseline_till[, baseline_share := numeric()]} + baseline_till[, baseline_share := numeric()] + } target_till = copy(till_targets_2045) - - target_till[, target_share := target_acres_raw / sum(target_acres_raw), by = .(county_safe, crop_state)] + target_till[, target_share := target_acres_raw / sum(target_acres_raw), + by = .(county_safe, crop_state)] all_counties = unique(target_till$county_safe) all_crops = unique(target_till$crop_state) all_till_states = unique(c(baseline_till$till_state, target_till$till_state)) - annual_till_targets = CJ(county_safe = all_counties, crop_state = all_crops, till_state = all_till_states, - year = seq(start_year + 1L, end_year)) + annual_till_targets = CJ(county_safe = all_counties, crop_state = all_crops, + till_state = all_till_states, year = seq(start_year + 1L, end_year)) - annual_till_targets = merge(annual_till_targets, baseline_till[, .(county_safe, crop_state, till_state, baseline_share)], - by = c("county_safe", "crop_state", "till_state"), all.x = TRUE) + annual_till_targets = merge(annual_till_targets, + baseline_till[, .(county_safe, crop_state, till_state, baseline_share)], + by = c("county_safe", "crop_state", "till_state"), all.x = TRUE) - annual_till_targets = merge(annual_till_targets, target_till[, .(county_safe, crop_state, till_state, target_share)], - by = c("county_safe", "crop_state", "till_state"), all.x = TRUE) + annual_till_targets = merge(annual_till_targets, + target_till[, .(county_safe, crop_state, till_state, target_share)], + by = c("county_safe", "crop_state", "till_state"), all.x = TRUE) annual_till_targets[is.na(baseline_share), baseline_share := 0] annual_till_targets[is.na(target_share), target_share := 0] @@ -783,14 +650,14 @@ build_annual_till_targets = function(all_data, till_targets_2045, start_year, en annual_till_targets[ , share_sum := sum(till_share, na.rm = TRUE), by = .(county_safe, crop_state, year)] - + annual_till_targets[share_sum > 0, till_share := till_share / share_sum] - - annual_till_targets[share_sum <= 0, till_share := 1 / .N, by = .(county_safe, crop_state, year)] + annual_till_targets[share_sum <= 0, till_share := 1 / .N, + by = .(county_safe, crop_state, year)] annual_till_targets[, share_sum := NULL] - return(annual_till_targets)} + return(annual_till_targets[])} assign_till_by_targets = function(future_landiq, annual_till_targets) { @@ -853,14 +720,12 @@ assign_till_by_targets = function(future_landiq, annual_till_targets) { } else { - #randomize row order so the same parcels are not always assigned first temp[, rand_order := runif(.N)] setorder(temp, rand_order) targets[, target_acres := till_share * group_total_acres] targets[, upper_acres := cumsum(target_acres)] - #use cumulative parcel acreage midpoint for assignment temp[, cum_mid_acres := cumsum(ACRES) - ACRES / 2] idx = findInterval(temp$cum_mid_acres, targets$upper_acres) + 1L @@ -883,131 +748,154 @@ assign_till_by_targets = function(future_landiq, annual_till_targets) { out[, row_id := NULL] return(out[])} -##--------run predictions--------- -crop_mats = load_crop_matrices(crop_matrix_dir = crop_matrix_dir, crop_matrix_pattern = crop_matrix_pattern, - run_scenario = run_scenario) - -missing_mats = setdiff(unique(all_data$county_safe), names(crop_mats)) - -if (length(missing_mats) > 0) { - message("Counties in all_data without crop matrices: ", length(missing_mats))} -crop_targets_2045 = load_crop_targets( - target_dir = target_dir, - run_scenario = run_scenario, - end_year = end_year -) -future_crop = predict_grouped_markov_to_targets( - year_states = all_data, - transition_mats = crop_mats, - crop_targets_2045 = crop_targets_2045, - group_col = "county_safe", - start_year = start_year, - end_year = end_year, - state_col = "crop_class" -) - -if (nrow(future_crop) == 0) {stop("No future crop predictions were generated.")} - -##--------add parcel metadata--------- -parcel_meta = all_data[year <= start_year, - .SD[which.max(year)], by = parcel_id -][ - , - .( - parcel_id, county, county_geoid, county_safe, ACRES, till_state_2023 = till_state)] +##--------run predictions--------- +#Load shared optimized crop matrices once. These are NOT scenario-nested anymore. +crop_mats = load_crop_matrices(crop_matrix_dir = crop_matrix_dir) -future_landiq = merge(future_crop, parcel_meta, by = c("parcel_id", "county_safe"), all.x = TRUE) +missing_mats = setdiff(unique(all_data$county_safe), names(crop_mats)) +if (length(missing_mats) > 0) {message("Counties in all_data without crop matrices: ", length(missing_mats))} -##--------assign till states--------- -if (use_till_targets) { +run_one_prediction_scenario = function(run_scenario) { + + scenario_safe = safe_county_name(run_scenario) + scenario_prediction_dir = file.path(prediction_root_dir, scenario_safe) + target_dir = file.path(tillage_target_root, scenario_safe) + + dir.create(scenario_prediction_dir, recursive = TRUE, showWarnings = FALSE) + + message("\n==============================") + message("Running prediction scenario: ", run_scenario) + message("Scenario-safe folder name: ", scenario_safe) + message("Reading tillage/crop targets from: ", target_dir) + message("Writing predictions to: ", scenario_prediction_dir) + message("==============================") + + if (!dir.exists(target_dir)) { + stop("Scenario target folder not found: ", target_dir)} + + ##The actual crop-class prediction is driven by the shared optimized crop matrix. + crop_targets_2045 = load_crop_targets(target_dir = target_dir, + scenario_safe = scenario_safe, end_year = end_year) + + future_crop = predict_grouped_markov_to_targets( + year_states = all_data, transition_mats = crop_mats, + crop_targets_2045 = crop_targets_2045, group_col = "county_safe", + start_year = start_year, end_year = end_year, state_col = "crop_class") + + if (nrow(future_crop) == 0) {stop("No future crop predictions were generated for scenario: ", run_scenario)} + + ##--------add parcel metadata--------- + parcel_meta = all_data[ + year <= start_year, + .SD[which.max(year)], + by = parcel_id + ][ + , + .(parcel_id, county, county_geoid, county_safe, ACRES, till_state_2023 = till_state) + ] + + future_landiq = merge( + future_crop, + parcel_meta, + by = c("parcel_id", "county_safe"), + all.x = TRUE + ) + ##--------assign till states from scenario-specific tillage targets--------- message("Assigning till states using ", run_scenario, " crop-by-till targets.") - till_targets_2045 = load_crop_by_till_targets(target_dir = target_dir, run_scenario = run_scenario, end_year = end_year) + till_targets_2045 = load_crop_by_till_targets( + target_dir = target_dir, + scenario_safe = scenario_safe, + end_year = end_year + ) + + annual_till_targets = build_annual_till_targets( + all_data = all_data, + till_targets_2045 = till_targets_2045, + start_year = start_year, + end_year = end_year + ) + + future_landiq = assign_till_by_targets( + future_landiq = future_landiq, + annual_till_targets = annual_till_targets + ) + + ##--------final formatting--------- + future_landiq[, season := NA_integer_] + future_landiq[, source := "predicted"] + future_landiq[, scenario := run_scenario] + future_landiq[, scenario_safe := scenario_safe] + + future_landiq = future_landiq[year >= start_year + 1L & year <= end_year] + + future_landiq = assign_predicted_subclass( + future_landiq = future_landiq, + subclass_obs = with_subclass, + lookup_subclass = lookup_subclass, + crop_col = "CLASS", + group_col = "county_safe", + start_year = start_year + ) + + if ("till_state_2023" %in% names(future_landiq)) { + future_landiq[, till_state_2023 := NULL] + } - annual_till_targets = build_annual_till_targets(all_data = all_data, till_targets_2045 = till_targets_2045, - start_year = start_year, end_year = end_year) + optional_cols = c("SUBCLASS", "CLASS_desc", "SUBCLASS_desc", "PFT") - future_landiq = assign_till_by_targets(future_landiq = future_landiq, annual_till_targets = annual_till_targets) + for (cc in optional_cols) { + if (!(cc %in% names(future_landiq))) { + future_landiq[, (cc) := NA_character_] + } + } -} else { + final_cols = c("parcel_id", "county", "county_geoid", "county_safe", "year", "season", + "CLASS", "SUBCLASS", "CLASS_desc", "SUBCLASS_desc", "PFT", + "till_state", "prob_crop_class", "prob_till_state", "ACRES", "source", + "scenario") - message(run_scenario, " has no crop-by-till target files. Carrying forward 2023 observed till_state.") + future_landiq = future_landiq[, ..final_cols] + setorder(future_landiq, county_safe, parcel_id, year) - future_landiq[, till_state := till_state_2023] - future_landiq[, prob_till_state := 1]} - -##--------final formatting--------- -future_landiq[, season := NA_integer_] -future_landiq[, source := "predicted"] -future_landiq[, scenario := run_scenario] - -future_landiq = future_landiq[year >= start_year + 1L & year <= end_year] - -#add subclass + lookup fields -future_landiq = assign_predicted_subclass( future_landiq = future_landiq, - subclass_obs = with_subclass, - lookup_subclass = lookup_subclass, - crop_col = "CLASS", group_col = "county_safe", - start_year = start_year) - -#clean optional helper columns -if ("till_state_2023" %in% names(future_landiq)) {future_landiq[, till_state_2023 := NULL]} - -#make sure optional columns exist before selecting -optional_cols = c("SUBCLASS", "CLASS_desc", "SUBCLASS_desc", "PFT") - -for (cc in optional_cols) {if (!(cc %in% names(future_landiq))) { - future_landiq[, (cc) := NA_character_]}} - -final_cols = c("parcel_id", "county", "county_geoid", "county_safe", "year", "season", "CLASS", "SUBCLASS", "CLASS_desc", - "SUBCLASS_desc", "PFT", "till_state", "prob_crop_class", "prob_till_state", "ACRES", "source", "scenario") - -future_landiq = future_landiq[, ..final_cols] - -setorder(future_landiq, county_safe, parcel_id, year) - - -pred_2045_check = future_landiq[ - year == end_year, - .(predicted_acres = sum(ACRES, na.rm = TRUE)), - by = .(county_safe, CLASS) -] - -target_2045_check = copy(crop_targets_2045) -setnames(target_2045_check, "crop_state", "CLASS") - -check_2045 = merge( - pred_2045_check, - target_2045_check, - by = c("county_safe", "CLASS"), - all = TRUE -) - -check_2045[is.na(predicted_acres), predicted_acres := 0] -check_2045[is.na(target_acres), target_acres := 0] -check_2045[, diff_acres := predicted_acres - target_acres] -check_2045[, abs_diff_acres := abs(diff_acres)] - -check_path = file.path( - scenario_prediction_dir, - paste0("prediction_target_check_", run_scenario, "_", end_year, ".csv") -) - -fwrite(check_2045, check_path) - -message("Wrote 2045 target check to: ", check_path) - - - - -##--------write county-separated outputs--------- -for (cty in unique(na.omit(future_landiq$county_safe))) { + ##--------write county-separated outputs--------- + for (cty in unique(na.omit(future_landiq$county_safe))) { + + out_path = file.path(scenario_prediction_dir, paste0(cty, "_predicted_2024_", end_year, ".csv")) + + fwrite(future_landiq[county_safe == cty], out_path)} - out_path = file.path(scenario_prediction_dir, paste0(cty, "_predicted_2024_", end_year, ".csv")) + scenario_manifest = future_landiq[, .( + n_rows = .N, + n_parcels = uniqueN(parcel_id), + total_acres = sum(ACRES, na.rm = TRUE) + ), by = .(scenario, county_safe)] - fwrite(future_landiq[county_safe == cty], out_path)} + manifest_path = file.path(scenario_prediction_dir, paste0("prediction_manifest_", scenario_safe, ".csv")) + fwrite(scenario_manifest, manifest_path) + + message("Finished writing predictions to: ", scenario_prediction_dir) + return(scenario_manifest[])} + +all_prediction_manifests = rbindlist(lapply(run_scenarios, function(scen) { + tryCatch( + run_one_prediction_scenario(scen), + error = function(e) { + data.table( + scenario = scen, + scenario_safe = safe_county_name(scen), + run_status = "error", + error_message = conditionMessage(e) + ) + } + ) +}), fill = TRUE) + +all_prediction_manifest_path = file.path(prediction_root_dir, "all_prediction_manifests.csv") +dir.create(prediction_root_dir, recursive = TRUE, showWarnings = FALSE) +fwrite(all_prediction_manifests, all_prediction_manifest_path) -message("Finished writing predictions to: ", scenario_prediction_dir) \ No newline at end of file +print(all_prediction_manifests)