Script: Harbor Approach Chart Generation

#!/usr/bin/env python3
"""
Fetch harbour bathymetry data from MPI ArcGIS REST API for multiple NZ ports.

The MPI proxy does not support f=geojson, so we fetch f=json (Esri JSON)
and convert the polyline paths to GeoJSON MultiLineString features,
filtering to each port's bounding box.

Output: scripts/data/{port}-bathymetry.geojson

Usage:
    python fetch_bathymetry.py [port|all]

    port: wellington | auckland | lyttelton | tauranga | otago
    Default: all (fetches bathymetry for all ports)

Source:
    MPI Marine Bathymetry ArcGIS REST service, Layer 0 (contour interval 10m < 200m)
    https://maps.mpi.govt.nz/wss/service/ags-relay/arcgis1/gu...
    rest/services/MARINE/MARINE_Bathymetry/MapServer/0
"""

import urllib.request
import json
import os
import ssl
import sys


# ---------------------------------------------------------------------------
# Port bounding boxes (WGS-84) for bathymetry queries
# ---------------------------------------------------------------------------
PORT_BBOXES = {
    "wellington": {
        "name": "Wellington Harbour",
        "xmin": 174.7,  "xmax": 174.95,
        "ymin": -41.35, "ymax": -41.2,
    },
    "auckland": {
        "name": "Auckland (Waitemata Harbour)",
        "xmin": 174.6,  "xmax": 175.0,
        "ymin": -36.95, "ymax": -36.7,
    },
    "lyttelton": {
        "name": "Lyttelton Harbour",
        "xmin": 172.5,  "xmax": 172.85,
        "ymin": -43.7,  "ymax": -43.5,
    },
    "tauranga": {
        "name": "Tauranga Harbour",
        "xmin": 176.1,  "xmax": 176.45,
        "ymin": -37.7,  "ymax": -37.5,
    },
    "otago": {
        "name": "Otago Harbour (Port Chalmers)",
        "xmin": 170.4,  "xmax": 170.8,
        "ymin": -45.85, "ymax": -45.7,
    },
}

# MPI Marine Bathymetry service — Layer 0 (contour interval 10m in depths < 200m)
# Layer 0 provides more contour lines for shallow harbour waters than Layer 2.
BASE_URL = (
    "https://maps.mpi.govt.nz/wss/service/ags-relay/arcgis1/gu...
    "rest/services/MARINE/MARINE_Bathymetry/MapServer/0/query"
)


def fetch_port(port_key, script_dir):
    """Fetch bathymetry for a single port and save as GeoJSON."""
    cfg = PORT_BBOXES[port_key]
    xmin, xmax = cfg["xmin"], cfg["xmax"]
    ymin, ymax = cfg["ymin"], cfg["ymax"]

    output_path = os.path.join(script_dir, "data", f"{port_key}-bathymetry.geojson")

    print(f"\n{'='*60}")
    print(f"Fetching bathymetry for {cfg['name']}")
    print(f"  Bounding box: [{xmin}, {ymin}] to [{xmax}, {ymax}]")
    print(f"{'='*60}")

    # Query parameters — use f=json since the MPI proxy blocks f=geojson
    params = {
        "where": "1=1",
        "geometry": json.dumps({
            "xmin": xmin, "ymin": ymin,
            "xmax": xmax, "ymax": ymax
        }),
        "geometryType": "esriGeometryEnvelope",
        "inSR": "4326",
        "outSR": "4326",
        "outFields": "Depth",
        "f": "json",
        "returnGeometry": "true",
    }

    query_string = "&".join(
        f"{k}={urllib.request.quote(str(v), safe='')}"
        for k, v in params.items()
    )
    url = f"{BASE_URL}?{query_string}"

    print(f"  Querying MPI bathymetry service...")
    print(f"  URL: {url[:120]}...")

    ctx = ssl.create_default_context()
    req = urllib.request.Request(url)
    req.add_header("User-Agent", "RecoveryLibrary/1.0")

    resp = urllib.request.urlopen(req, timeout=120, context=ctx)
    data = resp.read().decode("utf-8")
    obj = json.loads(data)

    if "error" in obj:
        print(f"  ERROR from API: {obj['error']}")
        return False

    if "features" not in obj:
        print(f"  Unexpected response keys: {list(obj.keys())}")
        print(f"  First 500 chars: {data[:500]}")
        return False

    print(f"  Received {len(obj['features'])} raw features")

    # -------------------------------------------------------------------
    # Convert Esri JSON to GeoJSON, filtering paths to port bbox
    # -------------------------------------------------------------------
    geojson = {
        "type": "FeatureCollection",
        "features": []
    }

    for feat in obj["features"]:
        depth = feat["attributes"]["Depth"]
        paths = feat["geometry"]["paths"]

        # Filter paths: keep only segments with at least one point in bbox
        port_paths = []
        for path in paths:
            for pt in path:
                if xmin <= pt[0] <= xmax and ymin <= pt[1] <= ymax:
                    port_paths.append(path)
                    break

        if port_paths:
            feature = {
                "type": "Feature",
                "properties": {"Depth": depth},
                "geometry": {
                    "type": "MultiLineString",
                    "coordinates": port_paths
                }
            }
            geojson["features"].append(feature)
            print(f"  Depth {depth}m: {len(port_paths)} paths in {port_key} area")

    # Save
    os.makedirs(os.path.dirname(output_path), exist_ok=True)
    with open(output_path, "w") as f:
        json.dump(geojson, f)

    print(f"  Saved {len(geojson['features'])} features to {output_path}")
    print(f"  File size: {os.path.getsize(output_path)} bytes")
    return True


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

    # Parse command-line argument
    if len(sys.argv) > 1:
        port_arg = sys.argv[1].lower().strip()
    else:
        port_arg = "all"

    if port_arg == "all":
        ports = list(PORT_BBOXES.keys())
    elif port_arg in PORT_BBOXES:
        ports = [port_arg]
    else:
        print(f"Unknown port: {port_arg}")
        print(f"Valid ports: {', '.join(PORT_BBOXES.keys())}, all")
        sys.exit(1)

    print(f"Fetching bathymetry for: {', '.join(ports)}")

    success = 0
    for port_key in ports:
        try:
            if fetch_port(port_key, script_dir):
                success += 1
        except Exception as e:
            print(f"  FAILED for {port_key}: {e}")

    print(f"\nComplete: {success}/{len(ports)} ports fetched successfully.")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Fetch detailed coastline data from OpenStreetMap Overpass API for NZ harbour charts.

