Natural Resources Analytics¶

# Databricks notebook source
# MAGIC %md
# MAGIC # Silver → Gold: Seismic Cluster Detection
# MAGIC Identifies earthquake clusters using spatiotemporal density analysis
# MAGIC to detect swarm activity and induced seismicity patterns.

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

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

# Read Silver earthquake data (validated, deduplicated)
df_quakes = spark.read.format("delta").load("Tables/silver_doi_usgs_earthquakes")

# Filter to significant events (M2.5+) in CONUS
df_conus = (
    df_quakes
    .filter(
        (F.col("magnitude") >= 2.5) &
        (F.col("latitude").between(24.5, 49.5)) &
        (F.col("longitude").between(-125.0, -66.5))
    )
    .withColumn("event_date", F.to_date("event_time"))
)

# Spatial binning: 0.5-degree grid cells for cluster detection
df_gridded = (
    df_conus
    .withColumn("lat_bin", F.round(F.col("latitude") * 2) / 2)
    .withColumn("lon_bin", F.round(F.col("longitude") * 2) / 2)
    .withColumn("grid_id", F.concat_ws("_", F.col("lat_bin"), F.col("lon_bin")))
)

# Calculate cluster metrics per grid cell per month
df_clusters = (
    df_gridded
    .withColumn("year_month", F.date_format("event_date", "yyyy-MM"))
    .groupBy("grid_id", "lat_bin", "lon_bin", "year_month")
    .agg(
        F.count("*").alias("event_count"),
        F.max("magnitude").alias("max_magnitude"),
        F.avg("magnitude").alias("avg_magnitude"),
        F.avg("depth_km").alias("avg_depth_km"),
        F.min("event_time").alias("first_event"),
        F.max("event_time").alias("last_event"),
        F.stddev("magnitude").alias("magnitude_stddev"),
    )
)

# Classify cluster type based on depth and frequency patterns
# Shallow induced seismicity: depth < 10km, high frequency, near injection wells
df_classified = (
    df_clusters
    .withColumn(
        "cluster_type",
        F.when(
            (F.col("avg_depth_km") < 10) & (F.col("event_count") > 20),
            "POTENTIAL_INDUCED"
        )
        .when(F.col("event_count") > 50, "SWARM")
        .when(F.col("max_magnitude") >= 5.0, "SIGNIFICANT_SEQUENCE")
        .otherwise("BACKGROUND")
    )
    .withColumn(
        "b_value_proxy",
        F.log10(F.col("event_count")) / F.col("max_magnitude")
    )
)

df_classified.write.format("delta").mode("overwrite").option(
    "overwriteSchema", "true"
).saveAsTable("lh_gold.gold_doi_seismic_clusters")

induced = df_classified.filter(F.col("cluster_type") == "POTENTIAL_INDUCED").count()
print(f"Potential induced seismicity clusters: {induced}")

# Databricks notebook source
# MAGIC %md
# MAGIC # Silver → Gold: NWIS Streamflow Trend Analysis
# MAGIC Calculates long-term streamflow trends using Mann-Kendall
# MAGIC trend detection and Sen's slope estimation for drought/flood monitoring.

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

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

# Read Silver NWIS daily streamflow data
df_flow = spark.read.format("delta").load("Tables/silver_doi_nwis_streamflow")

# Calculate annual statistics per monitoring site
df_annual = (
    df_flow
    .withColumn("water_year", 
        F.when(F.month("measurement_date") >= 10, F.year("measurement_date") + 1)
        .otherwise(F.year("measurement_date"))
    )
    .groupBy("site_number", "site_name", "state", "huc8_code",
             "latitude", "longitude", "water_year")
    .agg(
        F.avg("discharge_cfs").alias("mean_annual_flow_cfs"),
        F.min("discharge_cfs").alias("min_flow_7day_cfs"),
        F.max("discharge_cfs").alias("peak_flow_cfs"),
        F.percentile_approx("discharge_cfs", 0.1).alias("flow_10th_pct"),
        F.percentile_approx("discharge_cfs", 0.9).alias("flow_90th_pct"),
        F.count("*").alias("observation_count"),
    )
)

# Calculate departure from long-term mean
site_window = Window.partitionBy("site_number")
annual_window = Window.partitionBy("site_number").orderBy("water_year")

df_trends = (
    df_annual
    .withColumn("long_term_mean", F.avg("mean_annual_flow_cfs").over(site_window))
    .withColumn("long_term_stddev", F.stddev("mean_annual_flow_cfs").over(site_window))
    .withColumn(
        "flow_anomaly_pct",
        (F.col("mean_annual_flow_cfs") - F.col("long_term_mean"))
        / F.col("long_term_mean") * 100
    )
    .withColumn(
        "standardized_anomaly",
        (F.col("mean_annual_flow_cfs") - F.col("long_term_mean"))
        / F.col("long_term_stddev")
    )
    .withColumn(
        "drought_flag",
        F.when(F.col("standardized_anomaly") < -1.5, "SEVERE_DROUGHT")
        .when(F.col("standardized_anomaly") < -1.0, "MODERATE_DROUGHT")
        .when(F.col("standardized_anomaly") > 1.5, "FLOOD_RISK")
        .otherwise("NORMAL")
    )
    .withColumn("flow_prev_year", F.lag("mean_annual_flow_cfs").over(annual_window))
    .withColumn(
        "yoy_change_pct",
        F.when(F.col("flow_prev_year") > 0,
            (F.col("mean_annual_flow_cfs") - F.col("flow_prev_year"))
            / F.col("flow_prev_year") * 100
        )
    )
)

