Script: Geological Data and Mineral Map Generator

#!/usr/bin/env python3
"""
Generate geological reference data for the Recovery Library.

Produces:
  - tables-geological.md  (NZ mineral deposits, critical gaps, aggregates,
                           geothermal resources, volcanic hazard zones)
  - site/images/geological-minerals.png  (mineral deposit map)

Uses: matplotlib, json, standard library.

Usage:
    python generate_geological_data.py [output_md] [output_png]
    Defaults:
        ../tables-geological.md
        ../site/images/geological-minerals.png

Data sources: Crown Minerals / MBIE, GNS Science.

Run with:
    scripts/.venv/bin/python generate_geological_data.py
"""

import json
import os
import sys
from datetime import datetime

# ---------------------------------------------------------------------------
# Data definitions
# ---------------------------------------------------------------------------

# NZ Mineral Deposits
# Fields: mineral, location, lon, lat, resource_estimate, extraction_method,
#         current_status, feasibility
MINERAL_DEPOSITS = [
    # [data omitted for brevity — 23 entries covering ironsand, coal, gold,
    #  limestone, silica sand, clay/kaolin, aggregate, pounamu/nephrite,
    #  bentonite, zeolite, phosphate, and serpentinite deposits]
]

# Critical Mineral Gap Analysis
# Fields: mineral, nz_source, annual_consumption_est, strategic_importance,
#         substitution_options, recommended_stockpile
CRITICAL_GAPS = [
    # [data omitted for brevity — 12 entries covering copper, tin, zinc, lead,
    #  aluminium, chromium, manganese, nickel, lithium, rare earths, sulfur,
    #  and potash]
]

# Aggregate and Building Stone
# Fields: region, location, rock_type, suitability
AGGREGATES = [
    # [data omitted for brevity — 17 entries covering quarry locations from
    #  Northland to Southland]
]

# Geothermal Resources
# Fields: name, lon, lat, temperature_C, current_use, expansion_potential
GEOTHERMAL = [
    # [data omitted for brevity — 12 entries covering TVZ and Northland
    #  geothermal fields]
]

# Volcanic Hazard Zones
# Fields: volcano, lon, lat, type, hazard_radius_km, last_eruption, notes
VOLCANIC_HAZARDS = [
    # [data omitted for brevity — 11 entries covering active and dormant
    #  NZ volcanoes]
]

# Deposit types for the map (lon, lat, marker_color, marker_style, label)
MAP_POINTS = [
    # [data omitted for brevity — 21 entries mapping mineral deposit locations
    #  with marker styles by type: ironsand, coal, gold, limestone, silica,
    #  clay/kaolin, pounamu, zeolite, bentonite, phosphate, serpentinite]
]


# ---------------------------------------------------------------------------
# Markdown generation
# ---------------------------------------------------------------------------