OSM coastline data is MUCH more detailed than GADM for harbour-scale rendering.
This script fetches coastline ways within a bounding box for each port and
converts them to GeoJSON LineString features.

Output: scripts/data/{port}-coastline.geojson

Usage:
    python fetch_osm_coastline.py [port|all]

    port: wellington | auckland | lyttelton | tauranga | otago
    Default: all

Source:
    OpenStreetMap Overpass API — natural=coastline ways
    https://overpass-api.de/api/interpreter
"""

import urllib.request
import json
import os
import ssl
import sys
import time


# ---------------------------------------------------------------------------
# Port bounding boxes (WGS-84) — same as bathymetry fetcher
# ---------------------------------------------------------------------------
PORT_BBOXES = {
    "wellington": {
        "name": "Wellington Harbour",
        "south": -41.40, "west": 174.70,
        "north": -41.20, "east": 174.95,
    },
    "auckland": {
        "name": "Auckland (Waitemata Harbour)",
        "south": -36.95, "west": 174.60,
        "north": -36.70, "east": 175.00,
    },
    "lyttelton": {
        "name": "Lyttelton Harbour",
        "south": -43.70, "west": 172.50,
        "north": -43.50, "east": 172.85,
    },
    "tauranga": {
        "name": "Tauranga Harbour",
        "south": -37.70, "west": 176.10,
        "north": -37.50, "east": 176.45,
    },
    "otago": {
        "name": "Otago Harbour (Port Chalmers)",
        "south": -45.85, "west": 170.40,
        "north": -45.70, "east": 170.80,
    },
}

OVERPASS_URL = "https://overpass-api.de/api/interpreter"


def fetch_coastline(port_key, script_dir):
    """Fetch OSM coastline ways for a single port and save as GeoJSON."""
    cfg = PORT_BBOXES[port_key]
    south, west = cfg["south"], cfg["west"]
    north, east = cfg["north"], cfg["east"]

    output_path = os.path.join(script_dir, "data", f"{port_key}-coastline.geojson")

    print(f"\n{'='*60}")
    print(f"Fetching OSM coastline for {cfg['name']}")
    print(f"  Bounding box: [{west}, {south}] to [{east}, {north}]")
    print(f"{'='*60}")

    # Overpass QL query — fetch coastline ways with full geometry
    query = (
        f'[out:json][bbox:{south},{west},{north},{east}];'
        f'way["natural"="coastline"];'
        f'out geom;'
    )

    # POST the query to the Overpass API
    post_data = f"data={urllib.request.quote(query)}".encode("utf-8")

    print(f"  Querying Overpass API...")
    print(f"  Query: {query}")

    ctx = ssl.create_default_context()
    req = urllib.request.Request(
        OVERPASS_URL,
        data=post_data,
        method="POST",
    )
    req.add_header("User-Agent", "RecoveryLibrary/1.0 (harbour chart generation)")
    req.add_header("Content-Type", "application/x-www-form-urlencoded")

    try:
        resp = urllib.request.urlopen(req, timeout=120, context=ctx)
    except urllib.error.HTTPError as e:
        if e.code == 429:
            print(f"  Rate limited (429). Waiting 30 seconds...")
            time.sleep(30)
            resp = urllib.request.urlopen(req, timeout=120, context=ctx)
        else:
            raise

    data = resp.read().decode("utf-8")
    obj = json.loads(data)

    elements = obj.get("elements", [])
    print(f"  Received {len(elements)} coastline ways")

    if not elements:
        print(f"  WARNING: No coastline ways found in bounding box")
        return False

    # -------------------------------------------------------------------
    # Convert Overpass response to GeoJSON FeatureCollection of LineStrings
    # -------------------------------------------------------------------
    geojson = {
        "type": "FeatureCollection",
        "features": [],
    }

    total_points = 0
    for element in elements:
        if element.get("type") != "way":
            continue
        geometry = element.get("geometry", [])
        if len(geometry) < 2:
            continue

        # Convert {lat, lon} objects to [lon, lat] coordinate pairs
        coordinates = [[pt["lon"], pt["lat"]] for pt in geometry]
        total_points += len(coordinates)

        feature = {
            "type": "Feature",
            "properties": {
                "osm_id": element.get("id"),
                "natural": "coastline",
            },
            "geometry": {
                "type": "LineString",
                "coordinates": coordinates,
            },
        }
        geojson["features"].append(feature)

    print(f"  Converted to {len(geojson['features'])} LineString features")
    print(f"  Total coordinate points: {total_points}")

    # Save
    os.makedirs(os.path.dirname(output_path), exist_ok=True)
    with open(output_path, "w") as f:
        json.dump(geojson, f)

    file_size = os.path.getsize(output_path)
    print(f"  Saved to {output_path}")
    print(f"  File size: {file_size:,} bytes")
    return True


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

    # Parse command-line argument
    if len(sys.argv) > 1:
        port_arg = sys.argv[1].lower().strip()
    else:
        port_arg = "all"

    if port_arg == "all":
        ports = list(PORT_BBOXES.keys())
    elif port_arg in PORT_BBOXES:
        ports = [port_arg]
    else:
        print(f"Unknown port: {port_arg}")
        print(f"Valid ports: {', '.join(PORT_BBOXES.keys())}, all")
        sys.exit(1)

    print(f"Fetching OSM coastline for: {', '.join(ports)}")

    success = 0
    for i, port_key in enumerate(ports):
        try:
            if fetch_coastline(port_key, script_dir):
                success += 1
        except Exception as e:
            print(f"  FAILED for {port_key}: {e}")

        # Rate limiting: pause between ports if fetching multiple
        if i < len(ports) - 1:
            print(f"\n  Pausing 5 seconds (rate limiting)...")
            time.sleep(5)

    print(f"\nComplete: {success}/{len(ports)} ports fetched successfully.")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Generate harbour approach charts for major New Zealand ports.

Produces: site/images/{port}-approach.png for each configured port.

Uses GADM NZ Level 1 coastline data and MPI Marine Bathymetry contours
(depth contour polylines from ArcGIS REST service).

Data preparation:
    1. Run fetch_bathymetry.py to download raw Esri JSON from MPI API
    2. This script loads the pre-processed GeoJSON in scripts/data/

Usage:
    python generate_harbor_charts.py [port]
    port: wellington | auckland | lyttelton | tauranga | otago | all (default)

Run with:
    scripts/.venv/bin/python generate_harbor_charts.py

Sources:
    - MPI Marine Bathymetry ArcGIS REST service (Layer 2: depth contours)
    - GADM v4.1 administrative boundaries (NZL Level 1)
    - LINZ hydrographic charts (navigation features, positions — approximate)
    - NZ Pilot (LINZ) — approach channel descriptions
"""