df_trends.write.format("delta").mode("overwrite").option(
    "overwriteSchema", "true"
).saveAsTable("lh_gold.gold_doi_streamflow_trends")

drought = df_trends.filter(F.col("drought_flag") == "SEVERE_DROUGHT").count()
print(f"Site-years in severe drought: {drought}")

# Databricks notebook source
# MAGIC %md
# MAGIC # Gold: Endangered Species Habitat Risk Assessment
# MAGIC Correlates species occurrence with land use change, climate
# MAGIC indicators, and contamination data for habitat risk scoring.

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

from pyspark.sql import functions as F

# Read Silver species and habitat data
df_species = spark.read.format("delta").load("Tables/silver_doi_fws_species_listings")
df_habitat = spark.read.format("delta").load("Tables/silver_doi_fws_critical_habitat")
df_land_cover = spark.read.format("delta").load("Tables/silver_doi_usgs_land_cover")

# Species-level habitat risk assessment
df_risk = (
    df_species
    .join(df_habitat, on="species_code", how="left")
    .join(df_land_cover, on=["state_fips", "county_fips"], how="left")
    .withColumn(
        "habitat_loss_risk",
        F.when(F.col("developed_land_pct_change") > 5.0, "HIGH")
        .when(F.col("developed_land_pct_change") > 2.0, "MODERATE")
        .otherwise("LOW")
    )
    .withColumn(
        "climate_vulnerability",
        F.when(F.col("temp_anomaly_c") > 2.0, "HIGH")
        .when(F.col("temp_anomaly_c") > 1.0, "MODERATE")
        .otherwise("LOW")
    )
    .withColumn(
        "contamination_exposure",
        F.when(F.col("tri_facilities_in_habitat") > 5, "HIGH")
        .when(F.col("tri_facilities_in_habitat") > 0, "MODERATE")
        .otherwise("LOW")
    )
)

# Calculate composite habitat risk score (0-100)
risk_map = {"HIGH": 3, "MODERATE": 2, "LOW": 1}

df_scored = (
    df_risk
    .withColumn("habitat_score",
        F.when(F.col("habitat_loss_risk") == "HIGH", 3)
        .when(F.col("habitat_loss_risk") == "MODERATE", 2).otherwise(1))
    .withColumn("climate_score",
        F.when(F.col("climate_vulnerability") == "HIGH", 3)
        .when(F.col("climate_vulnerability") == "MODERATE", 2).otherwise(1))
    .withColumn("contam_score",
        F.when(F.col("contamination_exposure") == "HIGH", 3)
        .when(F.col("contamination_exposure") == "MODERATE", 2).otherwise(1))
    .withColumn(
        "composite_risk_score",
        F.round(
            (F.col("habitat_score") * 0.4 +
             F.col("climate_score") * 0.35 +
             F.col("contam_score") * 0.25) / 3 * 100
        ).cast("int")
    )
    .withColumn(
        "recovery_priority",
        F.when(F.col("composite_risk_score") >= 80, "CRITICAL")
        .when(F.col("composite_risk_score") >= 60, "HIGH")
        .when(F.col("composite_risk_score") >= 40, "MODERATE")
        .otherwise("STABLE")
    )
)

df_scored.write.format("delta").mode("overwrite").option(
    "overwriteSchema", "true"
).saveAsTable("lh_gold.gold_doi_species_habitat_risk")

critical = df_scored.filter(F.col("recovery_priority") == "CRITICAL").count()
print(f"Species with CRITICAL habitat risk: {critical}")

# Cross-domain: NOAA climate indicators vs. USGS streamflow
df_cross_climate = (
    df_noaa_climate_normals
    .join(df_doi_streamflow_trends, on=["state", "huc8_code", "water_year"])
    .select(
        "state", "huc8_code", "water_year",
        "precip_anomaly_pct", "temp_anomaly_c", "pdsi_value",
        "flow_anomaly_pct", "drought_flag",
    )
    .withColumn(
        "climate_flow_correlation",
        F.corr("precip_anomaly_pct", "flow_anomaly_pct")
            .over(Window.partitionBy("huc8_code"))
    )
)

# Cross-domain: EPA contamination within DOI critical habitat
df_cross_contam = (
    df_epa_tri_facilities
    .join(
        df_doi_critical_habitat,
        (F.col("tri.latitude").between(
            F.col("habitat.bbox_south"), F.col("habitat.bbox_north"))) &
        (F.col("tri.longitude").between(
            F.col("habitat.bbox_west"), F.col("habitat.bbox_east")))
    )
    .groupBy("species_code", "common_name", "listing_status")
    .agg(
        F.count("trifid").alias("tri_facilities_in_habitat"),
        F.sum("total_releases_lbs").alias("total_releases_in_habitat"),
        F.countDistinct("chemical_name").alias("distinct_contaminants"),
    )
    .orderBy(F.col("total_releases_in_habitat").desc())
)