jeffrey1963/Cattle-Elk-R5

1def initialize_parcels(parcel_map, cluster_labels):2    parcel_dict = {}3    for i in range(parcel_map.shape[0]):4        for j in range(parcel_map.shape[1]):5            cluster_id = parcel_map[i, j]6            land_type = cluster_labels.get(cluster_id, "Unknown")7 8            parcel_dict[(i, j)] = {9                "land_type": land_type,10                "forage": None,11                "health": 1.0,12                "degraded": False,13                "cattle_grazing": {14                    "conservative": land_type == "Productive Grass",15                    "moderate": land_type in ["Productive Grass", "Pasture/Desert"],16                    "aggressive": land_type not in ["Water"]17                },18                "elk_grazing": {19                    "default": land_type in ["Riparian Sensitive Zone", "Productive Grass"]20                }21            }22    return parcel_dict23 24def get_land_forage_rates():25    return {26        'Productive Grass': 1.0,27        'Pasture/Desert': 0.4,28        'Riparian Sensitive Zone': 1.2,29        'Rocky Area': 0.2,30        'Water': 0.031    }32 33 34def assign_initial_forage(parcel_dict, land_forage_rates):35    for parcel in parcel_dict.values():36        rate = land_forage_rates.get(parcel["land_type"], 0.0)37        parcel["forage"] = rate * 100  # initial AUMs38 39import numpy as np40import numpy as np41import matplotlib.pyplot as plt42from collections import Counter43from matplotlib.colors import ListedColormap44import matplotlib.patches as mpatches45 46def simulate_period(parcel_dict, grazing_strategy="moderate", cattle_stocking_rate=100000, elk_pressure=3000):47    print(f"\n🟢 Running LP-based simulation using grazing strategy: **{grazing_strategy.upper()}**")48 49    import numpy as np50    from scipy.optimize import linprog51 52    keys = list(parcel_dict.keys())53    n_rows = max(i for i, _ in keys) + 154    n_cols = max(j for _, j in keys) + 155    num_cells = n_rows * n_cols56 57    # Step 1: Regrowth58    land_growth_rates = {59        'Productive Grass': 1.0,60        'Pasture/Desert': 0.4,61        'Riparian Sensitive Zone': 1.2,62        'Rocky Area': 0.2,63        'Water': 0.064    }65    for parcel in parcel_dict.values():66        base_growth = land_growth_rates.get(parcel["land_type"], 0.0)67        weather = np.random.normal(1.0, 0.15)68        regrowth = base_growth * weather * 169        regrowth *= parcel["health"]70        parcel["forage"] = min(parcel["forage"] + regrowth, 100)71 72    # Step 2: Subtract uniform elk grazing from all parcels73    elk_grazing_per_parcel = elk_pressure / num_cells74    for parcel in parcel_dict.values():75        parcel["forage"] -= elk_grazing_per_parcel76        parcel["forage"] = max(parcel["forage"], 0.0)77 78    # Step 3: LP for cattle79    cost = []80    bounds = []81    eligible_keys = []82    for i in range(n_rows):83        for j in range(n_cols):84            p = parcel_dict[(i, j)]85            if not p["cattle_grazing"].get(grazing_strategy, False):86                cost.append(0)87                bounds.append((0, 0))88                continue89            if grazing_strategy in {"conservative", "moderate"} and p["land_type"] == "Riparian Sensitive Zone":90                cost.append(0)91                bounds.append((0, 0))92                continue93            cost.append((i + j) * 0.02)94            bounds.append((0, p["forage"]))95            eligible_keys.append((i, j))96 97    A_eq = [1.0 if b[1] > 0 else 0.0 for b in bounds]98    b_eq = [cattle_stocking_rate]99 100    result = linprog(c=cost, A_eq=[A_eq], b_eq=b_eq, bounds=bounds, method="highs")101    if not result.success:102        raise RuntimeError("Grazing LP failed: " + result.message)103 104    grazing_values = result.x105 106    # Step 4: Apply grazing and your specified health rule107    for idx, ((i, j), x) in enumerate(zip(parcel_dict.keys(), grazing_values)):108        parcel = parcel_dict[(i, j)]109        parcel["forage"] -= x110 111        # ✅ Your health rule (fully time-dynamic)112        if parcel["forage"] <= 0:113            parcel["health"] = max(parcel["health"] - 0.25, 0.0)114        elif parcel["forage"] < 20:115            parcel["health"] = max(parcel["health"] - 0.1, 0.0)116        else:117            parcel["health"] = min(parcel["health"] + 0.02, 1.0)118 119def simulate_periodold(parcel_dict, grazing_strategy="moderate", cattle_stocking_rate=5000, elk_pressure=3000):120    print(f"\n🟢 Running simulation using cattle grazing strategy: **{grazing_strategy.upper()}**")121 122    land_growth_rates = {123        'Productive Grass': 1.0,124        'Pasture/Desert': 0.4,125        'Riparian Sensitive Zone': 1.2,126        'Rocky Area': 0.2,127        'Water': 0.0128    }129 130    # 1. Simulate forage regrowth for 8 months131    for parcel in parcel_dict.values():132        base_growth = land_growth_rates.get(parcel["land_type"], 0.0)133        weather = np.random.normal(1.0, 0.15)134        regrowth = base_growth * weather * 8  # ← 8 months, as you said135        regrowth *= parcel["health"]  # degrade means slower regrowth136        parcel["forage"] = min(parcel["forage"] + regrowth, 100)137 138    # 2. Count eligible parcels139    total_grazed_parcels = sum(140        1 for parcel in parcel_dict.values() if parcel["cattle_grazing"].get(grazing_strategy, False)141    )142    if total_grazed_parcels == 0:143        print("⚠️ No parcels match the selected grazing strategy.")144        return145 146    cattle_grazing_per_parcel = cattle_stocking_rate / total_grazed_parcels147    elk_grazing_per_parcel = elk_pressure / len(parcel_dict)148 149    # 3. Simulate grazing and degradation150    for parcel in parcel_dict.values():151        if not parcel["cattle_grazing"].get(grazing_strategy, False):152            continue153 154        total_grazing = cattle_grazing_per_parcel + elk_grazing_per_parcel155 156        if total_grazing > parcel["forage"]:157            parcel["degraded"] = True158            parcel["health"] = max(parcel["health"] - 0.1, 0.0)159        else:160            parcel["health"] = min(parcel["health"] + 0.02, 1.0)161 162        parcel["forage"] = max(parcel["forage"] - total_grazing, 0)163 164def get_forage_map(parcel_dict, n_rows, n_cols):165    return np.array([[parcel_dict[(i, j)]["forage"] for j in range(n_cols)] for i in range(n_rows)])166 167def get_health_map(parcel_dict, n_rows, n_cols):168    """169    Returns a 2D numpy array representing the health of each parcel.170    """171    return np.array([[parcel_dict[(i, j)]["health"] for j in range(n_cols)] for i in range(n_rows)])172 173 174def plot_health_map(health_map, title="Parcel Health Levels", save_path=None):175    """176    Plots a heatmap of the parcel health values.177    """178    import matplotlib.pyplot as plt179 180    plt.figure(figsize=(8, 6))181    plt.imshow(health_map, cmap='RdYlGn', origin='upper', vmin=0, vmax=1)182    plt.colorbar(label="Health Index (0–1)")183    plt.title(title)184    plt.axis('off')185    plt.tight_layout()186    if save_path:187        plt.savefig(save_path)188    plt.show()189 190 191 192def plot_forage_map(forage_map, title="Parcel Forage Levels"):193    plt.figure(figsize=(8, 6))194    plt.imshow(forage_map, cmap='YlGn', origin='upper')195    plt.colorbar(label="Forage AUMs")196    plt.title(title)197    plt.axis('off')198    plt.tight_layout()199    plt.show()200 201def run_full_simulation(parcel_map, cluster_labels, n_rows, n_cols, strategy="moderate"):202    parcel_dict = initialize_parcels(parcel_map, cluster_labels)203    land_forage_rates = get_land_forage_rates()204    assign_initial_forage(parcel_dict, land_forage_rates)205    simulate_period(parcel_dict, grazing_strategy=strategy)206    return parcel_dict207 208