import json
import os
import sys

# ---------------------------------------------------------------------------
# Depth color mapping: lighter = shallow, darker = deep
# ---------------------------------------------------------------------------
DEPTH_COLORS = {
    10: "#8ecae6",   # light blue — shallow
    20: "#219ebc",   # medium blue
    30: "#126782",   # dark blue
    40: "#023047",   # very dark blue — deep
    50: "#012a4a",   # darkest blue
}

DEPTH_LINEWIDTHS = {
    10: 1.0,
    20: 0.9,
    30: 0.8,
    40: 0.7,
    50: 0.6,
}

# Marker styles for navigation features
MARKER_STYLES = {
    "headland": {"marker": "^", "color": "#2d6a4f", "size": 60},
    "reef":     {"marker": "x", "color": "#d62828", "size": 70},
    "rock":     {"marker": "x", "color": "#d62828", "size": 60},
    "island":   {"marker": "s", "color": "#40916c", "size": 55},
    "port":     {"marker": "D", "color": "#1d3557", "size": 55},
    "beacon":   {"marker": "o", "color": "#e9c46a", "size": 40},
    "channel":  {"marker": "o", "color": "#d62828", "size": 35},
}

# Wider extent for coastline (so we don't clip coastline polygons)
COAST_MARGIN = 0.05

