Yield map from KBS LTER harvest data

Python Code

"""Precision Agriculture — Yield map from KBS LTER harvest data."""


import numpy as np
import pandas as pd
import plotly.graph_objects as go
import io, requests, warnings


URL = "https://lter.kbs.msu.edu/datatables/80.csv"

LON_C, LAT_C = -85.37, 42.39


def fetch_real():
    warnings.filterwarnings("ignore")
    requests.packages.urllib3.disable_warnings()
    r = requests.get(URL, verify=False, timeout=60)
    r.raise_for_status()
    df = pd.read_csv(io.StringIO(r.text), comment="#", low_memory=False)
    df.columns = df.columns.str.strip().str.lower().str.replace(" ", "_")
    return df


def make_synthetic():
    RNG = np.random.default_rng(3)
    n = 5000
    lons = LON_C + RNG.uniform(-0.003, 0.003, n)
    lats = LAT_C + RNG.uniform(-0.002, 0.002, n)
    base = 45 + 8 * (lats - LAT_C) / 0.002 + 6 * (lons - LON_C) / 0.003
    yields = np.clip(base + RNG.normal(0, 4, n), 20, 75)
    return pd.DataFrame({"longitude": lons, "latitude": lats, "yield_bu_ac": yields})


def generate():
    print("building yield map …")
    df = None
    try:
        raw = fetch_real()
        col_map = {}
        for c in raw.columns:
            if "lat" in c: col_map["latitude"] = c
            if "lon" in c or "lng" in c: col_map["longitude"] = c
        if len(col_map) == 2:
            raw = raw.rename(columns={v: k for k, v in col_map.items()})
        if "latitude" in raw.columns and "longitude" in raw.columns:
            flow_cols = [c for c in raw.columns
                         if "flow" in c or "yield" in c or "crop" in c]
            if flow_cols:
                raw["yield_bu_ac"] = pd.to_numeric(raw[flow_cols[0]], errors="coerce")
                raw = raw.dropna(subset=["latitude", "longitude", "yield_bu_ac"])
                raw = raw[raw["yield_bu_ac"] > 0]
                if len(raw) > 100:
                    if len(raw) > 5000:
                        raw = raw.sample(5000, random_state=42)
                    df = raw[["latitude", "longitude", "yield_bu_ac"]]
                    print(f"  real: {len(df)} points")
    except Exception as e:
        print(f"  fetch failed ({e})")

    if df is None:
        df = make_synthetic()
        print(f"  synthetic: {len(df)} points")

    center_lat = float(df["latitude"].mean())
    center_lon = float(df["longitude"].mean())
    q02 = float(df["yield_bu_ac"].quantile(0.02))
    q98 = float(df["yield_bu_ac"].quantile(0.98))

    fig = go.Figure(go.Scattermapbox(
        lat=df["latitude"],
        lon=df["longitude"],
        mode="markers",
        marker=dict(
            size=3,
            color=df["yield_bu_ac"],
            colorscale="YlOrRd",
            cmin=q02,
            cmax=q98,
            colorbar=dict(title="Yield (bu/ac)", thickness=14),
            opacity=0.9,
        ),
        hovertemplate="Lon: %{lon:.4f}<br>Lat: %{lat:.4f}<br>Yield: %{marker.color:.1f} bu/ac<extra></extra>",
    ))

    fig.update_layout(
        mapbox=dict(
            style="white-bg",
            layers=[{
                "below": "traces",
                "sourcetype": "raster",
                "sourceattribution": "Tiles © Esri",
                "source": ["https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}"],
            }],
            center=dict(lat=center_lat, lon=center_lon),
            zoom=16,
        ),
        margin=dict(t=0, b=0, l=0, r=0),
    )
    return fig


if __name__ == "__main__":
    generate()

Made with Plotly