#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Parsivel monthly TXT → daily NetCDF (full variables, minimal deps)

- Dependencies: only numpy + netCDF4 (everything else is stdlib)
- Reconstructs 32×32 spectra from V#D# / D#V# columns → raw_data[time, D, V]
- Computes fieldN (sum over velocity) and fieldV (velocity-weighted mean)
- Writes optional series if present (temperature, visibility, laser amplitude, codes, etc.)
- Skips future days automatically (America/Toronto if available; otherwise system local)
- Fast "already up-to-date" checks (avoid re-writing days)
- CLI flags for force/verify/skip behavior and compression level

Typical runs:
  ./parsivel_daily.py \
    --in-dir /instruments/UQAM_PK/PK_Rooftop/parsivel \
    --out-dir /station/instruments/disdrodb/PAR001_UQAM_PK \
    --fast-skip

  # Rebuild only Aug 2025 folder:
  ./parsivel_daily.py \
    --in-dir /instruments/UQAM_PK/PK_Rooftop/parsivel/2025_08 \
    --pattern 'master/*.txt' \
    --out-dir /station/instruments/disdrodb/PAR001_UQAM_PK \
    --no-skip-existing
"""

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; CLI can override)
# ============================
IN_DIR   = Path("/instruments/UQAM_PK/PK_Rooftop/parsivel")
PATTERN  = "**/master/*.txt"
OUT_DIR  = Path("/station/instruments/disdrodb/PAR001_UQAM_PK")

SITE_NAME = "UQAM-PK Rooftop"
LAT, LON, ALT_M = 45.508667, -73.568889, 69.0

# 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:
        # Python 3.9+: zoneinfo (if present on system)
        from zoneinfo import ZoneInfo  # type: ignore
        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:
    # Robustly detect delimiter by sampling the first ~10 lines
    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
            # ignore commas/semicolons inside quoted strings (roughly)
            # very light heuristic: remove text between pairs of quotes
            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(";")
    # default to semicolon if tied, as many Parsivel exports use ';'
    if semis >= commas:
        return ";"
    return ","


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]:
    """Return (counts_by_YYYY-MM-DD, header, time_col_idx)."""
    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()  # YYYY-MM-DD (UTC date)
            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)  # 'YYYYMMDD'
        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  # YYYY_MM
        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)

        # Variable-length UTF-8 strings for timestamps
        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)

        # Extras (numbers or strings)
        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 attributes
        ds.setncattr("title", "OTT Parsivel2 disdrometer data")
        ds.setncattr("institution", "Université du Québec à Montréal")
        ds.setncattr("source", "CANADA")
        ds.setncattr("history", f"produced on {datetime.utcnow().strftime('%Y-%m-%d')} by UQAM pipeline")
        ds.setncattr("sensor_type", "OTT Hydromet Parsivel2 optical disdrometer")
        ds.setncattr("site_name", SITE_NAME)
        ds.setncattr("sensor_name", id_tag)
        ds.setncattr("logging_software", "Custom UQAM pipeline")
        ds.setncattr("project_name", "UQAM-PK Weather Station")
        ds.setncattr("contributors", "Amelie Sauvageau, Hadleigh Thompson")
        ds.setncattr("data_license", "CC BY 4.0 https://creativecommons.org/licenses/by/4.0/")

# ============================
# 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)"),
    "mor_visibility":      ("mor_visibility","MOR Visibility (m)"),
    "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"),
    "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: list of (col_idx, d, v)
    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 using counts
    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, fill buffers
    with open(month_fp, "r", encoding="utf-8", errors="ignore", newline="") as f:
        r = csv.reader(f, delimiter=sep)
        try:
            _hdr = next(r)  # skip header
        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 buffers in chunks of 256 records if counts underestimated
                grow = 256
                # epochs
                new_epochs = np.empty(n + grow, dtype=buf["epochs"].dtype)
                new_epochs[:n] = buf["epochs"]
                buf["epochs"] = new_epochs
                # tstr
                new_tstr = np.empty(n + grow, dtype=object)
                new_tstr[:n] = buf["tstr"]
                buf["tstr"] = new_tstr
                # raw
                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
                # extras
                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
                        buf["extras"][ex] = new
                    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"]

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

            # fill spectra grid for this row
            if buf["raw"] is not None:
                grid = buf["raw"][i]
                for (col,d,v) in spec_cols:
                    if col < len(row):
                        grid[d,v] = to_float(row[col])

            # fill extras
            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
    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: UQAM-PK",
        f"station_name: {SITE_NAME}",
        f"sensor_name: PARSIVEL2",
        f"reader: LOCAL/UQAM_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:.6f}",
        f"longitude: {LON:.6f}",
        f"altitude: {ALT_M:.1f}",
        f"title: UQAM-PK Rooftop Disdrometer (OTT Parsivel²)",
        f"description: Parsivel² at the UQAM-PK Weather Station (69 m MSL) measuring drop size/velocity distributions and precipitation.",
        f"project_name: UQAM-PK Weather Station",
        "keywords: [Parsivel2, disdrometer, Montreal, UQAM, PK, DSD]",
        "summary: Station at Université du Québec à Montréal, PK building roof.",
        "location: UQAM-PK Weather Station, Montréal, QC, Canada",
        "country: Canada",
        "continent: North America",
        "sensor_long_name: OTT Parsivel²",
        "sensor_manufacturer: OTT Hydromet",
        "sensor_beam_length: 180",
        "sensor_beam_width: 30",
        "authors: [Amelie Sauvageau, Hadleigh Thompson, Julie M. Thériault]",
        "contact: [Amelie Sauvageau]",
        "contact_information: sauvageau.amelie@uqam.ca",
        "institution: UQAM",
        "license: CC-BY-4.0",
        f"id_tag: {id_tag}",
        "disdrodb_data_url: 'https://graph.uqam.ca/instruments/disdrodb/'",
        ""
    ]
    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 []

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

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

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

        # Restrict forced days to this month, and to today-or-earlier
        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}
            dropped = candidate_forced - forced_here
            for d in sorted(dropped):
                warn(f"ignoring future forced day: {d}")
            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 (rich variables, minimal deps)")
    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)))

    # parse force-days into set of 'YYYY-MM-DD', clamp to today
    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)
    dropped = set(fd) - set(args.force_days)
    for d in sorted(dropped):
        warn(f"ignoring future forced day: {d}")

    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 / "UQAM-PK_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)