# ---------------------------------------------------------------------------
# Port configurations
# ---------------------------------------------------------------------------
PORT_CONFIGS = {
    "wellington": {
        "title": "WELLINGTON HARBOUR (Port Nicholson)\nApproach Chart — Prototype",
        "chart_ref": "LINZ Chart NZ 4633",
        "chart_extent": {
            "xmin": 174.72, "xmax": 174.94,
            "ymin": -41.38, "ymax": -41.20,
        },
        "coastline_regions": {
            "Wellington", "Marlborough", "Nelson", "Waikato",
            "Manawatu-Wanganui", "Taranaki", "Hawke'sBay",
        },
        "nav_features": [
            ("Pencarrow Head", 174.8553, -41.3648, "headland",
             "Eastern entrance headland; Lt Fl(2) 12s 37m 17M"),
            ("Baring Head", 174.8667, -41.4028, "headland",
             "South of entrance; atmospheric CO2 station"),
            ("Point Halswell", 174.8283, -41.2856, "headland",
             "Inner harbour — Miramar Peninsula"),
            ("Ngauranga", 174.7803, -41.2428, "headland",
             "Northern shore; motorway junction"),
            ("Barrett Reef", 174.8400, -41.3306, "reef",
             "Major danger; TEV Wahine wreck 1968; breaks in heavy swell"),
            ("Steeple Rock", 174.8181, -41.3458, "rock",
             "Isolated rock S of entrance channel; lt Fl.R 4s"),
            ("Matiu / Somes Island", 174.8664, -41.2564, "island",
             "Inner harbour; historic quarantine station; lt Fl 6s"),
            ("Makaro / Ward Island", 174.8572, -41.2728, "island",
             "Small islet S of Matiu"),
            ("CentrePort Wellington", 174.7817, -41.2794, "port",
             "Main commercial port; container terminal; ferry berths"),
            ("Burnham Wharf (RNZN)", 174.7936, -41.3044, "port",
             "Naval base; Shelly Bay area"),
            ("Falcon Shoal", 174.8328, -41.3194, "reef",
             "Shoal near entrance; marks channel eastern boundary"),
        ],
        "label_offsets": {
            "Pencarrow Head": (7, -8),
            "Baring Head": (7, -8),
            "Point Halswell": (-70, 6),
            "Ngauranga": (-55, 6),
            "Barrett Reef": (-65, -5),
            "Steeple Rock": (-70, 5),
            "Matiu / Somes Island": (7, 5),
            "Makaro / Ward Island": (7, -8),
            "CentrePort Wellington": (-100, -10),
            "Burnham Wharf (RNZN)": (-90, -10),
            "Falcon Shoal": (-55, 8),
        },
        "approach_channel": [
            # Source: CentrePort Wellington Pilotage Information
            # Leading light bearing 016.5°; alter 005° at Steeple Beacon;
            # alter 310° at Jerningham/Halswell transit
            (174.835, -41.380),   # Outer approach (chart edge, on 016.5° line)
            (174.845, -41.350),   # Approaching entrance on 016.5° lead
            (174.853, -41.326),   # Waypoint D: abeam Steeple Light
            (174.855, -41.305),   # Mid harbour heading 005°
            (174.858, -41.283),   # Waypoint E: Jerningham/Halswell transit
            (174.842, -41.272),   # Waypoint F: city wharves approach (310°)
            (174.810, -41.276),   # Heading toward CentrePort
            (174.782, -41.279),   # CentrePort berths
        ],
        "geo_labels": [
            ("Wellington\nCity", 174.770, -41.295, 8, "#333333", "italic"),
            ("Evans Bay", 174.810, -41.302, 6, "#555555", "italic"),
            ("Lyall Bay", 174.790, -41.335, 6, "#555555", "italic"),
            ("Oriental\nBay", 174.795, -41.288, 5.5, "#555555", "italic"),
            ("Miramar\nPeninsula", 174.825, -41.305, 6, "#555555", "italic"),
            ("Cook Strait", 174.86, -41.365, 7, "#2a6496", "italic"),
            ("Entrance\nChannel", 174.836, -41.342, 5.5, "#aa3333", "normal"),
            ("Lambton\nHarbour", 174.788, -41.285, 5, "#555555", "italic"),
            ("Hutt Valley", 174.905, -41.215, 6, "#555555", "italic"),
        ],
    },

    "auckland": {
        "title": "AUCKLAND (Waitemata Harbour)\nApproach Chart — Prototype",
        "chart_ref": "LINZ Chart NZ 5322",
        "chart_extent": {
            "xmin": 174.68, "xmax": 174.92,
            "ymin": -36.89, "ymax": -36.71,
        },
        "coastline_regions": {
            "Auckland", "Waikato", "Northland",
        },
        "nav_features": [
            ("North Head", 174.8117, -36.8317, "headland",
             "Entrance headland; Devonport side; historic fortification"),
            ("Bastion Point", 174.8108, -36.8533, "headland",
             "Eastern shore; Tamaki Drive area"),
            ("Rangitoto Island", 174.8597, -36.7894, "island",
             "Volcanic island; prominent landmark 259m"),
            ("Bean Rock", 174.8075, -36.8275, "reef",
             "Reef at harbour entrance; lt Fl(2) 12s"),
            ("Ports of Auckland", 174.7650, -36.8422, "port",
             "Main container port; Bledisloe/Fergusson terminals"),
            ("Devonport Naval Base", 174.7969, -36.8264, "port",
             "RNZN base; North Shore"),
            ("Waiheke Island", 175.0, -36.8, "island",
             "Large island — Hauraki Gulf (label only)"),
        ],
        "label_offsets": {
            "North Head": (7, -8),
            "Bastion Point": (7, -8),
            "Rangitoto Island": (7, 5),
            "Bean Rock": (-55, -8),
            "Ports of Auckland": (-90, -10),
            "Devonport Naval Base": (-95, 8),
            "Waiheke Island": (7, 5),
        },
        "approach_channel": [
            # Source: POAL official passage plan CSV (2024-11)
            # Route: Bravo PS to Fergusson (FX)
            # Downloaded from poal.co.nz/customer-centre/forms-and-resources
            (174.859, -36.720),   # BRAVO PS (pilot boarding, N of Rangitoto)
            (174.802, -36.783),   # ST LEONARD (hdg 215° — W of Rangitoto)
            (174.842, -36.824),   # N LDG LT (North Leading Light, hdg 142°)
            (174.820, -36.834),   # 17 BUOY (hdg 241°)
            (174.808, -36.838),   # 19 BUOY (hdg 245°)
            (174.781, -36.838),   # COMPORT (hdg 270° — due W)
            (174.784, -36.845),   # FX — Fergusson terminal
        ],
        "geo_labels": [
            ("Auckland\nCBD", 174.760, -36.855, 8, "#333333", "italic"),
            ("Devonport", 174.795, -36.835, 6, "#555555", "italic"),
            ("Hauraki\nGulf", 174.875, -36.800, 7, "#2a6496", "italic"),
            ("Rangitoto\nChannel", 174.840, -36.810, 5.5, "#aa3333", "normal"),
            ("Waitemata\nHarbour", 174.800, -36.850, 6, "#2a6496", "italic"),
        ],
    },

    "lyttelton": {
        "title": "LYTTELTON HARBOUR\nApproach Chart — Prototype",
        "chart_ref": "LINZ Chart NZ 6321",
        "chart_extent": {
            "xmin": 172.58, "xmax": 172.88,
            "ymin": -43.66, "ymax": -43.54,
        },
        "coastline_regions": {
            "Canterbury", "West Coast", "Otago",
        },
        "nav_features": [
            ("Godley Head", 172.7742, -43.5853, "headland",
             "Northern entrance headland; lighthouse"),
            ("Adderley Head", 172.7653, -43.6042, "headland",
             "Southern entrance headland"),
            ("Quail Island\n(Otamahua)", 172.6908, -43.6192, "island",
             "Inner harbour island; historic quarantine station"),
            ("Lyttelton Port", 172.7172, -43.6058, "port",
             "Main Canterbury port; coal and container terminal"),
            ("Diamond Harbour", 172.7261, -43.6275, "headland",
             "Southern shore settlement"),
        ],
        "label_offsets": {
            "Godley Head": (7, 5),
            "Adderley Head": (-75, -8),
            "Quail Island\n(Otamahua)": (-75, 5),
            "Lyttelton Port": (-70, -10),
            "Diamond Harbour": (7, -8),
        },
        "approach_channel": [
            # Source: Lyttelton Port Company passage plan
            # PS BRAVO at 43°34.91'S 172°51.22'E, bearing 241° to Camp Bay
            # Camp Bay at 43°36.255'S 172°47.82'E, bearing 261° to Cashin Quay
            (172.854, -43.582),   # PS BRAVO pilot boarding (43°34.91'S 172°51.22'E)
            (172.820, -43.592),   # Between heads, heading 241°
            (172.797, -43.604),   # CAMP BAY waypoint (43°36.255'S 172°47.82'E)
            (172.760, -43.609),   # Mid harbour heading 261°
            (172.728, -43.613),   # CASHIN QUAY (43°36.75'S 172°43.7'E)
        ],
        "geo_labels": [
            ("Lyttelton", 172.720, -43.600, 7, "#333333", "italic"),
            ("Banks\nPeninsula", 172.740, -43.570, 7, "#555555", "italic"),
            ("Canterbury\nBight", 172.620, -43.570, 6, "#2a6496", "italic"),
            ("Harbour\nEntrance", 172.770, -43.595, 5.5, "#aa3333", "normal"),
        ],
    },

    "tauranga": {
        "title": "TAURANGA HARBOUR\nApproach Chart — Prototype",
        "chart_ref": "LINZ Chart NZ 5411",
        "chart_extent": {
            "xmin": 176.12, "xmax": 176.38,
            "ymin": -37.70, "ymax": -37.58,
        },
        "coastline_regions": {
            "BayofPlenty", "Waikato",
        },
        "nav_features": [
            ("Mount Maunganui", 176.1756, -37.6314, "headland",
             "The Mount — 232m; prominent entrance landmark"),
            ("Matakana Island", 176.18, -37.62, "island",
             "Barrier island; northern entrance side"),
            ("Pilot Bay", 176.1708, -37.6389, "headland",
             "Sheltered bay east of The Mount"),
            ("Port of Tauranga", 176.1792, -37.6508, "port",
             "NZ largest export port by volume; container/log terminal"),
            ("Entrance Channel", 176.17, -37.64, "channel",
             "Dredged entrance; strong tidal currents"),
        ],
        "label_offsets": {
            "Mount Maunganui": (7, 5),
            "Matakana Island": (7, -8),
            "Pilot Bay": (-55, -8),
            "Port of Tauranga": (7, -8),
            "Entrance Channel": (-70, 8),
        },
        "approach_channel": [
            # Source: Port of Tauranga — Port Information for Ships' Masters
            # (May 2024) official passage plan waypoints
            (176.198, -37.577),   # PS ZULU (37°34.61'S 176°11.87'E)
            (176.180, -37.602),   # A BEACON (37°36.12'S 176°10.81'E)
            (176.164, -37.626),   # NO.2 REACH (37°37.55'S 176°9.81'E)
            (176.164, -37.637),   # TANEA (37°38.21'S 176°9.81'E)
            (176.179, -37.643),   # ROADS (37°38.59'S 176°10.75'E)
            (176.181, -37.649),   # LINK SPAN (37°38.94'S 176°10.87'E)
            (176.178, -37.658),   # STELLA PASSAGE (37°39.50'S 176°10.69'E)
        ],
        "geo_labels": [
            ("Tauranga", 176.300, -37.680, 8, "#333333", "italic"),
            ("Bay of\nPlenty", 176.200, -37.590, 7, "#2a6496", "italic"),
            ("Matakana\nIsland", 176.200, -37.610, 6, "#555555", "italic"),
        ],
    },

    "otago": {
        "title": "OTAGO HARBOUR (Port Chalmers)\nApproach Chart — Prototype",
        "chart_ref": "LINZ Chart NZ 6612",
        "chart_extent": {
            "xmin": 170.50, "xmax": 170.76,
            "ymin": -45.84, "ymax": -45.72,
        },
        "coastline_regions": {
            "Otago", "Southland", "Canterbury",
        },
        "nav_features": [
            ("Taiaroa Head", 170.7264, -45.7675, "headland",
             "Harbour entrance; albatross colony; lt Fl(2) 30s"),
            ("Aramoana", 170.7000, -45.7725, "headland",
             "Northern entrance headland; mole"),
            ("Heyward Point", 170.6556, -45.7858, "headland",
             "Mid-harbour; narrows area"),
            ("Port Chalmers", 170.6242, -45.8131, "port",
             "Main Otago port; container terminal; cruise berths"),
            ("Goat Island", 170.6444, -45.7917, "island",
             "Mid-harbour island; navigation landmark"),
        ],
        "label_offsets": {
            "Taiaroa Head": (7, 5),
            "Aramoana": (-55, -8),
            "Heyward Point": (-70, 5),
            "Port Chalmers": (-75, -10),
            "Goat Island": (7, -8),
        },
        "approach_channel": [
            # Source: Port Otago Pilotage Guide Rev 7 (September 2022)
            # Inbound PBG to Port Chalmers — official passage plan
            (170.729, -45.727),   # WP001 PBG BRAVO (45°43.61'S 170°43.71'E)
            (170.725, -45.751),   # WP002 Fairway Beacon (45°45.06'S 170°43.48'E)
            (170.724, -45.760),   # WP003 North Approach (45°45.60'S 170°43.46'E)
            (170.722, -45.778),   # WP004 Pilots Beach (45°46.68'S 170°43.30'E)
            (170.721, -45.783),   # WP005 Harrington Point (45°47.00'S 170°43.29'E)
            (170.719, -45.793),   # WP007 Harrington Bend (45°47.57'S 170°43.16'E)
            (170.700, -45.796),   # WP008 Cross Channel East (45°47.76'S 170°42.00'E)
            (170.669, -45.791),   # WP009 Cross Channel West (45°47.45'S 170°40.11'E)
            (170.655, -45.798),   # WP010 Pulling Point (45°47.89'S 170°39.31'E)
            (170.641, -45.799),   # WP011 Acheron Point (45°47.96'S 170°38.49'E)
            (170.633, -45.803),   # WP012 Deborah Bend (45°48.16'S 170°37.96'E)
            (170.630, -45.809),   # WP013 Multi-Purpose N End (45°48.52'S 170°37.79'E)
            (170.631, -45.814),   # WP014 Beach Street N End (45°48.81'S 170°37.88'E)
        ],
        "geo_labels": [
            ("Otago\nHarbour", 170.660, -45.790, 7, "#2a6496", "italic"),
            ("Port\nChalmers", 170.625, -45.820, 7, "#333333", "italic"),
            ("Dunedin", 170.510, -45.830, 7, "#555555", "italic"),
            ("Aramoana\nSpit", 170.690, -45.765, 5.5, "#555555", "italic"),
        ],
    },
}


