#!/usr/bin/env python3

# Merge BYU statistical iceberg tracking CSVs and create an EPSG:3031 shapefile.
#
# Processing:
#   - Reads all CSVs from cwd/input/statistical
#   - Adds iceberg_id from each CSV filename
#   - Merges all CSVs while retaining all original attributes
#   - Creates point geometry from the reported lat/lon fields
#   - Writes the full merged CSV to cwd/output/statistical
#   - Writes spatially valid records as EPSG:3031 shapefile
#
# Install requirements:
#   pip install pandas geopandas shapely

from pathlib import Path
import sys

import geopandas as gpd
import pandas as pd
from shapely.geometry import Point


# =============================================================================
# CONFIGURATION
# =============================================================================

BASE_DIR = Path.cwd()

INPUT_DIR = BASE_DIR / "input" / "statistical"
OUTPUT_DIR = BASE_DIR / "output" / "statistical"

OUTPUT_CSV = OUTPUT_DIR / "iceberg_tracking_database_v7-1_statistical.csv"
OUTPUT_SHP = OUTPUT_DIR / "iceberg_tracking_database_v7-1_statistical.shp"

LAT_COL = "lat"
LON_COL = "lon"

SRC_EPSG = 4326
DST_EPSG = 3031


# =============================================================================
# HELPERS
# =============================================================================

def fail(msg: str) -> None:
    print(f"ERROR: {msg}", file=sys.stderr)
    sys.exit(1)


def read_and_tag_csv(csv_path: Path) -> pd.DataFrame:
    df = pd.read_csv(csv_path)
    df.insert(0, "iceberg_id", csv_path.stem)
    return df


# =============================================================================
# SPATIAL OUTPUT
# =============================================================================

def build_geometry(df: pd.DataFrame) -> gpd.GeoDataFrame:
    if LAT_COL not in df.columns or LON_COL not in df.columns:
        fail(f"Required columns '{LAT_COL}' and '{LON_COL}' not found")

    work_df = df.copy()

    work_df[LON_COL] = pd.to_numeric(work_df[LON_COL], errors="coerce")
    work_df[LAT_COL] = pd.to_numeric(work_df[LAT_COL], errors="coerce")

    valid = work_df[LON_COL].notna() & work_df[LAT_COL].notna()

    if valid.sum() == 0:
        fail("No valid coordinates found")

    spatial_df = work_df.loc[valid].copy()

    spatial_df["geometry"] = [
        Point(xy)
        for xy in zip(spatial_df[LON_COL], spatial_df[LAT_COL])
    ]

    gdf = gpd.GeoDataFrame(
        spatial_df,
        geometry="geometry",
        crs=f"EPSG:{SRC_EPSG}",
    )

    return gdf.to_crs(epsg=DST_EPSG)


# =============================================================================
# MAIN
# =============================================================================

def main() -> None:
    if not INPUT_DIR.exists():
        fail(f"Input directory not found: {INPUT_DIR}")

    OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

    csv_files = sorted(INPUT_DIR.glob("*.csv"))
    if not csv_files:
        fail("No CSV files found in input directory")

    print(f"Found {len(csv_files)} CSV files")

    frames = []

    for csv_path in csv_files:
        print(f"Reading {csv_path.name}")
        frames.append(read_and_tag_csv(csv_path))

    merged_df = pd.concat(frames, ignore_index=True, sort=False)

    # Keep iceberg_id as the first column.
    cols = merged_df.columns.tolist()
    cols = ["iceberg_id"] + [c for c in cols if c != "iceberg_id"]
    merged_df = merged_df[cols]

    print("Writing merged CSV")
    merged_df.to_csv(OUTPUT_CSV, index=False)

    print("Creating shapefile")
    gdf = build_geometry(merged_df)

    print("Writing shapefile")
    gdf.to_file(OUTPUT_SHP, driver="ESRI Shapefile")

    print("Done")
    print(f"CSV: {OUTPUT_CSV}")
    print(f"SHP: {OUTPUT_SHP}")


if __name__ == "__main__":
    main()