Agricultural Analytics¶

# Databricks notebook source
# MAGIC %md
# MAGIC # USDA Crop Production Trend Analysis
# MAGIC Analyze crop yield, acreage, and production trends from NASS QuickStats data.

# COMMAND ----------

from pyspark.sql import functions as F
from pyspark.sql.window import Window

# Load NASS crop production data from bronze layer
df_crops = spark.read.format("delta").load(
    "Tables/bronze_usda_crop_production"
)

# Focus on major field crops at state level
major_crops = ["CORN", "SOYBEANS", "WHEAT", "COTTON", "RICE"]

df_crop_trends = (
    df_crops
    .filter(F.col("commodity_desc").isin(major_crops))
    .filter(F.col("agg_level_desc") == "STATE")
    .filter(F.col("statisticcat_desc") == "YIELD")
    .filter(F.col("unit_desc").contains("BU / ACRE"))
    .filter(F.col("year").between(2000, 2025))
    .withColumn("value_numeric", F.col("value").cast("double"))
    .filter(F.col("value_numeric").isNotNull())
)

# Calculate state-level yield trends with moving averages
window_3yr = Window.partitionBy("state_name", "commodity_desc").orderBy("year").rowsBetween(-2, 0)
window_trend = Window.partitionBy("state_name", "commodity_desc").orderBy("year")

df_yield_analysis = (
    df_crop_trends
    .groupBy("year", "state_name", "commodity_desc")
    .agg(F.avg("value_numeric").alias("avg_yield_bu_acre"))
    .withColumn("moving_avg_3yr", F.round(F.avg("avg_yield_bu_acre").over(window_3yr), 2))
    .withColumn("yield_yoy_change", F.round(
        F.col("avg_yield_bu_acre") - F.lag("avg_yield_bu_acre").over(window_trend), 2
    ))
    .withColumn("yield_yoy_pct", F.round(
        F.col("yield_yoy_change") / F.lag("avg_yield_bu_acre").over(window_trend) * 100, 2
    ))
)

# Identify states with declining yield trends (potential climate stress)
df_declining_yields = (
    df_yield_analysis
    .filter(F.col("year") >= 2020)
    .groupBy("state_name", "commodity_desc")
    .agg(
        F.avg("yield_yoy_pct").alias("avg_annual_yield_change_pct"),
        F.min("avg_yield_bu_acre").alias("min_yield"),
        F.max("avg_yield_bu_acre").alias("max_yield")
    )
    .filter(F.col("avg_annual_yield_change_pct") < 0)
    .orderBy("avg_annual_yield_change_pct")
)

display(df_declining_yields)

# Write to silver layer
df_yield_analysis.write.format("delta").mode("overwrite").saveAsTable(
    "lh_silver.silver_usda_crop_yield_trends"
)

# Databricks notebook source
# MAGIC %md
# MAGIC # FSIS Food Safety Recall Pattern Detection
# MAGIC Analyze recall frequency, severity, reason, and establishment patterns.

# COMMAND ----------

from pyspark.sql import functions as F

# Load FSIS recall data
df_recalls = spark.read.format("delta").load(
    "Tables/bronze_usda_fsis_recalls"
)

# Parse and enrich recall records
df_recalls_enriched = (
    df_recalls
    .withColumn("recall_date", F.to_date("recall_date_str", "yyyy-MM-dd"))
    .withColumn("recall_year", F.year("recall_date"))
    .withColumn("recall_month", F.month("recall_date"))
    .withColumn("pounds_recalled_num", F.col("pounds_recalled").cast("double"))
    .filter(F.col("recall_year") >= 2015)
)

# Recall trends by reason and classification
df_recall_trends = (
    df_recalls_enriched
    .groupBy("recall_year", "recall_classification", "reason_for_recall")
    .agg(
        F.count("*").alias("recall_count"),
        F.sum("pounds_recalled_num").alias("total_pounds_recalled"),
        F.countDistinct("establishment_name").alias("unique_establishments"),
        F.avg("pounds_recalled_num").alias("avg_pounds_per_recall")
    )
    .orderBy("recall_year", F.desc("recall_count"))
)

# Identify repeat offender establishments
df_repeat_establishments = (
    df_recalls_enriched
    .groupBy("establishment_name", "establishment_number", "state")
    .agg(
        F.count("*").alias("total_recalls"),
        F.sum(F.when(F.col("recall_classification") == "Class I", 1).otherwise(0)).alias("class_i_recalls"),
        F.sum("pounds_recalled_num").alias("total_pounds_recalled"),
        F.min("recall_date").alias("first_recall"),
        F.max("recall_date").alias("latest_recall"),
        F.collect_set("reason_for_recall").alias("recall_reasons")
    )
    .filter(F.col("total_recalls") >= 3)
    .orderBy(F.desc("total_recalls"))
)

# Allergen vs. pathogen recall ratio over time
df_allergen_vs_pathogen = (
    df_recalls_enriched
    .withColumn("recall_category", F.when(
        F.col("reason_for_recall").contains("allergen") | F.col("reason_for_recall").contains("undeclared"),
        "Undeclared Allergen"
    ).when(
        F.col("reason_for_recall").rlike("(?i)salmonella|listeria|e\\.\\s*coli|clostridium"),
        "Pathogen Contamination"
    ).otherwise("Other"))
    .groupBy("recall_year", "recall_category")
    .agg(F.count("*").alias("recall_count"))
)