# ---------------------------------------------------------------------------
# Data loading functions
# ---------------------------------------------------------------------------

def load_coastline(gadm_path, port_key):
    """Load GADM Level 1 coastline polygons for the given port area."""
    cfg = PORT_CONFIGS[port_key]
    extent = cfg["chart_extent"]
    target_regions = cfg["coastline_regions"]

    xmin = extent["xmin"]
    xmax = extent["xmax"]
    ymin = extent["ymin"]
    ymax = extent["ymax"]

    with open(gadm_path, "r") as f:
        gadm = json.load(f)

    polygons = []
    for feature in gadm["features"]:
        region = feature["properties"].get("NAME_1", "")
        if region not in target_regions:
            continue

        geom = feature["geometry"]
        if geom["type"] == "MultiPolygon":
            for polygon in geom["coordinates"]:
                exterior = polygon[0]
                lons = [pt[0] for pt in exterior]
                lats = [pt[1] for pt in exterior]
                if (max(lons) >= xmin - COAST_MARGIN and
                    min(lons) <= xmax + COAST_MARGIN and
                    max(lats) >= ymin - COAST_MARGIN and
                    min(lats) <= ymax + COAST_MARGIN):
                    polygons.append((lons, lats))
        elif geom["type"] == "Polygon":
            exterior = geom["coordinates"][0]
            lons = [pt[0] for pt in exterior]
            lats = [pt[1] for pt in exterior]
            if (max(lons) >= xmin - COAST_MARGIN and
                min(lons) <= xmax + COAST_MARGIN and
                max(lats) >= ymin - COAST_MARGIN and
                min(lats) <= ymax + COAST_MARGIN):
                polygons.append((lons, lats))

    return polygons


