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 diff --git a/modules/data.remote/inst/predict_and_store.R b/modules/data.remote/inst/predict_and_store.R index 5bb6f6750b8..e978b2c415f 100644 --- a/modules/data.remote/inst/predict_and_store.R +++ b/modules/data.remote/inst/predict_and_store.R @@ -1,299 +1,901 @@ -##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 / NBS_Targets 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") +##-------scenario setup------- +#Crop classes are predicted from county_optimized_matrices/_crop_matrix.csv. +#Tillage states are assigned from county_tillage_targets//crop_by_till_targets_*.csv. + +run_scenarios = c("BAU Targets", "NBS Targets") -row_id_col = colnames(tmat_df)[1] -states = colnames(tmat_df)[-1] +prediction_root_dir = "county_landiq_predictions" +crop_matrix_dir = "county_optimized_matrices" +tillage_target_root = "county_tillage_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" +start_year = 2023L +end_year = 2045L -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] +##-------helper functions------- +safe_county_name = function(x) {gsub("[^A-Za-z0-9_]+", "_", x)} -#clean classes and set in same order as transition matrix -start_info[, dominant_crop := trimws(as.character(dominant_crop))] -states = rownames(tmat_final) -##specify how long into the future the user wants to predict to -end_year = 2050 +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)} -##--------design point predictions--------- -design_points = fread('/projectnb/dietzelab/ccmmf/management/design_points_landiq_2018-2023.csv') -design_points[, parcel_id := as.character(parcel_id)] -start_info[, parcel_id := as.character(parcel_id)] -year_states[, parcel_id := as.character(parcel_id)] +repair_transition_matrix = function(A, matrix_name = "matrix") { + + A[is.na(A)] = 0 + + 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} + + 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} + + 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: ", + paste(zero_rows, collapse = ", ")) + + for (s in zero_rows) {A[s, ] = 0 + A[s, s] = 1}} + + row_sums = rowSums(A) + A = sweep(A, 1, row_sums, "/") + + return(A)} -#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] -##function that pulls transition matrix probabilities -all_preds = start_info[, { +load_crop_matrices = function(crop_matrix_dir) { + + matrix_files = list.files(crop_matrix_dir, pattern = "_crop_matrix\\.csv$", + full.names = TRUE) - current_state = dominant_crop - years = seq(year + 1, end_year) + if (length(matrix_files) == 0) {stop("No optimized crop matrix files found in: ", crop_matrix_dir)} - preds = character(length(years)) - probs = numeric(length(years)) + transition_mats = list() - for (i in seq_along(years)) { - p = tmat_year[current_state, ] + for (f in matrix_files) {matrix_name = sub("_crop_matrix\\.csv$", "", basename(f)) - if (sum(p) == 0) { - preds[i] = NA_character_ - probs[i] = NA_real_ + A = read_tmat(f) + A = repair_transition_matrix(A, paste0("optimized crop matrix ", matrix_name)) + + transition_mats[[matrix_name]] = A} + + return(transition_mats)} + + +load_crop_targets = function(target_dir, scenario_safe, end_year) { + + 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, + " with pattern: ", pattern)} + + 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[])} + + +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)} + +##-------crop prediction functions------- +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]] + + 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]} + + 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 { + matrix_target_vec = X0_vec} + + #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_] + + remaining_target = copy(matrix_target_vec) + + for (to_state in states[order(-matrix_target_vec)]) { + + need = remaining_target[to_state] + + if (!is.finite(need) || need <= 0) { next } - idx = which.max(p) - next_state = states[idx] + candidates = copy(dt[is.na(final_CLASS)]) + + if (nrow(candidates) == 0) { + break + } - preds[i] = next_state - probs[i] = p[idx] + candidates[, prob_to := A_final[cbind(start_CLASS, rep(to_state, .N))]] + candidates[is.na(prob_to) | !is.finite(prob_to), prob_to := 0] - current_state = next_state - } + #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] + + if (length(take_ids) > 0) { + + dt[parcel_id %in% take_ids, `:=`( + final_CLASS = to_state, + prob_crop_class_final = A_final[cbind(start_CLASS, rep(to_state, .N))] + )] + + 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)}} - .( - year = years, - pred_class = preds, - pred_prob = probs - ) + #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)])] -}, 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 - ) -] - -#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) -] + #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))) { + + 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] + } + + pred_prob[n_years] = dt$prob_crop_class_final[i] + + } 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_prob[n_years] = dt$prob_crop_class_final[i] + } + + 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[])} -subclass_probs[, prob := N / sum(N), by = CLASS] -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") { - choices = subclass_probs[CLASS == class_name] + dt = copy(year_states) + all_preds = list() - if (nrow(choices) == 0) { - return(NA_character_) + 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 optimized scenario matrix 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))} +##--------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))] } -#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 -) +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]} + +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))] + +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.") -setorder(preds_subclass, parcel_id, year) +if (!"CLASS" %in% names(all_data)) {all_data[, CLASS := as.character(crop_class)]} -preds_subclass[, SUBCLASS := NA_character_] -#assigning subclass year by year -for (p in unique(preds_subclass$parcel_id)) { +assign_predicted_subclass = function(future_landiq, subclass_obs, lookup_subclass, + crop_col = "CLASS", group_col = "county_safe", + start_year = 2023L) { - idxs = which(preds_subclass$parcel_id == p) + dt = copy(future_landiq) + obs = copy(subclass_obs) - prev_class = preds_subclass$last_CLASS[idxs[1]] - prev_subclass = preds_subclass$last_SUBCLASS[idxs[1]] + dt[, orig_order := .I] - for (i in idxs) { + 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] + } + + 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))] + + 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_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]} + + 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_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_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) + + 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_] + + dt[ + run_id == 0 & + !is.na(CLASS) & + !is.na(last_CLASS) & + CLASS == last_CLASS & + !is.na(last_SUBCLASS), + SUBCLASS := last_SUBCLASS] + + run_table = unique(dt[!is.na(CLASS) & is.na(SUBCLASS), .(parcel_id, county_safe, run_id, CLASS)]) + + if (nrow(run_table) > 0) { - current_class = preds_subclass$CLASS[i] + run_table[, drawn_SUBCLASS := NA_character_] + 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))) { - } else if (!is.na(prev_class) && current_class == prev_class) { + cty = draw_groups$county_safe[ii] + cls = draw_groups$CLASS[ii] - #if a crop class does not transition, keep the subclass its currently in - preds_subclass$SUBCLASS[i] = prev_subclass + idx = which(run_table$county_safe == cty & run_table$CLASS == cls) - } else { + choices = group_probs[county_safe == cty & CLASS == cls] + + if (nrow(choices) == 0) { + choices = global_probs[CLASS == cls] + } + + if (nrow(choices) == 0) { + choices = lookup_probs[CLASS == cls] + } + + if (!("prob" %in% names(choices))) { + choices[, prob := NA_real_] + } + + choices = choices[ + !is.na(SUBCLASS) & + !is.na(prob) & + is.finite(prob) & + prob > 0] - #if crop class transitions, draw subclass from new class distribution - preds_subclass$SUBCLASS[i] = draw_subclass(current_class) + 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 + ) + } } - #update previous class/subclass for next predicted year - prev_class = current_class - prev_subclass = preds_subclass$SUBCLASS[i] + 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_] + } + + dt = merge(dt, lookup_subclass, by = c("CLASS", "SUBCLASS"), all.x = TRUE) + + setorder(dt, orig_order) + dt[, orig_order := NULL] + + return(dt[])} +##--------till target functions--------- +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, + " with pattern: ", pattern)} + + 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[])} + +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 { + + temp[, rand_order := runif(.N)] + setorder(temp, rand_order) + + targets[, target_acres := till_share * group_total_acres] + targets[, upper_acres := cumsum(target_acres)] + + 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[])} -#adding class/subclass descriptions from lookup table -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 = merge( - preds_subclass, - lookup_subclass, - by = c("CLASS", "SUBCLASS"), - all.x = TRUE -) - -#adding parcel information from original design points -parcel_meta = design_points[ - order(year, season), - .SD[.N], - by = parcel_id -][ - , .(parcel_id, site_id, lon, lat) -] - -preds_landiq = merge( - preds_subclass, - parcel_meta, - by = "parcel_id", - all.x = TRUE -) - -#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" -) - -fwrite( - design_points_predicted, - sprintf("predicted_2024_%s.csv", end_year) -) +##--------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) +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))} +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, + 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] + } + + 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) + + ##--------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)} + + scenario_manifest = future_landiq[, .( + n_rows = .N, + n_parcels = uniqueN(parcel_id), + total_acres = sum(ACRES, na.rm = TRUE) + ), by = .(scenario, county_safe)] + + 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) +print(all_prediction_manifests) 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