def generate_markdown():
    now = datetime.now().strftime("%Y-%m-%d")
    lines = []

    lines.append("---")
    lines.append('title: "Geological Atlas — NZ Mineral and Geothermal Data"')
    lines.append('subtitle: "Recovery Library — Doc \\#022"')
    lines.append(f'date: "{now}"')
    lines.append("---")
    lines.append("")
    lines.append("# Geological Atlas — NZ Mineral and Geothermal Data")
    lines.append("")
    lines.append("*Recovery Library — Doc #022*")
    lines.append("")
    lines.append(
        "**Sources:** Crown Minerals / MBIE mineral permit database; GNS Science "
        "geological maps and geothermal resource assessments; NZ Petroleum & Minerals "
        "annual reviews. Resource estimates are order-of-magnitude figures; consult "
        "original survey reports before operational decisions."
    )
    lines.append("")
    lines.append(
        "**Feasibility ratings:** [A] = recoverable with pre-collapse infrastructure; "
        "[B] = recoverable with moderate effort/investment; "
        "[C] = technically possible but significant barriers; "
        "[D] = not realistically recoverable under near-term recovery conditions."
    )
    lines.append("")
    lines.append("![NZ Mineral Deposits](images/geological-minerals.png)")
    lines.append("")
    lines.append("---")
    lines.append("")

    # -----------------------------------------------------------------------
    # Table 1: NZ Mineral Deposits
    # -----------------------------------------------------------------------
    lines.append("## 1. NZ Mineral Deposits")
    lines.append("")
    lines.append(
        "| Mineral | Location | Coordinates | Estimated Resource | "
        "Extraction Method | Current Status | Feasibility |"
    )
    lines.append(
        "|---------|----------|-------------|-------------------|"
        "-------------------|----------------|-------------|"
    )
    for row in MINERAL_DEPOSITS:
        mineral, location, lon, lat, resource, method, status, feas = row
        coord = f"{abs(lat):.2f}°S {lon:.2f}°E"
        lines.append(
            f"| {mineral} | {location} | {coord} | {resource} | "
            f"{method} | {status} | [{feas}] |"
        )
    lines.append("")
    lines.append(
        "*Footnote 1:* Ironsand resource figures from Crown Minerals NZ Minerals "
        "Programme 2023. Gold resource figures from NZX disclosures (OceanaGold) "
        "and GNS Science ore deposit database. Phosphate from Chatham Rock Phosphate "
        "prospectus and EPA decision (2015). Pounamu ownership vested in Ngāi Tahu "
        "under the Ngāi Tahu Claims Settlement Act 1998 — any recovery-era use "
        "requires Ngāi Tahu governance.[^1]"
    )
    lines.append("")
    lines.append("---")
    lines.append("")

    # -----------------------------------------------------------------------
    # Table 2: Critical Mineral Gap Analysis
    # -----------------------------------------------------------------------
    lines.append("## 2. Critical Mineral Gap Analysis")
    lines.append("")
    lines.append(
        "Minerals NZ does not produce in economically significant quantities and must "
        "stockpile or substitute. Annual consumption figures are approximate national "
        "estimates based on Statistics NZ trade data and industry reports."
    )
    lines.append("")
    lines.append(
        "| Mineral | Current NZ Source | Annual Consumption Est. | "
        "Strategic Importance | Substitution Options | Recommended Stockpile |"
    )
    lines.append(
        "|---------|-------------------|-------------------------|"
        "---------------------|----------------------|----------------------|"
    )
    for row in CRITICAL_GAPS:
        mineral, source, consumption, importance, subs, stockpile = row
        lines.append(
            f"| {mineral} | {source} | {consumption} | "
            f"{importance} | {subs} | {stockpile} |"
        )
    lines.append("")
    lines.append(
        "*Footnote 2:* Consumption estimates derived from Statistics NZ overseas "
        "trade data (HS commodity codes) and MBIE energy and minerals statistics. "
        "Potash and sulfur figures are particularly critical — NZ superphosphate "
        "industry depends entirely on imported sulfur and rock phosphate. "
        "Tiwai Point aluminium smelter closure (scheduled 2024) ends domestic "
        "primary aluminium production.[^2]"
    )
    lines.append("")
    lines.append("---")
    lines.append("")

    # -----------------------------------------------------------------------
    # Table 3: Aggregate and Building Stone
    # -----------------------------------------------------------------------
    lines.append("## 3. Aggregate and Building Stone — Major Quarry Locations")
    lines.append("")
    lines.append(
        "| Region | Location | Rock Type | Suitability |"
    )
    lines.append(
        "|--------|----------|-----------|-------------|"
    )
    for region, location, rock, suitability in AGGREGATES:
        lines.append(f"| {region} | {location} | {rock} | {suitability} |")
    lines.append("")
    lines.append(
        "*Footnote 3:* Quarry locations and rock types from GNS Science QMAP "
        "geological map series and regional council aggregate surveys. "
        "River gravel extraction from Waikato and Canterbury rivers is "
        "increasingly consent-constrained; alternative hard-rock quarry capacity "
        "should be maintained.[^3]"
    )
    lines.append("")
    lines.append("---")
    lines.append("")

    # -----------------------------------------------------------------------
    # Table 4: Geothermal Resources
    # -----------------------------------------------------------------------
    lines.append("## 4. Geothermal Resources")
    lines.append("")
    lines.append(
        "All major fields are within the Taupō Volcanic Zone (TVZ) except Ngāwha "
        "(Northland). Temperatures are approximate reservoir values; wellhead "
        "temperatures vary. Current installed capacity ~1,000 MWe nationally."
    )
    lines.append("")
    lines.append(
        "| Field | Coordinates | Reservoir Temp. (°C) | Current Use | Expansion Potential |"
    )
    lines.append(
        "|-------|-------------|----------------------|-------------|---------------------|"
    )
    for row in GEOTHERMAL:
        name, lon, lat, temp, use, potential = row
        coord = f"{abs(lat):.2f}°S {lon:.2f}°E"
        lines.append(f"| {name} | {coord} | {temp} | {use} | {potential} |")
    lines.append("")
    lines.append(
        "*Footnote 4:* Reservoir temperatures and installed capacity from GNS Science "
        "geothermal programme reports and generator company annual reports (Contact "
        "Energy, Mercury NZ, Top Energy). Geothermal generation contributes ~17% of "
        "NZ electricity; resilient to climate variability — a key post-disruption "
        "energy advantage.[^4]"
    )
    lines.append("")
    lines.append("---")
    lines.append("")

    # -----------------------------------------------------------------------
    # Table 5: Volcanic Hazard Zones
    # -----------------------------------------------------------------------
    lines.append("## 5. Volcanic Hazard Zones")
    lines.append("")
    lines.append(
        "Hazard radii indicate the primary risk zone (pyroclastic flows, lahars, "
        "ballistics, or surge depending on volcano type). Ashfall can extend "
        "hundreds of kilometres beyond these radii. The Taupō supervolcano "
        "hazard radius reflects potential national impact from a large (VEI 6+) "
        "eruption."
    )
    lines.append("")
    lines.append(
        "| Volcano | Coordinates | Type | Primary Hazard Radius (km) | "
        "Last Significant Eruption | Notes |"
    )
    lines.append(
        "|---------|-------------|------|---------------------------|"
        "--------------------------|-------|"
    )
    for row in VOLCANIC_HAZARDS:
        name, lon, lat, vtype, radius, last, notes = row
        coord = f"{abs(lat):.2f}°S {lon:.2f}°E"
        lines.append(
            f"| {name} | {coord} | {vtype} | {radius} | {last} | {notes} |"
        )
    lines.append("")
    lines.append(
        "*Footnote 5:* Hazard radii from GNS Science volcanic hazard models and "
        "National Hazards Centre reports. Auckland Volcanic Field represents a "
        "unique planning challenge — the next eruption could occur at any location "
        "within the ~360 km² field. Taupō has erupted at VEI 8 (Oruanui, ~26,500 BP) "
        "and VEI 7 (Hatepe, ~232 CE); a repeat would be a civilisation-scale event "
        "for the Southern Hemisphere.[^5]"
    )
    lines.append("")
    lines.append("---")
    lines.append("")

    # Footnotes block
    lines.append("[^1]: Crown Minerals Act 1991; GNS Science Mineral Deposit Database; "
                 "OceanaGold NZX disclosures 2023; Ngāi Tahu Claims Settlement Act 1998; "
                 "Chatham Rock Phosphate EPA Decision EPBC 2012/0635.")
    lines.append("")
    lines.append("[^2]: Statistics NZ Overseas Merchandise Trade, HS codes 26–28, 31, 75–81 "
                 "(2022 calendar year); MBIE Energy in NZ 2023; Rio Tinto Tiwai closure announcement.")
    lines.append("")
    lines.append("[^3]: GNS Science QMAP 1:250,000 geological map series; "
                 "Auckland Council aggregate demand study 2021; "
                 "Environment Canterbury aggregate resources review 2019.")
    lines.append("")
    lines.append("[^4]: GNS Science geothermal programme annual report 2023; "
                 "Contact Energy generation portfolio report 2023; "
                 "Mercury NZ annual report 2023; Top Energy Ngāwha geothermal project.")
    lines.append("")
    lines.append("[^5]: GNS Science volcanic hazard models (Volcanic Alert Levels); "
                 "National Emergency Management Agency volcanic hazard guides; "
                 "Wilson et al. (2009) Taupō eruptive history, *J. Volcanology*; "
                 "Auckland Volcanic Field hazard assessment, GNS Science 2020.")
    lines.append("")

    return "\n".join(lines)