def load_osm_coastline(coastline_path):
    """Load OSM coastline GeoJSON for harbour detail overlay."""
    if not os.path.exists(coastline_path):
        return []
    with open(coastline_path, "r") as f:
        gj = json.load(f)
    lines = []
    for feature in gj["features"]:
        geom = feature["geometry"]
        if geom["type"] == "LineString":
            coords = geom["coordinates"]
            lons = [c[0] for c in coords]
            lats = [c[1] for c in coords]
            lines.append((lons, lats))
    return lines


def load_bathymetry(bathy_path):
    """Load bathymetry contours from GeoJSON."""
    with open(bathy_path, "r") as f:
        gj = json.load(f)

    contours = []  # list of (depth, lines)
    for feature in gj["features"]:
        depth = feature["properties"]["Depth"]
        geom = feature["geometry"]
        lines = []
        if geom["type"] == "MultiLineString":
            for line in geom["coordinates"]:
                lines.append(line)
        elif geom["type"] == "LineString":
            lines.append(geom["coordinates"])
        contours.append((depth, lines))

    return contours


def clip_line_to_bbox(line, xmin, ymin, xmax, ymax, margin=0.02):
    """
    Filter a polyline to only include segments where at least one point
    is within the bounding box (with margin). Returns list of sub-lines.
    """
    xmin -= margin
    xmax += margin
    ymin -= margin
    ymax += margin

    result_lines = []
    current_segment = []

    for pt in line:
        in_box = (xmin <= pt[0] <= xmax and ymin <= pt[1] <= ymax)
        if in_box:
            current_segment.append(pt)
        else:
            if len(current_segment) >= 2:
                result_lines.append(current_segment)
            current_segment = []

    if len(current_segment) >= 2:
        result_lines.append(current_segment)

    return result_lines


# ---------------------------------------------------------------------------
# Chart generation
# ---------------------------------------------------------------------------

