From cd082f41a7a5200397deee3dfb3eae91c54cf223 Mon Sep 17 00:00:00 2001 From: Suma Battula Date: Mon, 27 Oct 2025 09:04:07 -0500 Subject: [PATCH 1/5] Add scripts for generating QKrige time series for CONUS(daily and range) --- Scripts/qkrig_ts_daily.py | 70 +++++++++++++++++++++++++++++++++++++++ Scripts/qkrig_ts_range.sh | 34 +++++++++++++++++++ 2 files changed, 104 insertions(+) create mode 100644 Scripts/qkrig_ts_daily.py create mode 100755 Scripts/qkrig_ts_range.sh diff --git a/Scripts/qkrig_ts_daily.py b/Scripts/qkrig_ts_daily.py new file mode 100644 index 0000000..65e6d20 --- /dev/null +++ b/Scripts/qkrig_ts_daily.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python3 +""" +Extract TS for a single date and save as CSV. +Usage: python qkrig_ts_daily.py YYYY-MM-DD /output_dir +""" + +import sys, os, datetime as dt +import numpy as np +import pandas as pd +import geopandas as gpd + +# --- USER CONFIG --- +GPKG_PATH = "~/.ngiab/hydrofabric/v2.2/conus_nextgen.gpkg" +EXPORT_DIR = "/mnt/disk1/usgskrig/exports/gridsize/200/conus/range/100km/" +LAYER = "divides" +ID_FIELD = "divide_id" +GRID_LON_KEY = "grid_lon" +GRID_LAT_KEY = "grid_lat" +GRID_VAL_KEY = "z_interp" + +# --- Load GDF --- +gdf = gpd.read_file(GPKG_PATH, layer=LAYER) +if gdf.crs.is_geographic: + gdf_proj = gdf.to_crs("EPSG:5070") +else: + gdf_proj = gdf +cent_proj = gdf_proj.geometry.centroid +gdf["centroid"] = gpd.GeoSeries(cent_proj, crs=gdf_proj.crs).to_crs(4326) + +# --- Helper functions --- +def npz_path_for_date(d: dt.date) -> str: + return os.path.join(EXPORT_DIR, f"interp_{d.isoformat()}.npz") + +def nearest_grid_value(lons, lats, vals, pt_lon, pt_lat): + ix = np.argmin(np.abs(lons - pt_lon)) + iy = np.argmin(np.abs(lats - pt_lat)) + return float(vals[iy, ix]) + +def grid_sample(npz_path: str, centroids: pd.Series, lon_key: str, lat_key: str, val_key: str) -> pd.Series: + with np.load(npz_path, allow_pickle=True) as z: + L = z[lon_key]; A = z[lat_key]; V = z[val_key] + return centroids.apply(lambda pt: nearest_grid_value(L, A, V, pt.x, pt.y)) + +# --- MAIN --- +if __name__ == "__main__": + if len(sys.argv) < 3: + print("Usage: python extract_ts_csv.py YYYY-MM-DD /output_dir") + sys.exit(1) + + date_str = sys.argv[1] + out_dir = sys.argv[2] + os.makedirs(out_dir, exist_ok=True) + + d = dt.date.fromisoformat(date_str) + npz_file = npz_path_for_date(d) + if not os.path.exists(npz_file): + print(f"No NPZ file for {d}, skipping") + sys.exit(0) + + ser = grid_sample(npz_file, gdf["centroid"], GRID_LON_KEY, GRID_LAT_KEY, GRID_VAL_KEY) + + df_day = pd.DataFrame({ + "date": pd.to_datetime(d), + "divide_id": gdf[ID_FIELD].values, + "value": ser.values + }) + + csv_path = os.path.join(out_dir, f"ts_{date_str}.csv") + df_day.to_csv(csv_path, index=False) + print(f"Saved {csv_path}") diff --git a/Scripts/qkrig_ts_range.sh b/Scripts/qkrig_ts_range.sh new file mode 100755 index 0000000..30529c0 --- /dev/null +++ b/Scripts/qkrig_ts_range.sh @@ -0,0 +1,34 @@ +#!/bin/bash +set -e + +START_DATE="2020-01-01" +END_DATE="2025-08-08" +CSV_DIR="/mnt/disk1/qkrig/conus/" +JOBS=16 # number of parallel jobs + +mkdir -p "$CSV_DIR" + +# Generate list of dates +DATES=$(seq 0 $(( ($(date -d "$END_DATE" +%s) - $(date -d "$START_DATE" +%s)) / 86400 )) | \ + xargs -I{} date -I -d "$START_DATE + {} days") + +# Run extraction in parallel using xargs +echo "$DATES" | xargs -n 1 -P $JOBS -I{} python3 qkrig_ts_daily.py {} "$CSV_DIR" + +echo "All daily CSVs saved in $CSV_DIR" + +# Optional: combine all CSVs into a single file +FINAL_CSV="/mnt/disk1/qkrig/conus/TS_full.csv" +echo "Merging daily CSVs into $FINAL_CSV..." +python3 - < Date: Wed, 5 Nov 2025 13:52:45 -0600 Subject: [PATCH 2/5] Update time series scripts for CONUS hydrofabric --- Scripts/qkrig_ts_daily.py | 58 +++++++++++++++++++++++++-------------- Scripts/qkrig_ts_range.sh | 29 +++++--------------- 2 files changed, 45 insertions(+), 42 deletions(-) mode change 100644 => 100755 Scripts/qkrig_ts_daily.py diff --git a/Scripts/qkrig_ts_daily.py b/Scripts/qkrig_ts_daily.py old mode 100644 new mode 100755 index 65e6d20..51c2bb3 --- a/Scripts/qkrig_ts_daily.py +++ b/Scripts/qkrig_ts_daily.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 """ -Extract TS for a single date and save as CSV. +Extract TS for a single date and save into water-year-wise catchment CSVs. Usage: python qkrig_ts_daily.py YYYY-MM-DD /output_dir """ @@ -10,13 +10,14 @@ import geopandas as gpd # --- USER CONFIG --- -GPKG_PATH = "~/.ngiab/hydrofabric/v2.2/conus_nextgen.gpkg" -EXPORT_DIR = "/mnt/disk1/usgskrig/exports/gridsize/200/conus/range/100km/" -LAYER = "divides" -ID_FIELD = "divide_id" +GPKG_PATH = "~/.ngiab/hydrofabric/v2.2/conus_nextgen.gpkg" +EXPORT_DIR = "/mnt/disk1/usgskrig/exports/gridsize/200/conus/range/100km/all_guages/" +LAYER = "divides" +ID_FIELD = "divide_id" GRID_LON_KEY = "grid_lon" GRID_LAT_KEY = "grid_lat" GRID_VAL_KEY = "z_interp" +GRID_VAR_KEY = "kriging_variance" # variance key confirmed # --- Load GDF --- gdf = gpd.read_file(GPKG_PATH, layer=LAYER) @@ -36,35 +37,52 @@ def nearest_grid_value(lons, lats, vals, pt_lon, pt_lat): iy = np.argmin(np.abs(lats - pt_lat)) return float(vals[iy, ix]) -def grid_sample(npz_path: str, centroids: pd.Series, lon_key: str, lat_key: str, val_key: str) -> pd.Series: +def grid_sample_both(npz_path, centroids): with np.load(npz_path, allow_pickle=True) as z: - L = z[lon_key]; A = z[lat_key]; V = z[val_key] - return centroids.apply(lambda pt: nearest_grid_value(L, A, V, pt.x, pt.y)) + L, A = z[GRID_LON_KEY], z[GRID_LAT_KEY] + V, VV = z[GRID_VAL_KEY], z.get(GRID_VAR_KEY, None) + if VV is None: + return centroids.apply(lambda pt: (nearest_grid_value(L, A, V, pt.x, pt.y), np.nan)) + return centroids.apply(lambda pt: ( + nearest_grid_value(L, A, V, pt.x, pt.y), + nearest_grid_value(L, A, VV, pt.x, pt.y) + )) + +def water_year(d: dt.date) -> int: + """Return the water year for a given date (Oct 1 - Sep 30).""" + return d.year + 1 if d.month >= 10 else d.year # --- MAIN --- if __name__ == "__main__": if len(sys.argv) < 3: - print("Usage: python extract_ts_csv.py YYYY-MM-DD /output_dir") + print("Usage: python qkrig_ts_daily.py YYYY-MM-DD /output_dir") sys.exit(1) date_str = sys.argv[1] out_dir = sys.argv[2] - os.makedirs(out_dir, exist_ok=True) - d = dt.date.fromisoformat(date_str) npz_file = npz_path_for_date(d) + if not os.path.exists(npz_file): print(f"No NPZ file for {d}, skipping") sys.exit(0) - ser = grid_sample(npz_file, gdf["centroid"], GRID_LON_KEY, GRID_LAT_KEY, GRID_VAL_KEY) + # --- Sample daily values (both mean + variance) --- + vals = grid_sample_both(npz_file, gdf["centroid"]) + ser_val = vals.apply(lambda x: x[0]) + ser_var = vals.apply(lambda x: x[1]) + + wy = water_year(d) + wy_dir = os.path.join(out_dir, f"WY{wy}") + os.makedirs(wy_dir, exist_ok=True) - df_day = pd.DataFrame({ - "date": pd.to_datetime(d), - "divide_id": gdf[ID_FIELD].values, - "value": ser.values - }) + # --- Append to per-catchment CSVs --- + for cat_id, val, var_val in zip(gdf[ID_FIELD].values, ser_val.values, ser_var.values): + cat_file = os.path.join(wy_dir, f"{cat_id}.csv") + df_day = pd.DataFrame({"date": [d], "qkrig": [val], "variance": [var_val]}) + if os.path.exists(cat_file): + df_day.to_csv(cat_file, mode='a', header=False, index=False) + else: + df_day.to_csv(cat_file, mode='w', header=True, index=False) - csv_path = os.path.join(out_dir, f"ts_{date_str}.csv") - df_day.to_csv(csv_path, index=False) - print(f"Saved {csv_path}") + print(f"Appended daily values + variance for {date_str} to WY{wy} catchment files") diff --git a/Scripts/qkrig_ts_range.sh b/Scripts/qkrig_ts_range.sh index 30529c0..a9200cb 100755 --- a/Scripts/qkrig_ts_range.sh +++ b/Scripts/qkrig_ts_range.sh @@ -1,34 +1,19 @@ #!/bin/bash set -e -START_DATE="2020-01-01" -END_DATE="2025-08-08" -CSV_DIR="/mnt/disk1/qkrig/conus/" -JOBS=16 # number of parallel jobs +START_DATE="2011-06-01" +END_DATE="2019-09-30" +OUT_DIR="/mnt/disk1/qkrig/conus/" +JOBS=4 # number of parallel jobs -mkdir -p "$CSV_DIR" +mkdir -p "$OUT_DIR" # Generate list of dates DATES=$(seq 0 $(( ($(date -d "$END_DATE" +%s) - $(date -d "$START_DATE" +%s)) / 86400 )) | \ xargs -I{} date -I -d "$START_DATE + {} days") # Run extraction in parallel using xargs -echo "$DATES" | xargs -n 1 -P $JOBS -I{} python3 qkrig_ts_daily.py {} "$CSV_DIR" +echo "$DATES" | xargs -n 1 -P $JOBS -I{} python3 qkrig_ts_daily.py {} "$OUT_DIR" -echo "All daily CSVs saved in $CSV_DIR" +echo "✅ All water-year catchment CSVs updated in $OUT_DIR" -# Optional: combine all CSVs into a single file -FINAL_CSV="/mnt/disk1/qkrig/conus/TS_full.csv" -echo "Merging daily CSVs into $FINAL_CSV..." -python3 - < Date: Wed, 5 Nov 2025 14:50:19 -0600 Subject: [PATCH 3/5] update bash script for streamflow kriging across a range of dates --- Scripts/qkrig_ts_range.sh | 1 + 1 file changed, 1 insertion(+) diff --git a/Scripts/qkrig_ts_range.sh b/Scripts/qkrig_ts_range.sh index a9200cb..d6c813a 100755 --- a/Scripts/qkrig_ts_range.sh +++ b/Scripts/qkrig_ts_range.sh @@ -17,3 +17,4 @@ echo "$DATES" | xargs -n 1 -P $JOBS -I{} python3 qkrig_ts_daily.py {} "$OUT_DIR" echo "✅ All water-year catchment CSVs updated in $OUT_DIR" + From e6dc8ae1ef0d09ca0e954e23761bab95ab21ad4b Mon Sep 17 00:00:00 2001 From: sumabattula Date: Mon, 12 Jan 2026 19:50:37 -0600 Subject: [PATCH 4/5] fixvariance --- src/core/base_krig.py | 78 ++++++++++++++++++++++++++++++++++++++----- 1 file changed, 70 insertions(+), 8 deletions(-) diff --git a/src/core/base_krig.py b/src/core/base_krig.py index 1fc7b7e..6520a8a 100644 --- a/src/core/base_krig.py +++ b/src/core/base_krig.py @@ -6,6 +6,7 @@ from pyproj import Geod from pykrige.ok import OrdinaryKriging from typing import Optional, Tuple, Dict +from scipy.optimize import curve_fit class BaseKrig: def __init__(self, data, config_path, year, month, day): @@ -62,16 +63,77 @@ def __init__(self, data, config_path, year, month, day): # --------------------------------------------------------------------- # Core computations # --------------------------------------------------------------------- + def _spherical_model(self, h, sill, rng, nugget): + """Spherical variogram model.""" + h = np.asarray(h, dtype=float) + gamma = np.where( + h <= rng, + nugget + sill * (1.5 * (h / rng) - 0.5 * (h / rng) ** 3), + nugget + sill, + ) + return gamma + + def fit_variogram_from_empirical(self, bins: Optional[int] = None): + """ + Fit a spherical variogram model to the empirical semivariogram. + Returns sill, range_km, nugget. + """ + if not self.semivariogram_ready(bins): + self.compute_semivariogram(bins=bins) + + h_km, gamma = self._semivar_cache + + mask = np.isfinite(gamma) + h_km = h_km[mask] + gamma = gamma[mask] + + if len(h_km) < 3: + raise RuntimeError("Not enough variogram points to fit model.") + + # --- Initial guesses (robust defaults) + nugget0 = self.config["kriging"].get("nugget", 0.0) + sill0 = np.nanmax(gamma) + range0 = self.config["kriging"].get("range", np.nanmax(h_km)) + + p0 = [sill0, range0, nugget0] + + bounds = ( + (0.0, 1e-6, 0.0), # lower + (np.inf, np.inf, np.inf), # upper + ) + + popt, _ = curve_fit( + self._spherical_model, + h_km, + gamma, + p0=p0, + bounds=bounds, + maxfev=10_000, + ) + + sill, range_km, nugget = popt + return float(sill), float(range_km), float(nugget) + + + + + + + + + def compute_kriging(self): kcfg = self.config.get("kriging", {}) or {} - variogram_params = None - - if kcfg.get("range"): - variogram_params = { - "sill": kcfg.get("sill", None), - "range": float(kcfg["range"]) / 111.0, # km -> degrees (approx) - "nugget": kcfg.get("nugget", 0.0), - } + # --- Fit variogram from empirical data (PER TIMESTEP) + sill, range_km, nugget = self.fit_variogram_from_empirical( + bins=kcfg.get("variogram_bins") + ) + + variogram_params = { + "sill": sill, + "range": range_km / 111.0, # km → degrees + "nugget": nugget, + } ok = OrdinaryKriging( self.lons, self.lats, self.values, From 2555787ad764a758944a3b47ea01d0074b7c6244 Mon Sep 17 00:00:00 2001 From: sumabattula Date: Wed, 14 Jan 2026 15:28:33 -0600 Subject: [PATCH 5/5] daily_variogram_parameters --- src/core/base_krig.py | 150 +++++++++++++++++++++++++----------------- 1 file changed, 88 insertions(+), 62 deletions(-) diff --git a/src/core/base_krig.py b/src/core/base_krig.py index 6520a8a..71c5abe 100644 --- a/src/core/base_krig.py +++ b/src/core/base_krig.py @@ -6,7 +6,6 @@ from pyproj import Geod from pykrige.ok import OrdinaryKriging from typing import Optional, Tuple, Dict -from scipy.optimize import curve_fit class BaseKrig: def __init__(self, data, config_path, year, month, day): @@ -59,84 +58,107 @@ def __init__(self, data, config_path, year, month, day): # Semivariogram cache self._semivar_cache: Optional[Tuple[np.ndarray, np.ndarray]] = None # (bin_centers_km, semi_variance) self._semivar_bins_used: Optional[int] = None - - # --------------------------------------------------------------------- - # Core computations - # --------------------------------------------------------------------- - def _spherical_model(self, h, sill, rng, nugget): - """Spherical variogram model.""" - h = np.asarray(h, dtype=float) - gamma = np.where( - h <= rng, - nugget + sill * (1.5 * (h / rng) - 0.5 * (h / rng) ** 3), - nugget + sill, - ) - return gamma - def fit_variogram_from_empirical(self, bins: Optional[int] = None): + def fit_daily_variogram(self): """ - Fit a spherical variogram model to the empirical semivariogram. - Returns sill, range_km, nugget. + Fit empirical variogram for this day. + Returns sill, nugget, range in DEGREES (PyKrige-ready). """ - if not self.semivariogram_ready(bins): - self.compute_semivariogram(bins=bins) - - h_km, gamma = self._semivar_cache - - mask = np.isfinite(gamma) - h_km = h_km[mask] - gamma = gamma[mask] - - if len(h_km) < 3: - raise RuntimeError("Not enough variogram points to fit model.") + num_points = len(self.lons) + if num_points < 5: + raise ValueError("Too few points for variogram fitting") + + distances = [] + semivariances = [] + + for i in range(num_points): + for j in range(i + 1, num_points): + _, _, d_m = self.geod.inv( + self.lons[i], self.lats[i], + self.lons[j], self.lats[j] + ) + distances.append(d_m / 1000.0) # km + semivariances.append(0.5 * (self.values[i] - self.values[j]) ** 2) + + distances = np.asarray(distances) + semivariances = np.asarray(semivariances) + + n_bins = self.variogram_bins + max_dist = np.percentile(distances, 90) + bins = np.linspace(0, max_dist, n_bins + 1) + bin_ids = np.digitize(distances, bins) + + gamma = [] + bin_centers = [] + + for k in range(1, len(bins)): + mask = bin_ids == k + if np.any(mask): + gamma.append(np.mean(semivariances[mask])) + bin_centers.append(0.5 * (bins[k] + bins[k - 1])) - # --- Initial guesses (robust defaults) - nugget0 = self.config["kriging"].get("nugget", 0.0) - sill0 = np.nanmax(gamma) - range0 = self.config["kriging"].get("range", np.nanmax(h_km)) + gamma = np.asarray(gamma) + bin_centers = np.asarray(bin_centers) - p0 = [sill0, range0, nugget0] + nugget = float(np.min(gamma)) + sill = float(np.percentile(gamma, 90)) - bounds = ( - (0.0, 1e-6, 0.0), # lower - (np.inf, np.inf, np.inf), # upper - ) + # First distance where semivariance reaches sill + idx = np.argmax(gamma >= 0.95 * sill) + range_km = bin_centers[idx] + range_deg = range_km / 111.0 - popt, _ = curve_fit( - self._spherical_model, - h_km, - gamma, - p0=p0, - bounds=bounds, - maxfev=10_000, - ) + return { + "nugget": nugget, + "sill": sill, + "range": range_deg + } - sill, range_km, nugget = popt - return float(sill), float(range_km), float(nugget) - + # --------------------------------------------------------------------- + # Core computations + # --------------------------------------------------------------------- + def compute_kriging(self): + kcfg = self.config.get("kriging", {}) or {} + use_daily = kcfg.get("fit_daily_variogram", False) + variogram_params = None + if use_daily: + try: + vg = self.fit_daily_variogram() - + if vg["range"] <= 0 or not np.isfinite(vg["sill"]): + raise ValueError("Invalid daily variogram") + variogram_params = { + "sill": vg["sill"], + "range": vg["range"], + "nugget": vg["nugget"], + } + print( + f"[{self.year}-{self.month:02d}-{self.day:02d}] " + f"Daily variogram | " + f"sill={vg['sill']:.3f}, " + f"range={vg['range']*111:.1f} km, " + f"nugget={vg['nugget']:.3f}" + ) - def compute_kriging(self): - kcfg = self.config.get("kriging", {}) or {} - # --- Fit variogram from empirical data (PER TIMESTEP) - sill, range_km, nugget = self.fit_variogram_from_empirical( - bins=kcfg.get("variogram_bins") - ) + except Exception as e: + print(f"⚠️ Daily variogram failed, using config values: {e}") - variogram_params = { - "sill": sill, - "range": range_km / 111.0, # km → degrees - "nugget": nugget, - } + if variogram_params is None and kcfg.get("range"): + variogram_params = { + "sill": kcfg.get("sill", None), + "range": float(kcfg["range"]) / 111.0, + "nugget": kcfg.get("nugget", 0.0), + } ok = OrdinaryKriging( - self.lons, self.lats, self.values, + self.lons, + self.lats, + self.values, variogram_model=self.variogram_model, exact_values=kcfg.get("exact_values", True), nlags=kcfg.get("nlags", 12), @@ -144,7 +166,11 @@ def compute_kriging(self): variogram_parameters=variogram_params, ) - self.z_interp, self.kriging_variance = ok.execute("grid", self.grid_lon, self.grid_lat) + self.z_interp, self.kriging_variance = ok.execute( + "grid", self.grid_lon, self.grid_lat + ) + + def compute_semivariogram(self, bins: Optional[int] = None) -> Tuple[np.ndarray, np.ndarray]: """