display(df_allergen_vs_pathogen)

# Write to silver layer
df_recall_trends.write.format("delta").mode("overwrite").saveAsTable(
    "lh_silver.silver_usda_fsis_recall_trends"
)
df_repeat_establishments.write.format("delta").mode("overwrite").saveAsTable(
    "lh_silver.silver_usda_fsis_repeat_establishments"
)

# Databricks notebook source
# MAGIC %md
# MAGIC # SNAP Benefit Distribution Geographic Analysis
# MAGIC Analyze SNAP participation and benefits by state with food insecurity indicators.

# COMMAND ----------

from pyspark.sql import functions as F
from pyspark.sql.window import Window

# Load SNAP participation data
df_snap = spark.read.format("delta").load(
    "Tables/bronze_usda_snap_participation"
)

# Load state population estimates for per-capita analysis
df_population = spark.read.format("delta").load(
    "Tables/bronze_census_state_population"
)

# Monthly SNAP analysis by state
df_snap_analysis = (
    df_snap
    .filter(F.col("fiscal_year") >= 2018)
    .withColumn("avg_persons_num", F.col("avg_participation_persons").cast("double"))
    .withColumn("total_benefits_num", F.col("total_benefits_issued").cast("double"))
    .withColumn("avg_benefit_per_person",
        F.round(F.col("total_benefits_num") / F.col("avg_persons_num"), 2)
    )
)

# State-level annual aggregation with per-capita participation rate
df_state_snap = (
    df_snap_analysis
    .groupBy("fiscal_year", "state_name")
    .agg(
        F.avg("avg_persons_num").alias("avg_monthly_participants"),
        F.sum("total_benefits_num").alias("annual_benefits_issued"),
        F.avg("avg_benefit_per_person").alias("avg_monthly_benefit_per_person")
    )
    .join(
        df_population.select("year", "state_name", "population"),
        (F.col("fiscal_year") == F.col("year")) & (df_snap_analysis["state_name"] == df_population["state_name"]),
        "left"
    )
    .withColumn("participation_rate_pct",
        F.round(F.col("avg_monthly_participants") / F.col("population") * 100, 2)
    )
)

# Identify food insecurity hotspots (high participation rate states)
df_food_insecurity = (
    df_state_snap
    .filter(F.col("fiscal_year") == 2023)
    .orderBy(F.desc("participation_rate_pct"))
)

# Trend analysis — participation change since pre-pandemic baseline
window_trend = Window.partitionBy("state_name").orderBy("fiscal_year")
df_snap_trends = (
    df_state_snap
    .withColumn("participation_change_pct", F.round(
        (F.col("avg_monthly_participants") - F.first("avg_monthly_participants").over(
            Window.partitionBy("state_name").orderBy("fiscal_year").rowsBetween(Window.unboundedPreceding, Window.unboundedPreceding)
        )) / F.first("avg_monthly_participants").over(
            Window.partitionBy("state_name").orderBy("fiscal_year").rowsBetween(Window.unboundedPreceding, Window.unboundedPreceding)
        ) * 100, 2
    ))
)

display(df_food_insecurity)

# Write to gold layer
df_state_snap.write.format("delta").mode("overwrite").saveAsTable(
    "lh_gold.gold_usda_snap_state_analysis"
)

# Cross-domain: USDA × EPA agricultural runoff analysis
df_crop_acreage = spark.read.format("delta").load("Tables/bronze_usda_crop_production")
df_water_quality = spark.read.format("delta").load("Tables/bronze_epa_water_quality")

# Calculate agricultural intensity per county (acres planted / total land area)
df_ag_intensity = (
    df_crop_acreage
    .filter(F.col("statisticcat_desc") == "AREA PLANTED")
    .filter(F.col("agg_level_desc") == "COUNTY")
    .groupBy("year", "state_name", "county_name", "state_fips_code", "county_code")
    .agg(F.sum("value_numeric").alias("total_acres_planted"))
)

# Join with EPA water quality for downstream nitrate measurements
df_runoff_correlation = (
    df_ag_intensity
    .join(
        df_water_quality
        .filter(F.col("parameter_name") == "Nitrate")
        .groupBy("year", "county_name", "state_name")
        .agg(F.avg("result_value").alias("avg_nitrate_mg_l")),
        on=["year", "state_name", "county_name"],
        how="inner"
    )
)

correlation = df_runoff_correlation.stat.corr("total_acres_planted", "avg_nitrate_mg_l")
print(f"Agricultural Intensity-Nitrate Correlation: {correlation:.4f}")

# Cross-domain: USDA × NOAA weather-crop yield correlation
df_yields = spark.read.format("delta").load("Tables/silver_usda_crop_yield_trends")
df_drought = spark.read.format("delta").load("Tables/bronze_noaa_drought_index")

# Join crop yields with growing season (Apr-Sep) drought severity
df_weather_yield = (
    df_yields
    .filter(F.col("commodity_desc") == "CORN")
    .join(
        df_drought
        .filter(F.col("month").between(4, 9))
        .groupBy("year", "state_name")
        .agg(F.avg("pdsi").alias("growing_season_avg_pdsi")),
        on=["year", "state_name"],
        how="inner"
    )
)

# Quantify yield impact per unit PDSI change
correlation = df_weather_yield.stat.corr("growing_season_avg_pdsi", "avg_yield_bu_acre")
print(f"PDSI-Corn Yield Correlation: {correlation:.4f}")