def generate_chart(port_key, gadm_path, bathy_path, output_png):
    """Generate a harbour approach chart for the given port."""
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    import matplotlib.ticker as ticker
    from matplotlib.lines import Line2D
    import numpy as np

    cfg = PORT_CONFIGS[port_key]
    extent = cfg["chart_extent"]
    chart_xmin = extent["xmin"]
    chart_xmax = extent["xmax"]
    chart_ymin = extent["ymin"]
    chart_ymax = extent["ymax"]

    print(f"\n  Loading coastline for {port_key}...")
    coastline_polys = load_coastline(gadm_path, port_key)
    print(f"    {len(coastline_polys)} polygons in chart area")

    print(f"  Loading bathymetry...")
    contours = load_bathymetry(bathy_path)
    for depth, lines in contours:
        print(f"    Depth {depth}m: {len(lines)} line segments")

    # -------------------------------------------------------------------
    # Create figure
    # -------------------------------------------------------------------
    fig, ax = plt.subplots(figsize=(10, 12), dpi=200)
    fig.patch.set_facecolor("white")

    # Sea background
    ax.set_facecolor("#c8dff0")

    # -------------------------------------------------------------------
    # Layer 1: Coastline polygons (zorder=2)
    # -------------------------------------------------------------------
    for lons, lats in coastline_polys:
        ax.fill(
            lons, lats,
            facecolor="#f5eed6",
            edgecolor="#8a7d5a",
            linewidth=0.3,
            zorder=2,
        )

    # -------------------------------------------------------------------
    # Layer 1b: OSM coastline detail overlay (zorder=2.5)
    # -------------------------------------------------------------------
    coastline_geojson = os.path.join(
        os.path.dirname(os.path.abspath(__file__)),
        "data", f"{port_key}-coastline.geojson"
    )
    osm_lines = load_osm_coastline(coastline_geojson)
    if osm_lines:
        print(f"  Drawing {len(osm_lines)} OSM coastline segments")
    for lons, lats in osm_lines:
        ax.plot(lons, lats, color="#5a4a3a", linewidth=0.8, zorder=2.5,
                solid_capstyle="round")

    # -------------------------------------------------------------------
    # Layer 2: Bathymetry contours (zorder=3)
    # -------------------------------------------------------------------
    label_positions = {}  # depth -> (x, y) for label placement

    for depth, lines in sorted(contours, key=lambda x: x[0]):
        color = DEPTH_COLORS.get(depth, "#666666")
        lw = DEPTH_LINEWIDTHS.get(depth, 0.7)

        for line in lines:
            clipped = clip_line_to_bbox(
                line, chart_xmin, chart_ymin, chart_xmax, chart_ymax
            )
            for seg in clipped:
                xs = [pt[0] for pt in seg]
                ys = [pt[1] for pt in seg]
                ax.plot(xs, ys, color=color, linewidth=lw, zorder=3,
                        solid_capstyle="round")

                seg_len = len(seg)
                if depth not in label_positions or seg_len > label_positions[depth][2]:
                    mid = seg_len // 2
                    label_positions[depth] = (seg[mid][0], seg[mid][1], seg_len)

    # -------------------------------------------------------------------
    # Layer 2b: Depth labels on contours (zorder=4)
    # -------------------------------------------------------------------
    for depth, (lx, ly, _) in label_positions.items():
        if (chart_xmin <= lx <= chart_xmax and chart_ymin <= ly <= chart_ymax):
            ax.annotate(
                f"{depth}m",
                (lx, ly),
                fontsize=6.5,
                fontweight="bold",
                color=DEPTH_COLORS.get(depth, "#333333"),
                backgroundcolor="white",
                ha="center",
                va="center",
                zorder=4,
                bbox=dict(
                    boxstyle="round,pad=0.15",
                    facecolor="white",
                    edgecolor=DEPTH_COLORS.get(depth, "#333333"),
                    linewidth=0.5,
                    alpha=0.85,
                ),
            )

    # -------------------------------------------------------------------
    # Layer 3: Approach channel (zorder=5)
    # -------------------------------------------------------------------
    approach = cfg["approach_channel"]
    ch_lons = [pt[0] for pt in approach]
    ch_lats = [pt[1] for pt in approach]
    ax.plot(
        ch_lons, ch_lats,
        color="#d62828",
        linewidth=1.8,
        linestyle="--",
        zorder=5,
        label="Approach channel",
    )
    # Arrow at end (toward port)
    ax.annotate(
        "",
        xy=(ch_lons[-1], ch_lats[-1]),
        xytext=(ch_lons[-2], ch_lats[-2]),
        arrowprops=dict(
            arrowstyle="->",
            color="#d62828",
            lw=1.8,
        ),
        zorder=5,
    )

    # -------------------------------------------------------------------
    # Layer 4: Navigation features (zorder=6-7)
    # -------------------------------------------------------------------
    nav_features = cfg["nav_features"]
    offsets = cfg.get("label_offsets", {})

    for name, lon, lat, mtype, desc in nav_features:
        style = MARKER_STYLES.get(mtype, {"marker": "o", "color": "#333", "size": 40})

        # Only plot if within chart bounds
        if not (chart_xmin <= lon <= chart_xmax and chart_ymin <= lat <= chart_ymax):
            continue

        scatter_kw = dict(
            c=style["color"],
            marker=style["marker"],
            s=style["size"],
            linewidths=0.5,
            zorder=6,
        )
        # 'x' markers are unfilled — edgecolors not applicable
        if mtype not in ("reef", "rock"):
            scatter_kw["edgecolors"] = "black"
        ax.scatter(lon, lat, **scatter_kw)

        offset = offsets.get(name, (6, 4))
        ax.annotate(
            name,
            (lon, lat),
            xytext=offset,
            textcoords="offset points",
            fontsize=6,
            fontweight="bold",
            color="#222222",
            zorder=7,
            arrowprops=dict(
                arrowstyle="-",
                color="#888888",
                linewidth=0.4,
            ) if abs(offset[0]) > 20 or abs(offset[1]) > 20 else None,
        )

    # -------------------------------------------------------------------
    # Geographic labels
    # -------------------------------------------------------------------
    for text, lon, lat, fsize, color, style in cfg.get("geo_labels", []):
        if (chart_xmin <= lon <= chart_xmax and chart_ymin <= lat <= chart_ymax):
            ax.text(
                lon, lat, text,
                fontsize=fsize,
                color=color,
                fontstyle=style,
                ha="center",
                va="center",
                zorder=8,
            )

    # -------------------------------------------------------------------
    # Chart furniture
    # -------------------------------------------------------------------
    ax.set_xlim(chart_xmin, chart_xmax)
    ax.set_ylim(chart_ymin, chart_ymax)
    ax.set_aspect("auto")

    # Correct aspect ratio for latitude
    mid_lat = (chart_ymin + chart_ymax) / 2
    lat_factor = np.cos(np.radians(abs(mid_lat)))
    ax.set_aspect(1.0 / lat_factor)

    # Title block
    ax.set_title(
        cfg["title"],
        fontsize=14,
        fontweight="bold",
        pad=25,
    )
    # Prominent NOT FOR NAVIGATION warning below title
    fig.text(
        0.5, 0.91,
        "NOT FOR NAVIGATION \u2014 PROTOTYPE ONLY",
        ha="center",
        fontsize=10,
        fontweight="bold",
        color="#cc0000",
        bbox=dict(
            boxstyle="round,pad=0.3",
            facecolor="#fff0f0",
            edgecolor="#cc0000",
            linewidth=1.5,
        ),
    )

    # Axis labels
    ax.set_xlabel("Longitude (\u00b0E)", fontsize=9, labelpad=8)
    ax.set_ylabel("Latitude (\u00b0S)", fontsize=9, labelpad=8)

    # Format tick labels
    ax.tick_params(axis="both", labelsize=7)
    ax.yaxis.set_major_formatter(ticker.FuncFormatter(
        lambda val, pos: f"{abs(val):.2f}\u00b0S"
    ))
    ax.xaxis.set_major_formatter(ticker.FuncFormatter(
        lambda val, pos: f"{val:.2f}\u00b0E"
    ))

    # Grid
    ax.grid(True, linestyle=":", linewidth=0.4, color="#aaaaaa", alpha=0.6, zorder=1)

    # -------------------------------------------------------------------
    # Legend
    # -------------------------------------------------------------------
    legend_elements = [
        Line2D([0], [0], color="#8ecae6", linewidth=2, label="10m contour"),
        Line2D([0], [0], color="#219ebc", linewidth=2, label="20m contour"),
        Line2D([0], [0], color="#126782", linewidth=2, label="30m contour"),
        Line2D([0], [0], color="#023047", linewidth=2, label="40m contour"),
        Line2D([0], [0], color="#d62828", linewidth=1.8, linestyle="--",
               label="Approach channel"),
        Line2D([0], [0], marker="^", color="w", markerfacecolor="#2d6a4f",
               markeredgecolor="black", markersize=8, label="Headland"),
        Line2D([0], [0], marker="x", color="#d62828", markeredgecolor="#d62828",
               markersize=8, linestyle="none", label="Reef / Rock"),
        Line2D([0], [0], marker="s", color="w", markerfacecolor="#40916c",
               markeredgecolor="black", markersize=7, label="Island"),
        Line2D([0], [0], marker="D", color="w", markerfacecolor="#1d3557",
               markeredgecolor="black", markersize=7, label="Port facility"),
    ]

    ax.legend(
        handles=legend_elements,
        loc="lower left",
        fontsize=6.5,
        framealpha=0.92,
        title="Chart Legend",
        title_fontsize=7.5,
        edgecolor="#888888",
    )

    # -------------------------------------------------------------------
    # Compass rose (simple N arrow)
    # -------------------------------------------------------------------
    arrow_x = chart_xmax - 0.018
    arrow_y = chart_ymax + 0.012
    ax.annotate(
        "N",
        xy=(arrow_x, arrow_y + 0.008),
        fontsize=10,
        fontweight="bold",
        ha="center",
        va="bottom",
        zorder=10,
    )
    ax.annotate(
        "",
        xy=(arrow_x, arrow_y + 0.007),
        xytext=(arrow_x, arrow_y - 0.007),
        arrowprops=dict(
            arrowstyle="->",
            color="black",
            lw=2.0,
        ),
        zorder=10,
    )

    # -------------------------------------------------------------------
    # Scale bar (approximate)
    # -------------------------------------------------------------------
    # 1 degree longitude = 111.32 * cos(lat) km
    # 1 nautical mile = 1.852 km
    km_per_deg = 111.32 * np.cos(np.radians(abs(mid_lat)))
    nm_per_deg = km_per_deg / 1.852
    one_nm_deg = 1.0 / nm_per_deg

    # Position scale bar at bottom-center (away from legend in lower-left)
    scale_x = (chart_xmin + chart_xmax) / 2 - 1.5 * one_nm_deg
    scale_y = chart_ymin + 0.008
    for i in range(3):
        color = "black" if i % 2 == 0 else "white"
        ax.plot(
            [scale_x + i * one_nm_deg, scale_x + (i + 1) * one_nm_deg],
            [scale_y, scale_y],
            color=color,
            linewidth=4,
            solid_capstyle="butt",
            zorder=10,
        )
        if color == "white":
            ax.plot(
                [scale_x + i * one_nm_deg, scale_x + (i + 1) * one_nm_deg],
                [scale_y, scale_y],
                color="black",
                linewidth=4.5,
                solid_capstyle="butt",
                zorder=9,
            )
    # Labels
    ax.text(
        scale_x, scale_y - 0.004, "0",
        fontsize=5.5, ha="center", va="top", zorder=10,
    )
    ax.text(
        scale_x + 3 * one_nm_deg, scale_y - 0.004, "3 nm",
        fontsize=5.5, ha="center", va="top", zorder=10,
    )
    ax.text(
        scale_x + 1.5 * one_nm_deg, scale_y + 0.004, "Scale (nautical miles)",
        fontsize=5, ha="center", va="bottom", zorder=10,
    )

    # -------------------------------------------------------------------
    # Source note and caution
    # -------------------------------------------------------------------
    fig.text(
        0.5, 0.005,
        (
            "PROTOTYPE \u2014 NOT FOR NAVIGATION. "
            "Bathymetry: MPI Marine Bathymetry Service. Crown Copyright Reserved. "
            "Coastline: GADM v4.1 (CC BY 4.0) + OpenStreetMap (\u00a9 OSM contributors, ODbL). "
            "Positions approximate. "
            f"Refer to {cfg['chart_ref']} for official navigation."
        ),
        ha="center",
        fontsize=5.5,
        color="#aa0000",
        fontstyle="italic",
    )

    # -------------------------------------------------------------------
    # Save
    # -------------------------------------------------------------------
    plt.tight_layout(rect=[0, 0.02, 1, 0.98])

    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"  Chart saved to: {output_png}")


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

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

    # Parse command-line argument
    if len(sys.argv) > 1:
        port_arg = sys.argv[1].lower().strip()
    else:
        port_arg = "all"

    if port_arg == "all":
        ports = list(PORT_CONFIGS.keys())
    elif port_arg in PORT_CONFIGS:
        ports = [port_arg]
    else:
        print(f"Unknown port: {port_arg}")
        print(f"Valid ports: {', '.join(PORT_CONFIGS.keys())}, all")
        sys.exit(1)

    # Data paths
    gadm_path = os.path.join(script_dir, "data", "gadm41_NZL_1.json")

    if not os.path.exists(gadm_path):
        print(f"ERROR: GADM data not found at {gadm_path}")
        sys.exit(1)

    print(f"Generating charts for: {', '.join(ports)}")

    for port_key in ports:
        bathy_path = os.path.join(script_dir, "data", f"{port_key}-bathymetry.geojson")
        output_png = os.path.join(
            script_dir, "..", "site", "images", f"{port_key}-approach.png"
        )
        output_png = os.path.normpath(output_png)

        if not os.path.exists(bathy_path):
            print(f"\nWARNING: Bathymetry data not found for {port_key} at {bathy_path}")
            print(f"  Run: python fetch_bathymetry.py {port_key}")
            continue

        print(f"\n{'='*60}")
        print(f"Generating chart: {port_key}")
        print(f"  Output: {output_png}")
        print(f"{'='*60}")

        generate_chart(port_key, gadm_path, bathy_path, output_png)

    print("\nDone.")


if __name__ == "__main__":
    main()