# ---------------------------------------------------------------------------
# Map generation
# ---------------------------------------------------------------------------

def generate_map(gadm_path, output_png):
    try:
        import matplotlib
        matplotlib.use("Agg")
        import matplotlib.pyplot as plt
        from matplotlib.lines import Line2D
    except ImportError:
        print("ERROR: matplotlib not installed. Run: pip install matplotlib")
        sys.exit(1)

    # Load GADM Level 1 boundaries (GeoJSON FeatureCollection)
    with open(gadm_path, "r") as f:
        gadm = json.load(f)

    SKIP_REGIONS = {"ChathamIslands", "NorthernIslands", "SouthernIslands"}

    fig, ax = plt.subplots(figsize=(8, 11), dpi=200)
    fig.patch.set_facecolor("white")
    ax.set_facecolor("#d0e8f5")  # sea colour

    # --- Base layer: GADM region polygons (zorder=1) -----------------------
    for feature in gadm["features"]:
        region_name = feature["properties"].get("NAME_1", "")
        if region_name in SKIP_REGIONS:
            continue

        # MultiPolygon: coordinates → list of polygons → each polygon is
        # list of rings → first ring is exterior → [[lon, lat], ...]
        for polygon in feature["geometry"]["coordinates"]:
            exterior_ring = polygon[0]  # first ring is exterior boundary
            lons = [pt[0] for pt in exterior_ring]
            lats = [pt[1] for pt in exterior_ring]
            ax.fill(
                lons, lats,
                facecolor="#e8e0d0",
                edgecolor="#999999",
                linewidth=0.3,
                zorder=1,
            )

    # --- Mineral deposit points (zorder=5) and labels (zorder=6) -----------
    # Collect legend entries per mineral type
    legend_map = {}
    for lon, lat, color, marker, mtype, label in MAP_POINTS:
        size = 120 if marker == "*" else 80
        ax.scatter(
            lon, lat,
            c=color, marker=marker, s=size,
            edgecolors="black", linewidths=0.6, zorder=5,
        )
        # Offset labels slightly to reduce overlap; use fontsize 5.5
        ax.annotate(
            label, (lon, lat),
            xytext=(5, 4), textcoords="offset points",
            fontsize=5.5, color="#222222", zorder=6,
        )
        if mtype not in legend_map:
            legend_map[mtype] = (color, marker)

    # Build legend handles
    legend_handles = []
    for mtype, (color, marker) in legend_map.items():
        handle = Line2D(
            [0], [0],
            marker=marker, color="w",
            markerfacecolor=color,
            markeredgecolor="black",
            markeredgewidth=0.6,
            markersize=8,
            label=mtype,
        )
        legend_handles.append(handle)

    ax.legend(
        handles=legend_handles,
        loc="lower left",
        fontsize=7,
        framealpha=0.9,
        title="Mineral Type",
        title_fontsize=7.5,
    )

    ax.set_title("NZ Mineral Deposits", fontsize=14, fontweight="bold", pad=10)
    ax.set_xlabel("Longitude (°E)", fontsize=8)
    ax.set_ylabel("Latitude (°S)", fontsize=8)

    # Map limits
    ax.set_xlim(165.5, 179.0)
    ax.set_ylim(-47.5, -34.0)

    # Invert y-axis label sense (latitudes are negative but display positive)
    y_ticks = ax.get_yticks()
    ax.set_yticklabels([f"{abs(v):.0f}" for v in y_ticks], fontsize=7)
    ax.set_xticklabels([f"{v:.0f}" for v in ax.get_xticks()], fontsize=7)

    ax.grid(True, linestyle=":", linewidth=0.4, color="#aaaaaa", alpha=0.7)

    # Source note
    fig.text(
        0.5, 0.01,
        "Boundaries: GADM v4.1 (gadm.org), CC BY 4.0. "
        "Mineral data: GNS Science / Crown Minerals.",
        ha="center", fontsize=6, color="#555555", style="italic",
    )

    plt.tight_layout(rect=[0, 0.02, 1, 1])

    # Ensure output directory exists
    os.makedirs(os.path.dirname(output_png), exist_ok=True)
    plt.savefig(output_png, dpi=200, bbox_inches="tight", facecolor="white")
    plt.close()
    print(f"Map saved to: {output_png}")


# ---------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------

def main():
    script_dir = os.path.dirname(os.path.abspath(__file__))

    # Output paths
    if len(sys.argv) > 1:
        output_md = sys.argv[1]
    else:
        output_md = os.path.join(script_dir, "..", "tables-geological.md")

    if len(sys.argv) > 2:
        output_png = sys.argv[2]
    else:
        output_png = os.path.join(
            script_dir, "..", "site", "images", "geological-minerals.png"
        )

    output_md = os.path.normpath(output_md)
    output_png = os.path.normpath(output_png)

    # GADM NZ Level 1 boundaries (GeoJSON)
    gadm_path = os.path.join(script_dir, "data", "gadm41_NZL_1.json")

    print(f"Generating markdown: {output_md}")
    md_content = generate_markdown()
    os.makedirs(os.path.dirname(output_md), exist_ok=True)
    with open(output_md, "w", encoding="utf-8") as f:
        f.write(md_content)
    print(f"Markdown written: {output_md}")

    print(f"Generating map: {output_png}")
    generate_map(gadm_path, output_png)

    print("Done.")


if __name__ == "__main__":
    main()