#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Parsivel monthly TXT → daily NetCDF (Terrace UNBC, MECHE)
- Dependencies: numpy + netCDF4
- Reconstructs 32×32 spectra (V#D# / D#V#) → raw_data[time, D, V]
- Computes fieldN (sum over velocity) and fieldV (velocity-weighted mean)
- Writes extras when present (sensor_temperature, MOR, laser_amplitude, etc.)
- Skips future days automatically (America/Toronto if available)
- Only processes months >= 2023-12
"""

import argparse
import csv
import os
import re
import sys
from datetime import datetime, timezone, date
from pathlib import Path
from typing import Dict, List, Optional, Set, Tuple

import numpy as np
try:
    from netCDF4 import Dataset
except Exception as e:
    sys.stderr.write(f"[fatal] netCDF4 import failed: {e}\n")
    sys.exit(2)

# ============================
# Config (defaults for TERRACE / UNBC Rooftop)
# ============================
IN_DIR   = Path("/instruments/TERRACE/UNBC_Rooftop/parsivel")
PATTERN  = "**/master/*.txt"
OUT_DIR  = Path("/station/instruments/disdrodb/PAR_UNBC_TERRACE")  # id_tag = folder name

# Station metadata (MECHE / Terrace)
SITE_NAME = "UNBC Rooftop"
CAMPAIGN  = "MECHE"
PROJECT   = "Monitoring Extreme Climate and Hydrometeorological Events (MECHE)"
TITLE     = "UNBC Rooftop Disdrometer (OTT Parsivel²)"
DESC      = ("Parsivel² at the Terrace UNBC Weather Station (69 m MSL) measuring "
             "drop size/velocity distributions and precipitation.")
SUMMARY   = "Station at University of Northern British Columbia's Northwest campus in Terrace, BC, Canada."
LOCATION  = "UNBC rooftop weather station, Terrace, BC, Canada"
AUTHORS   = ["Amelie Sauvageau", "Hadleigh Thompson", "Julie M. Thériault", "Stephen Déry"]
CONTACT   = ["Stephen Déry"]
CONTACT_INFO = "stephen.dery@unbc.ca"
INSTITUTION  = "UNBC and UQAM"

# Geo
LAT, LON, ALT_M = 54.5127607, -128.6055111, 69.0

# Processing window: only months >= MIN_YEAR/MIN_MONTH (inclusive)
MIN_YEAR  = 2023
MIN_MONTH = 12

# Speed/behavior flags
SKIP_EXISTING  = True
VERIFY_SAMPLES = False
DEFAULT_COMP   = 4

# Class limits (32×32)
D_MID = np.array([0.062,0.187,0.312,0.437,0.562,0.687,0.812,0.937,1.062,1.187,1.375,1.625,1.875,2.125,2.375,2.750,3.250,3.750,4.250,4.750,5.500,6.500,7.500,8.500,9.500,11.000,13.000,15.000,17.000,19.000,21.500,24.500], dtype=np.float32)
D_SPR = np.array([0.125,0.125,0.125,0.125,0.125,0.125,0.125,0.125,0.125,0.125,0.250,0.250,0.250,0.250,0.250,0.500,0.500,0.500,0.500,0.500,1.000,1.000,1.000,1.000,1.000,2.000,2.000,2.000,2.000,2.000,3.000,3.000], dtype=np.float32)
V_MID = np.array([0.050,0.150,0.250,0.350,0.450,0.550,0.650,0.750,0.850,0.950,1.100,1.300,1.500,1.700,1.900,2.200,2.600,3.000,3.400,3.800,4.400,5.200,6.000,6.800,7.600,8.800,10.400,12.000,13.600,15.200,17.600,20.800], dtype=np.float32)
V_SPR = np.array([0.100,0.100,0.100,0.100,0.100,0.100,0.100,0.100,0.100,0.100,0.200,0.200,0.200,0.200,0.200,0.400,0.400,0.400,0.400,0.400,0.800,0.800,0.800,0.800,0.800,1.600,1.600,1.600,1.600,1.600,3.200,3.200], dtype=np.float32)
D_LOW = np.maximum(0.0, (D_MID - 0.5*D_SPR).astype(np.float32))
D_UP  = (D_MID + 0.5*D_SPR).astype(np.float32)
V_LOW = (V_MID - 0.5*V_SPR).astype(np.float32)
V_UP  = (V_MID + 0.5*V_SPR).astype(np.float32)

TIME_CANDIDATES = ("Timestamp","timestamp","time","DateTime","datetime","Date","UTC","Time_UTC")

# ============================
# Logging helpers
# ============================
def _ts() -> str:
    return datetime.utcnow().strftime("%Y-%m-%dT%H:%M:%SZ")

def info(m: str) -> None:
    print(f"[info  {_ts()}] {m}")

def warn(m: str) -> None:
    print(f"[warn  {_ts()}] {m}")

def error(m: str) -> None:
    print(f"[error {_ts()}] {m}", file=sys.stderr)

# ============================
# "Today" in America/Toronto (fallback: system local)
# ============================
def today_iso_toronto() -> str:
    try:
        from zoneinfo import ZoneInfo  # py>=3.9
        tz = ZoneInfo("America/Toronto")
        return datetime.now(tz).date().isoformat()
    except Exception:
        return date.today().isoformat()

def clamp_days_to_today_str(days: Set[str]) -> Set[str]:
    today = today_iso_toronto()
    return {d for d in days if d < today}

# ============================
# Basics
# ============================
def detect_sep(path: Path) -> str:
    commas = 0
    semis = 0
    with open(path, "r", encoding="utf-8", errors="ignore") as f:
        for _ in range(10):
            line = f.readline()
            if not line:
                break
            tmp = []
            in_q = False
            for ch in line:
                if ch == '"':
                    in_q = not in_q
                    continue
                if not in_q:
                    tmp.append(ch)
            s = ''.join(tmp)
            commas += s.count(",")
            semis  += s.count(";")
    return ";" if semis >= commas else ","

def discover_files(root: Path, pattern: str) -> List[Path]:
    return sorted(root.rglob(pattern))

def parse_utc_timestamp(s: str) -> Optional[datetime]:
    if not s:
        return None
    t = s.strip()
    if t.endswith("Z"):
        t = t[:-1] + "+00:00"
    try:
        dt = datetime.fromisoformat(t)
        if dt.tzinfo is None:
            dt = dt.replace(tzinfo=timezone.utc)
        else:
            dt = dt.astimezone(timezone.utc)
        return dt
    except Exception:
        pass
    fmts = ["%Y-%m-%d %H:%M:%S","%Y/%m/%d %H:%M:%S","%Y-%m-%d %H:%M","%Y-%m-%dT%H:%M:%S","%Y-%m-%dT%H:%M"]
    for fmt in fmts:
        try:
            return datetime.strptime(t, fmt).replace(tzinfo=timezone.utc)
        except Exception:
            continue
    return None

def canon(s: str) -> str:
    return re.sub(r"[^0-9a-z]+", "", s.lower())

def find_col(header: List[str], *aliases: str) -> Optional[int]:
    idx = {canon(c): i for i, c in enumerate(header)}
    for a in aliases:
        ca = canon(a)
        if ca in idx:
            return idx[ca]
    for a in aliases:
        ca = canon(a)
        for i, c in enumerate(header):
            if ca in canon(c):
                return i
    return None

# ============================
# Quick per-day counts
# ============================
def quick_days_and_counts(path: Path) -> Tuple[Dict[str,int], List[str], int]:
    sep = detect_sep(path)
    with open(path, "r", encoding="utf-8", errors="ignore", newline="") as f:
        r = csv.reader(f, delimiter=sep)
        try:
            header = next(r)
        except StopIteration:
            return {}, [], -1
        t_idx = find_col(header, *TIME_CANDIDATES)
        if t_idx is None:
            raise ValueError(f"time column not found in {path.name}")
        counts: Dict[str,int] = {}
        for row in r:
            if not row or t_idx >= len(row):
                continue
            dt = parse_utc_timestamp(row[t_idx])
            if dt is None:
                continue
            key = dt.date().isoformat()
            counts[key] = counts.get(key, 0) + 1
    return counts, header, t_idx

# ============================
# Existing daily scan for a month
# ============================
DATE_RE = re.compile(r".*_(\d{8})\.nc\Z", re.IGNORECASE)

def existing_days_for_month(out_root: Path, year: int, month: int, id_tag: str) -> Set[str]:
    ydir = out_root / f"{year:04d}" / f"{year:04d}{month:02d}"
    out: Set[str] = set()
    if not ydir.is_dir():
        return out
    for fp in ydir.glob(f"{id_tag}_*.nc"):
        m = DATE_RE.match(fp.name)
        if not m:
            continue
        ymd = m.group(1)
        out.add(f"{ymd[:4]}-{ymd[4:6]}-{ymd[6:8]}")
    return out

def month_from_file(fp: Path) -> Optional[Tuple[int,int]]:
    try:
        ym = fp.parent.parent.name  # e.g., "2022_03"
        m = re.match(r"^(\d{4})[_-](\d{2})$", ym)
        if m:
            return int(m.group(1)), int(m.group(2))
    except Exception:
        pass
    return None

# ============================
# Verify existing NetCDF time length
# ============================
def read_time_len_fast(nc_path: Path) -> Optional[int]:
    try:
        with Dataset(str(nc_path), "r") as ds:
            if "time" in ds.dimensions:
                return len(ds.dimensions["time"])
            if "epoch_time" in ds.variables:
                return int(ds.variables["epoch_time"].shape[0])
        return None
    except Exception:
        return None

# ============================
# Planner
# ============================
def plan_days(counts: Dict[str,int],
              out_root: Path,
              id_tag: str,
              fp: Path,
              skip_existing: bool,
              verify_samples: bool,
              forced_days: Set[str]) -> Set[str]:
    today = today_iso_toronto()
    counts = {d:c for d,c in counts.items() if d < today}
    forced_days = {d for d in forced_days if d < today}

    to_do: Set[str] = set()
    ym = month_from_file(fp)
    existing: Set[str] = set()
    if ym:
        y, m = ym
        existing = {d for d in existing_days_for_month(out_root, y, m, id_tag) if d < today}

    missing = set(counts.keys()) - existing

    if skip_existing and not verify_samples and not forced_days and not missing:
        return set()

    to_do.update(missing)

    if verify_samples and ym:
        y, m = ym
        for day in (set(counts.keys()) & existing):
            ymd = day.replace("-", "")
            out_fp = out_root / f"{y:04d}" / f"{y:04d}{m:02d}" / f"{id_tag}_{ymd}.nc"
            have = read_time_len_fast(out_fp)
            if have is None or have != counts[day]:
                to_do.add(day)

    to_do.update(forced_days)

    if skip_existing and not verify_samples:
        to_do = (missing | (forced_days - existing))

    return to_do

# ============================
# NetCDF writer (rich variables)
# ============================
def write_day_netcdf(out_fp: Path,
                     epochs: np.ndarray,
                     t_strings: List[str],
                     raw3d: Optional[np.ndarray],
                     fieldN: Optional[np.ndarray],
                     fieldV: Optional[np.ndarray],
                     extras: Dict[str, np.ndarray],
                     complevel: int,
                     id_tag: str) -> None:
    out_fp.parent.mkdir(parents=True, exist_ok=True)
    with Dataset(str(out_fp), "w", format="NETCDF4") as ds:
        n = int(len(epochs))
        ds.createDimension("time", n)
        ds.createDimension("diameter_classes", 32)
        ds.createDimension("velocity_classes", 32)

        ds.createVariable("latitude",  "f4")[...] = LAT
        ds.createVariable("longitude", "f4")[...] = LON
        ds.createVariable("altitude",  "f4")[...] = ALT_M

        ds.createVariable("lower_diameter_class_limits", "f4", ("diameter_classes",), zlib=True, complevel=complevel)[:] = D_LOW
        ds.createVariable("upper_diameter_class_limits", "f4", ("diameter_classes",), zlib=True, complevel=complevel)[:] = D_UP
        ds.createVariable("lower_velocity_class_limits", "f4", ("velocity_classes",), zlib=True, complevel=complevel)[:] = V_LOW
        ds.createVariable("upper_velocity_class_limits", "f4", ("velocity_classes",), zlib=True, complevel=complevel)[:] = V_UP

        tvar = ds.createVariable("epoch_time", "i4", ("time",), zlib=True, complevel=complevel, shuffle=True)
        tvar[:] = epochs.astype(np.int32, copy=False)

        vlen_str = str
        ts = ds.createVariable("time_as_string", vlen_str, ("time",))
        ts[:] = np.asarray(t_strings, dtype=object)

        if raw3d is not None:
            v = ds.createVariable("raw_data", "f8", ("time","diameter_classes","velocity_classes"),
                                  zlib=True, complevel=complevel, shuffle=True)
            v[:] = raw3d

        if fieldN is not None:
            v = ds.createVariable("fieldN", "f4", ("time","diameter_classes"),
                                  zlib=True, complevel=complevel, shuffle=True)
            v[:] = fieldN.astype(np.float32, copy=False)

        if fieldV is not None:
            v = ds.createVariable("fieldV", "f4", ("time","diameter_classes"),
                                  zlib=True, complevel=complevel, shuffle=True)
            v[:] = fieldV.astype(np.float32, copy=False)

        for name, arr in extras.items():
            if arr.dtype.kind in ("U","O"):
                vv = ds.createVariable(name, vlen_str, ("time",))
                vv[:] = np.asarray(arr, dtype=object)
            elif arr.dtype == np.float32:
                vv = ds.createVariable(name, "f4", ("time",), zlib=True, complevel=complevel, shuffle=True)
                vv[:] = arr
            else:
                vv = ds.createVariable(name, "f8", ("time",), zlib=True, complevel=complevel, shuffle=True)
                vv[:] = arr.astype(np.float64, copy=False)

        # Global attrs (MECHE / Terrace)
        ds.setncattr("title", TITLE)
        ds.setncattr("institution", INSTITUTION)
        ds.setncattr("source", "CANADA")
        ds.setncattr("history", f"produced on {datetime.utcnow().strftime('%Y-%m-%d')} by MECHE pipeline")
        ds.setncattr("sensor_type", "OTT Hydromet Parsivel2 optical disdrometer")
        ds.setncattr("site_name", SITE_NAME)
        ds.setncattr("sensor_name", "PARSIVEL2")
        ds.setncattr("logging_software", "Custom MECHE pipeline")
        ds.setncattr("project_name", PROJECT)
        ds.setncattr("campaign_name", CAMPAIGN)
        ds.setncattr("contributors", ", ".join(AUTHORS))
        ds.setncattr("contact", ", ".join(CONTACT))
        ds.setncattr("contact_information", CONTACT_INFO)
        ds.setncattr("data_license", "CC BY 4.0 https://creativecommons.org/licenses/by/4.0/")
        ds.setncattr("description", DESC)
        ds.setncattr("summary", SUMMARY)
        ds.setncattr("location", LOCATION)
        ds.setncattr("keywords", "Parsivel2, disdrometer, Terrace, UNBC, MECHE, DSD")

# ============================
# Monthly stream → per-day write
# ============================
_NUM_KEEP = re.compile(r"[^\d\.\+\-eE]")
SPEC_RE1 = re.compile(r"^V(\d{1,2})D(\d{1,2})$")
SPEC_RE2 = re.compile(r"^D(\d{1,2})V(\d{1,2})$")

EXTRA_MAP = {
    "rainfall_rate_32bit": ("rainfall_rate_32bit","Intensity of precipitation","Intensity of precipitation (mm/h)"),
    "reflectivity_32bit":  ("reflectivity_32bit","Radar reflectivity","Radar reflectivity (dBz)","Reflectivity"),
    "mor_visibility":      ("mor_visibility","MOR Visibility (m)","Visibility","MOR"),
    "laser_amplitude":     ("laser_amplitude","Signal amplitude of Laserband","Sensor voltage (V)","Laser amplitude","Laser intensity","Laser signal"),
    "number_particles_validated": ("number_particles_validated","Number of detected particles","num_particles"),
    "sensor_temperature":  ("sensor_temperature","Temperature in sensor (C)","Sensor_temperture","Sensor temperature","Temperature sensor","Internal temperature","Temp_sensor","SensorTemp"),
    "kinetic_energy":      ("kinetic_energy","Kinetic Energy"),
    "weather_code_synop_4680": ("weather_code_synop_4680","Weather code SYNOP WaWa"),
    "weather_code_synop_4677": ("weather_code_synop_4677",),
    "weather_code_metar_4678": ("weather_code_metar_4678","Weather code METAR/SPECI","Weather code METAR"),
    "weather_code_nws":    ("weather_code_nws","Weather code NWS"),
    "error_code":          ("error_code","Error code"),
}

def to_float(s: str) -> float:
    if s is None:
        return np.nan
    t = s.strip()
    if t == "" or t.lower() in ("nan","none","null"):
        return np.nan
    if re.search(r"\d,\d", t) and not re.search(r"\d\.\d", t):
        t = t.replace(".", "").replace(",", ".")
    else:
        t = _NUM_KEEP.sub("", t)
    try:
        return float(t)
    except Exception:
        return np.nan

def write_selected_days_from_month(month_fp: Path,
                                   counts: Dict[str,int],
                                   header: List[str],
                                   t_idx: int,
                                   days_to_do: Set[str],
                                   out_root: Path,
                                   complevel: int,
                                   id_tag: str,
                                   write_empty_forced: bool) -> List[Path]:
    written: List[Path] = []
    if not days_to_do:
        return written

    # map extras -> column index + dtype target
    extra_idx: Dict[str, Tuple[int, str]] = {}
    for canon_name, aliases in EXTRA_MAP.items():
        idx = find_col(header, *aliases)
        if idx is not None:
            if canon_name in ("rainfall_rate_32bit","reflectivity_32bit","kinetic_energy"):
                extra_idx[canon_name] = (idx, "f4")
            elif canon_name in ("weather_code_metar_4678","weather_code_nws"):
                extra_idx[canon_name] = (idx, "str")
            else:
                extra_idx[canon_name] = (idx, "f8")

    # detect spectra columns
    spec_cols: List[Tuple[int,int,int]] = []
    for i, raw in enumerate(header):
        if i == t_idx:
            continue
        s = re.sub(r"\s+", "", str(raw))
        m = SPEC_RE1.match(s)
        if m:
            v, d = int(m.group(1)), int(m.group(2))
            if 0 <= d < 32 and 0 <= v < 32:
                spec_cols.append((i,d,v)); continue
        m = SPEC_RE2.match(s)
        if m:
            d, v = int(m.group(1)), int(m.group(2))
            if 0 <= d < 32 and 0 <= v < 32:
                spec_cols.append((i,d,v)); continue

    sep = detect_sep(month_fp)

    # Preallocate per-day buffers
    day_buffers: Dict[str, Dict] = {}
    for day in days_to_do:
        n = int(counts.get(day, 0))
        if n == 0:
            day_buffers[day] = {"n": 0}
            continue
        buf = {
            "pos": 0, "n": n,
            "epochs": np.empty(n, dtype=np.int64),
            "tstr":   np.empty(n, dtype=object),
            "raw":    np.full((n, 32, 32), np.nan, dtype=np.float64) if spec_cols else None,
            "extras": {}
        }
        for name, (idx, dtype) in extra_idx.items():
            if dtype == "str":
                buf["extras"][name] = np.empty(n, dtype=object)
            elif dtype == "f4":
                buf["extras"][name] = np.full(n, np.nan, dtype=np.float32)
            else:
                buf["extras"][name] = np.full(n, np.nan, dtype=np.float64)
        day_buffers[day] = buf

    # Stream rows
    with open(month_fp, "r", encoding="utf-8", errors="ignore", newline="") as f:
        r = csv.reader(f, delimiter=sep)
        try:
            _hdr = next(r)
        except StopIteration:
            _hdr = []
        for row in r:
            if not row or t_idx >= len(row):
                continue
            dt = parse_utc_timestamp(row[t_idx])
            if dt is None:
                continue
            day = dt.date().isoformat()
            if day not in day_buffers:
                continue
            buf = day_buffers[day]
            n = buf.get("n", 0)
            if n == 0:
                continue
            i = buf["pos"]
            if i >= n:
                grow = 256
                new_epochs = np.empty(n + grow, dtype=buf["epochs"].dtype); new_epochs[:n] = buf["epochs"]; buf["epochs"] = new_epochs
                new_tstr   = np.empty(n + grow, dtype=object); new_tstr[:n] = buf["tstr"]; buf["tstr"] = new_tstr
                if buf["raw"] is not None:
                    old = buf["raw"]; new_raw = np.full((old.shape[0] + grow, 32, 32), np.nan, dtype=old.dtype)
                    new_raw[:old.shape[0]] = old; buf["raw"] = new_raw
                for ex, arr in list(buf["extras"].items()):
                    if arr.dtype.kind in ("U","O"):
                        new = np.empty(arr.shape[0] + grow, dtype=object); new[:arr.shape[0]] = arr
                    else:
                        new = np.full(arr.shape[0] + grow, np.nan, dtype=arr.dtype); new[:arr.shape[0]] = arr
                    buf["extras"][ex] = new
                buf["n"] = n + grow; n = buf["n"]

            buf["epochs"][i] = int(dt.timestamp())
            buf["tstr"][i]   = dt.strftime("%Y-%m-%d %H:%M:%S")

            if buf["raw"] is not None:
                grid = buf["raw"][i]
                for (col,d,v) in spec_cols:
                    if col < len(row):
                        # Allow empty → NaN
                        try:
                            grid[d,v] = to_float(row[col])
                        except Exception:
                            grid[d,v] = np.nan

            for name, (col, dtype) in extra_idx.items():
                if col >= len(row): continue
                if dtype == "str":
                    buf["extras"][name][i] = row[col]
                elif dtype == "f4":
                    buf["extras"][name][i] = np.float32(to_float(row[col]))
                else:
                    buf["extras"][name][i] = np.float64(to_float(row[col]))

            buf["pos"] = i + 1

    # Finalize each day and write
    written: List[Path] = []
    for day, buf in day_buffers.items():
        n = int(buf.get("pos", 0)) if buf.get("n", 0) else 0
        ymd = day.replace("-", "")
        y = int(ymd[:4]); m = int(ymd[4:6])
        out_dir = out_root / f"{y:04d}" / f"{y:04d}{m:02d}"
        out_fp = out_dir / f"{id_tag}_{ymd}.nc"

        if n == 0:
            if write_empty_forced and day in days_to_do:
                write_day_netcdf(out_fp, np.array([], dtype=np.int32), [], None, None, None, {}, complevel, id_tag)
                info(f"wrote (empty): {out_fp}"); written.append(out_fp)
            continue

        epochs = buf["epochs"][:n].astype(np.int64, copy=False)
        tstr   = list(buf["tstr"][:n])
        raw3d  = buf["raw"][:n] if buf["raw"] is not None else None

        fieldN = None; fieldV = None
        if raw3d is not None:
            fieldN = np.nansum(raw3d, axis=2).astype(np.float32)
            V = V_MID.astype(np.float64)[None, None, :]
            num = np.nansum(raw3d * V, axis=2)
            den = np.nansum(raw3d, axis=2)
            with np.errstate(invalid="ignore", divide="ignore"):
                fv = (num / den); fv[~np.isfinite(fv)] = np.nan
            fieldV = fv.astype(np.float32)

        extras_out: Dict[str, np.ndarray] = {}
        for name, arr in buf.get("extras", {}).items():
            extras_out[name] = arr[:n]

        write_day_netcdf(out_fp, epochs, tstr, raw3d, fieldN, fieldV, extras_out, complevel, id_tag)
        info(f"wrote: {out_fp}"); written.append(out_fp)

    return written

# ============================
# Metadata YAML (no pyyaml)
# ============================
def write_metadata_yaml(yaml_path: Path, id_tag: str) -> None:
    lines = [
        f"data_source: CANADA",
        f"campaign_name: {CAMPAIGN}",
        f"station_name: {SITE_NAME}",
        f"sensor_name: PARSIVEL2",
        f"reader: LOCAL/MECHE_PARSIVEL2",
        f"raw_data_format: text",
        f"raw_data_glob_pattern: '**/master/*.txt'",
        f"measurement_interval: 60",
        f"deployment_status: ongoing",
        f"deployment_mode: land",
        f"platform_type: fixed",
        f"latitude: {LAT:.7f}",
        f"longitude: {LON:.7f}",
        f"altitude: {ALT_M:.1f}",
        f"title: {TITLE}",
        f"description: {DESC}",
        f"project_name: {PROJECT}",
        "keywords: [Parsivel2, disdrometer, Terrace, UNBC, MECHE, DSD]",
        f"summary: {SUMMARY}",
        f"location: {LOCATION}",
        "country: Canada",
        "continent: North America",
        "sensor_long_name: OTT Parsivel²",
        "sensor_manufacturer: OTT Hydromet",
        "sensor_beam_length: 180",
        "sensor_beam_width: 30",
        f"authors: [{', '.join(AUTHORS)}]",
        f"contact: [{', '.join(CONTACT)}]",
        f"contact_information: {CONTACT_INFO}",
        f"institution: {INSTITUTION}",
        "license: CC-BY-4.0",
        f"id_tag: {id_tag}",
        ""
    ]
    yaml_path.parent.mkdir(parents=True, exist_ok=True)
    yaml_path.write_text("\n".join(lines), encoding="utf-8")
    info(f"wrote metadata: {yaml_path}")

# ============================
# Orchestration
# ============================
def run(
    in_dir: Path,
    out_dir: Path,
    pattern: str,
    skip_existing: bool,
    verify_samples: bool,
    complevel: int,
    force_days: Set[str],
    force_last_n_days: int,
    write_empty_forced: bool,
) -> List[Path]:
    info(f"IN_DIR  = {in_dir}")
    info(f"OUT_DIR = {out_dir}")
    id_tag = out_dir.name
    out_dir.mkdir(parents=True, exist_ok=True)

    files = discover_files(in_dir, pattern)
    if not files:
        warn(f"no files matching {pattern!r} under {in_dir}")
        return []

    # Filter to only months >= 2023-12
    fil2 = []
    for fp in files:
        ym = month_from_file(fp)
        if not ym:
            continue
        y, m = ym
        if (y > MIN_YEAR) or (y == MIN_YEAR and m >= MIN_MONTH):
            fil2.append(fp)
    if not fil2:
        warn(f"no monthly files at or after {MIN_YEAR}-{MIN_MONTH:02d}")
        return []

    created: List[Path] = []
    months_skipped = 0
    today = today_iso_toronto()

    for fp in fil2:
        try:
            counts, header, t_idx = quick_days_and_counts(fp)
        except Exception as e:
            warn(f"quick scan failed: {fp} ({e})")
            continue

        counts = {d:c for d,c in counts.items() if d < today}

        ym = month_from_file(fp)
        forced_here: Set[str] = set()
        if ym:
            y, m = ym
            candidate_forced = {d for d in force_days if (len(d)==10 and int(d[:4])==y and int(d[5:7])==m)}
            forced_here = {d for d in candidate_forced if d < today}
            if force_last_n_days > 0 and counts:
                last_sorted = sorted([d for d in counts.keys() if d < today])[-force_last_n_days:]
                forced_here |= set(last_sorted)

        to_do = plan_days(counts, out_dir, id_tag, fp, skip_existing, verify_samples, forced_here)

        if not to_do:
            months_skipped += 1
            info(f"up-to-date: {fp.name}")
            continue

        try:
            created += write_selected_days_from_month(
                month_fp=fp,
                counts=counts,
                header=header,
                t_idx=t_idx,
                days_to_do=to_do,
                out_root=out_dir,
                complevel=complevel,
                id_tag=id_tag,
                write_empty_forced=write_empty_forced,
            )
        except Exception as e:
            warn(f"write failed for {fp}: {e}")

    info(f"[summary] new_or_updated_days={len(created)} | months_skipped={months_skipped}")
    return created

# ============================
# CLI
# ============================
def parse_args(argv=None):
    p = argparse.ArgumentParser(description="Parsivel monthly TXT → daily NetCDF (Terrace UNBC, MECHE)")
    p.add_argument("--in-dir", type=Path, default=IN_DIR)
    p.add_argument("--out-dir", type=Path, default=OUT_DIR)
    p.add_argument("--pattern", default=PATTERN)
    p.add_argument("--skip-existing", action="store_true", default=SKIP_EXISTING)
    p.add_argument("--no-skip-existing", dest="skip_existing", action="store_false")
    p.add_argument("--verify-samples", action="store_true", default=VERIFY_SAMPLES)
    p.add_argument("--no-verify-samples", dest="verify_samples", action="store_false")
    p.add_argument("--fast-skip", action="store_true", help="Equivalent to --skip-existing --no-verify-samples")
    p.add_argument("--complevel", type=int, default=DEFAULT_COMP)
    p.add_argument("--force-days", type=str, default="", help="Comma-separated YYYY-MM-DD list")
    p.add_argument("--force-last-n-days", type=int, default=0)
    p.add_argument("--write-empty-forced", action="store_true", help="Write placeholder files for forced days with no samples")
    p.add_argument("--no-metadata", dest="write_metadata", action="store_false")
    p.add_argument("--write-metadata", dest="write_metadata", action="store_true", default=True)
    args, unknown = p.parse_known_args(argv)
    if unknown:
        warn(f"ignoring unknown args: {' '.join(unknown)}")
    if args.fast_skip:
        args.skip_existing = True
        args.verify_samples = False
    args.complevel = max(0, min(9, int(args.complevel)))
    fd: Set[str] = set()
    if args.force_days.strip():
        for token in args.force_days.split(","):
            d = token.strip()
            if re.match(r"^\d{4}-\d{2}-\d{2}$", d):
                fd.add(d)
            else:
                warn(f"ignoring invalid --force-days token: {token!r}")
    args.force_days = clamp_days_to_today_str(fd)
    return args

# ============================
# Main
# ============================
def main() -> int:
    py = sys.version_info
    info(f"python: {sys.executable} ({py.major}.{py.minor}.{py.micro})")
    if (py.major, py.minor) < (3,7):
        error("Python >= 3.7 required")
        return 2

    in_ipy = ('ipykernel' in sys.modules) or ('IPython' in sys.modules)
    args = parse_args([] if in_ipy else None)

    created = run(
        in_dir=args.in_dir,
        out_dir=args.out_dir,
        pattern=args.pattern,
        skip_existing=args.skip_existing,
        verify_samples=args.verify_samples,
        complevel=args.complevel,
        force_days=args.force_days,
        force_last_n_days=int(args.force_last_n_days),
        write_empty_forced=bool(args.write_empty_forced),
    )

    if args.write_metadata:
        meta_path = args.out_dir / "MECHE_UNBC_TERRACE_metadata.yml"
        write_metadata_yaml(meta_path, id_tag=args.out_dir.name)

    if created:
        info("==== Created/Updated daily files ====")
        for fp in created:
            print(fp)
    else:
        info("no new or updated daily files were written")
    return 0

if __name__ == "__main__":
    try:
        sys.exit(main())
    except Exception as exc:
        error(str(exc))
        sys.exit(1)