"""
kinder_morgan_kmls.py  — KML rectangle generator for Kinder Morgan
Canonical rule:
  - Path >= 50 miles: endpoints extend 1 mile PAST each site (away from other site)
  - Path < 50 miles:  endpoints placed 50 miles inward from each site
  - Rectangle is 2 miles wide (1 mile each side of centerline)
  - Polygon wound CCW for correct Shapely/GeoJSON spatial predicate behavior
"""

import os
import math
import pandas as pd

EARTH_RADIUS_MILES = 3958.8
HALF_WIDTH_MILES   = 1.0
EXTEND_MILES       = 50.0
BUFFER_MILES       = 1.0   # past-site buffer for paths >= 50 miles

CSV_FILE = 'Kinder_Morgan_Network_Definition_v2.csv'
OUT_DIR  = r'D:\KinderMorgan\kml_files'


def haversine_miles(lat1, lon1, lat2, lon2):
    R = EARTH_RADIUS_MILES
    phi1, phi2 = math.radians(lat1), math.radians(lat2)
    dphi  = math.radians(lat2 - lat1)
    dlam  = math.radians(lon2 - lon1)
    a = math.sin(dphi/2)**2 + math.cos(phi1)*math.cos(phi2)*math.sin(dlam/2)**2
    return 2 * R * math.asin(math.sqrt(a))


def initial_bearing(lat1, lon1, lat2, lon2):
    phi1, phi2 = math.radians(lat1), math.radians(lat2)
    dlam = math.radians(lon2 - lon1)
    x = math.sin(dlam) * math.cos(phi2)
    y = math.cos(phi1)*math.sin(phi2) - math.sin(phi1)*math.cos(phi2)*math.cos(dlam)
    return math.degrees(math.atan2(x, y)) % 360


def destination_point(lat, lon, bearing_deg, distance_miles):
    R   = EARTH_RADIUS_MILES
    d   = distance_miles
    b   = math.radians(bearing_deg)
    phi = math.radians(lat)
    lam = math.radians(lon)
    phi2 = math.asin(math.sin(phi)*math.cos(d/R) +
                     math.cos(phi)*math.sin(d/R)*math.cos(b))
    lam2 = lam + math.atan2(math.sin(b)*math.sin(d/R)*math.cos(phi),
                             math.cos(d/R) - math.sin(phi)*math.sin(phi2))
    return math.degrees(phi2), math.degrees(lam2)


def make_kml_rectangle(path_id, tx_lat, tx_lon, rx_lat, rx_lon):
    path_len     = haversine_miles(tx_lat, tx_lon, rx_lat, rx_lon)
    bearing_tx_rx = initial_bearing(tx_lat, tx_lon, rx_lat, rx_lon)
    bearing_rx_tx = initial_bearing(rx_lat, rx_lon, tx_lat, tx_lon)

    if path_len >= EXTEND_MILES:
        # Extend 1 mile past each site away from the other
        tx_end_lat, tx_end_lon = destination_point(tx_lat, tx_lon, bearing_rx_tx, BUFFER_MILES)
        rx_end_lat, rx_end_lon = destination_point(rx_lat, rx_lon, bearing_tx_rx, BUFFER_MILES)
    else:
        # Short path: extend 50 miles inward from each site
        tx_end_lat, tx_end_lon = destination_point(tx_lat, tx_lon, bearing_tx_rx, EXTEND_MILES)
        rx_end_lat, rx_end_lon = destination_point(rx_lat, rx_lon, bearing_rx_tx, EXTEND_MILES)

    # Perpendicular bearings
    perp_right = (bearing_tx_rx + 90)  % 360
    perp_left  = (bearing_tx_rx - 90) % 360

    # Four corners
    A = destination_point(tx_end_lat, tx_end_lon, perp_right, HALF_WIDTH_MILES)
    B = destination_point(tx_end_lat, tx_end_lon, perp_left,  HALF_WIDTH_MILES)
    C = destination_point(rx_end_lat, rx_end_lon, perp_left,  HALF_WIDTH_MILES)
    D = destination_point(rx_end_lat, rx_end_lon, perp_right, HALF_WIDTH_MILES)

    # Wind CCW: A(tx-right) → D(rx-right) → C(rx-left) → B(tx-left) → A
    # This is counter-clockwise when viewed in standard geographic orientation.
    coords = (
        f"{A[1]},{A[0]},0 "
        f"{D[1]},{D[0]},0 "
        f"{C[1]},{C[0]},0 "
        f"{B[1]},{B[0]},0 "
        f"{A[1]},{A[0]},0"
    )

    return f'''<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
  <Document>
    <name>{path_id}</name>
    <Style id="rectStyle">
      <LineStyle>
        <color>ff0000ff</color>
        <width>2</width>
      </LineStyle>
      <PolyStyle>
        <color>400000ff</color>
      </PolyStyle>
    </Style>
    <Placemark>
      <name>{path_id}</name>
      <styleUrl>#rectStyle</styleUrl>
      <Polygon>
        <outerBoundaryIs>
          <LinearRing>
            <coordinates>
              {coords}
            </coordinates>
          </LinearRing>
        </outerBoundaryIs>
      </Polygon>
    </Placemark>
  </Document>
</kml>'''


def main():
    df = pd.read_csv(CSV_FILE)
    os.makedirs(OUT_DIR, exist_ok=True)

    count = 0
    cw_count = 0
    from shapely.geometry import Polygon as ShapelyPolygon
    from xml.etree import ElementTree as ET

    for _, row in df.iterrows():
        path_id = str(row[PATH_ID_COL])
        kml = make_kml_rectangle(
            path_id,
            float(row[LAT_A_COL]), float(row[LON_A_COL]),
            float(row[LAT_B_COL]), float(row[LON_B_COL])
        )
        out_path = os.path.join(OUT_DIR, f"{path_id}.kml")
        with open(out_path, 'w') as f:
            f.write(kml)
        count += 1

        # Inline winding validation
        root = ET.fromstring(kml)
        coords_el = root.find(".//{http://www.opengis.net/kml/2.2}coordinates")
        if coords_el is not None:
            pts = []
            for token in coords_el.text.strip().split():
                parts = token.split(",")
                if len(parts) >= 2:
                    pts.append((float(parts[0]), float(parts[1])))
            if len(pts) >= 3:
                poly = ShapelyPolygon(pts)
                if not poly.exterior.is_ccw:
                    cw_count += 1
                    print(f"  WARNING: {path_id} still CW after fix!")

    print(f"Generated {count} KML files → {OUT_DIR}")
    if cw_count == 0:
        print(f"Winding validation: all {count} polygons are CCW ✓")
    else:
        print(f"Winding validation: {cw_count} polygons are still CW — check geometry logic")


if __name__ == "__main__":
    main()