diff --git a/data/gold/expert_code/harmonize_ess-dive_soilmoisture_data.py b/data/gold/expert_code/harmonize_ess-dive_soilmoisture_data.py deleted file mode 100644 index be4d517..0000000 --- a/data/gold/expert_code/harmonize_ess-dive_soilmoisture_data.py +++ /dev/null @@ -1,1239 +0,0 @@ -# %% -from __future__ import annotations - -import json -import re -from pathlib import Path - -import numpy as np -import pandas as pd -from pyproj import Transformer - -# %% -# ============================================================= -# Config / Paths -# ============================================================= - -HOME = Path.home() -BASE_DIR = Path(HOME, "ess-dive_wfsfa_soil_datasets") -OUT_DIR = Path("data", "processed", "ess-dive_wfsfa_soil_datasets") -MAP_JSON_PATH = OUT_DIR / "sm_data_harmonization_mapping.json" - -OUT_DIR.mkdir(parents=True, exist_ok=True) - -with MAP_JSON_PATH.open("r", encoding="utf-8") as f: - map_json: list[dict] = json.load(f) - -# map_json[0] is the reference dataset (lookup) -REF_IDX = 0 - -# %% -# ============================================================= -# Helpers -# ============================================================= - -def as_list(x): - if x is None: - return [] - return x if isinstance(x, list) else [x] - - -def dsid(idx: int) -> str: - return map_json[idx]["dataset_identifier"] - - -def ds_path(idx: int) -> Path: - return BASE_DIR / dsid(idx) - - -def read_ds_csv(idx: int, filename: str, encoding='utf-8', errors='ignore', **kwargs) -> pd.DataFrame: - return pd.read_csv(ds_path(idx) / filename, encoding=encoding, **kwargs) - - -def parse_local_to_utc(series: pd.Series, fmt: str | None, tz: str) -> pd.Series: - dt = pd.to_datetime(series, format=fmt, errors="coerce") - dt = dt.dt.tz_localize(tz, ambiguous="NaT", nonexistent="shift_forward") - return dt.dt.tz_convert("UTC") - - -def interval_min(s: pd.Series) -> pd.Series: - return s.diff().dt.total_seconds() / 60.0 - - -def utm32613_to_latlon(df: pd.DataFrame, e_col: str, n_col: str) -> pd.DataFrame: - tr = Transformer.from_crs("EPSG:32613", "EPSG:4326", always_xy=True) - e = pd.to_numeric(df[e_col], errors="coerce").values - n = pd.to_numeric(df[n_col], errors="coerce").values - lon, lat = tr.transform(e, n) - out = df.copy() - out["longitude"] = lon - out["latitude"] = lat - return out - - -def ensure_harmonized_cols(df: pd.DataFrame) -> pd.DataFrame: - cols = [ - "datetime_UTC", - "site_id", - "depth_m", - "replicate", - "is_timeseries", - "interval_min", - "volumetric_water_content_m3_m3", - "gravimetric_water_content_gH2O_gs", - "water_potential_kPa", - ] - for c in cols: - if c not in df.columns: - df[c] = np.nan - return df[cols] - - -def add_loc_qc(df: pd.DataFrame) -> pd.DataFrame: - if "qc_flag" not in df.columns: - df["qc_flag"] = np.where(df["latitude"].isna() | df["longitude"].isna(), "g2", np.nan) - return df - - -harmonized_data: list[pd.DataFrame] = [] -harmonized_ids: list[str] = [] -loc_data: list[pd.DataFrame] = [] - -# %% -# ============================================================= -# Dataset 1 -# ============================================================= -idx = 1 -f1 = as_list(map_json[idx]["data_payload_files"]) -m1 = as_list(map_json[idx]["location_metadata_files"])[0] - -ls1 = [read_ds_csv(idx, x) for x in f1] -mdf1 = read_ds_csv(idx, m1) - -ls1_h = [] -for i, d in enumerate(ls1): - cols = [c for c in d.columns if re.search(r"Water_Content|VWC|Potential", c)] - x = d[["Time"] + cols].copy() - - x["datetime_UTC"] = parse_local_to_utc(x["Time"], "%Y-%m-%d %H:%M:%S", "America/Denver") - x["interval_min"] = interval_min(x["datetime_UTC"]) - x = x.drop(columns=["Time"]) - - long = x.melt(id_vars=["datetime_UTC", "interval_min"], var_name="name", value_name="value") - m = long["name"].str.extract(r"^(.*)_at_(.*)$") - long["var_name"] = m[0] - long["depth"] = m[1] - long = long.dropna(subset=["depth"]) - - long["value"] = pd.to_numeric(long["value"], errors="coerce") - long.loc[long["value"] == -9999, "value"] = np.nan - long = long.dropna(subset=["value"]) - - long["dest_var"] = np.where( - long["var_name"].str.contains(r"Water_Content|VWC", regex=True, na=False), - "volumetric_water_content_m3_m3", - np.where( - long["var_name"].str.contains(r"Potential", regex=True, na=False), - "water_potential_kPa", - np.nan, - ), - ) - long = long.dropna(subset=["dest_var"]) - - long["replicate"] = long.groupby(["datetime_UTC", "depth", "dest_var"]).cumcount() + 1 - - wide = ( - long.pivot_table( - index=["datetime_UTC", "interval_min", "depth", "replicate"], - columns="dest_var", - values="value", - aggfunc="first", - ) - .reset_index() - ) - - wide["depth_m"] = pd.to_numeric(wide["depth"].str.replace("cm", "", regex=False), errors="coerce") / 100 - wide["site_id"] = re.sub(r"\.csv$", "", f1[i]) - wide["is_timeseries"] = True - wide["gravimetric_water_content_gH2O_gs"] = np.nan - - ls1_h.append(ensure_harmonized_cols(wide)) - -df1_harmonized = pd.concat(ls1_h, ignore_index=True) -harmonized_data.append(df1_harmonized) -harmonized_ids.append(dsid(idx)) - -loc1_tmp = utm32613_to_latlon(mdf1, "Easting", "Northing") -loc1 = pd.DataFrame( - { - "site_id": loc1_tmp["Name"], - "latitude": loc1_tmp["latitude"], - "longitude": loc1_tmp["longitude"], - } -) -loc1["source_dataset_id"] = dsid(idx) -loc1 = add_loc_qc(loc1) -loc_data.append(loc1) - -# %% -# ============================================================= -# Dataset 2 -# ============================================================= -idx = 2 -f2 = as_list(map_json[idx]["data_payload_files"]) -m2 = as_list(map_json[idx]["location_metadata_files"])[0] - -ls2 = [read_ds_csv(idx, x) for x in f2] -mdf2 = read_ds_csv(idx, m2) - -ls2_h = [] -for i, d in enumerate(ls2): - cols = [c for c in d.columns if re.search(r"Moisture", c)] - x = d[["DateTime"] + cols].copy() - x["datetime_UTC"] = parse_local_to_utc(x["DateTime"], "%m/%d/%Y %I:%M:%S %p", "America/Denver") - x["interval_min"] = interval_min(x["datetime_UTC"]) - x = x.drop(columns=["DateTime"]) - - long = x.melt( - id_vars=["datetime_UTC", "interval_min"], - var_name="name", - value_name="volumetric_water_content_m3_m3", - ) - long["depth_m"] = ( - pd.to_numeric(long["name"].str.split("_").str[-1].str.replace("cm", "", regex=False), errors="coerce") / 100 - ) - long["site_id"] = re.sub(r"\.csv$", "", f2[i]) - long["is_timeseries"] = True - long["water_potential_kPa"] = np.nan - long["gravimetric_water_content_gH2O_gs"] = np.nan - long["replicate"] = long.groupby(["datetime_UTC", "depth_m", "name"]).cumcount() + 1 - - ls2_h.append(ensure_harmonized_cols(long)) - -df2_harmonized = pd.concat(ls2_h, ignore_index=True) -harmonized_data.append(df2_harmonized) -harmonized_ids.append(dsid(idx)) - -loc2 = mdf2.rename(columns={"Name": "site_id", "Lat ": "latitude", "Lon": "longitude"})[ - ["site_id", "latitude", "longitude"] -].copy() -loc2["source_dataset_id"] = dsid(idx) -loc2 = add_loc_qc(loc2) -loc_data.append(loc2) - -# %% -# ============================================================= -# Dataset 3 -# ============================================================= -idx = 3 -f3 = as_list(map_json[idx]["data_payload_files"])[0] -m3 = as_list(map_json[idx]["location_metadata_files"])[0] - -df3 = read_ds_csv(idx, f3) -mdf3 = read_ds_csv(idx, m3) - -x = df3.copy() -x["datetime_UTC"] = parse_local_to_utc(x["TIMESTAMP"], "%Y-%m-%d %H:%M", "America/Denver") -x["interval_min"] = interval_min(x["datetime_UTC"]) - -mp_cols = [c for c in x.columns if "_MP" in c] -for c in mp_cols: - x[c] = pd.to_numeric(x[c], errors="coerce") - -long = x.melt(id_vars=["datetime_UTC", "interval_min"], value_vars=mp_cols, var_name="name", value_name="water_potential_kPa") -long["depth_m"] = pd.to_numeric( - long["name"].str.replace(r"[._]", " ", regex=True).str.split().str[1].str.replace("cm", "", regex=False), - errors="coerce", -) / 100 -long["site_id"] = "ER-LLN1" -long["is_timeseries"] = True -long["water_potential_kPa"] = long["water_potential_kPa"].where(~np.isnan(long["water_potential_kPa"]), np.nan) -long["volumetric_water_content_m3_m3"] = np.nan -long["gravimetric_water_content_gH2O_gs"] = np.nan -long["replicate"] = long.groupby(["datetime_UTC", "depth_m"]).cumcount() + 1 - -df3_harmonized = ensure_harmonized_cols(long) -harmonized_data.append(df3_harmonized) -harmonized_ids.append(dsid(idx)) - -loc3 = mdf3.rename(columns={"Location_ID": "site_id", "Latitude": "latitude", "Longitude": "longitude"})[ - ["site_id", "latitude", "longitude"] -].copy() -loc3["source_dataset_id"] = dsid(idx) -loc3 = add_loc_qc(loc3) -loc_data.append(loc3) - -# %% -# ============================================================= -# Dataset 4 -# ============================================================= -idx = 4 -f4 = as_list(map_json[idx]["data_payload_files"])[0] -m4 = as_list(map_json[idx]["location_metadata_files"])[0] -ddf4 = read_ds_csv(idx, f4) -mdf4 = read_ds_csv(idx, m4) - -x = pd.concat([ddf4.reset_index(drop=True), mdf4.reset_index(drop=True)], axis=1) -dt = pd.to_datetime(x["datetime"], errors="coerce").dt.tz_localize("America/Denver", ambiguous="NaT", nonexistent="shift_forward") -x["datetime_UTC"] = dt.dt.tz_convert("UTC") -x["site_id"] = x["site"] -x["depth_m"] = np.nan -x["replicate"] = 1 -x["is_timeseries"] = True -x["interval_min"] = 24 * 60 -x["volumetric_water_content_m3_m3"] = pd.to_numeric(x["swc"], errors="coerce") -x["water_potential_kPa"] = pd.to_numeric(x["swp"], errors="coerce") -x["gravimetric_water_content_gH2O_gs"] = np.nan -x = x.sort_values(["datetime_UTC", "site_id"]) -x = x[x["site_id"] != "tb"] - -df4_harmonized = ensure_harmonized_cols(x) -harmonized_data.append(df4_harmonized) -harmonized_ids.append(dsid(idx)) - -loc4 = pd.DataFrame( - { - "site_id": ["ph1", "ph2", "sg5"], - "latitude": [38.92, 38.922583, 38.926250], - "longitude": [-106.95, -106.947288, -106.98], - } -) -loc4["source_dataset_id"] = dsid(idx) -loc4 = add_loc_qc(loc4) -loc_data.append(loc4) - -# %% -# ============================================================= -# Dataset 5 -# ============================================================= -idx = 5 -f5 = as_list(map_json[idx]["data_payload_files"])[0] -m5 = as_list(map_json[REF_IDX]["location_metadata_files"])[0] - -ddf5 = read_ds_csv(idx, f5) -mdf5 = read_ds_csv(REF_IDX, m5) - -x = ddf5.copy() -dt_local = pd.to_datetime(x["Date Time"], errors="coerce").dt.tz_localize("America/Denver", ambiguous="NaT", nonexistent="shift_forward") -x["datetime_UTC"] = dt_local.dt.tz_convert("UTC") -x["site_id"] = x["SFA_Location"] -x["depth_m"] = np.nan -x["replicate"] = 1 -x["is_timeseries"] = True -x["interval_min"] = 60 -x["volumetric_water_content_m3_m3"] = pd.to_numeric(x["Measurement"], errors="coerce") -x["water_potential_kPa"] = np.nan -x["gravimetric_water_content_gH2O_gs"] = np.nan - -df5_harmonized = ensure_harmonized_cols(x) -harmonized_data.append(df5_harmonized) -harmonized_ids.append(dsid(idx)) - -loc5 = pd.DataFrame({"site_id": df5_harmonized["site_id"].dropna().unique()}) - -loc5a = loc5.merge( - mdf5[["Location_ID", "Latitude", "Longitude"]], - left_on="site_id", - right_on="Location_ID", - how="left", -).drop(columns=["Location_ID"]) - -loc5b = loc5.merge( - mdf5[["Parent_Location_ID", "Latitude", "Longitude"]], - left_on="site_id", - right_on="Parent_Location_ID", - how="left", -).drop(columns=["Parent_Location_ID"]) - -loc5m = loc5a.merge(loc5b, on="site_id", how="left", suffixes=(".x", ".y")) -loc5m["latitude"] = loc5m["Latitude.x"].combine_first(loc5m["Latitude.y"]) -loc5m["longitude"] = loc5m["Longitude.x"].combine_first(loc5m["Longitude.y"]) -loc5m["source_dataset_id"] = dsid(idx) - -loc5 = loc5m[["site_id", "latitude", "longitude", "source_dataset_id"]].drop_duplicates(subset=["site_id"]) -loc5["qc_flag"] = np.where(loc5["latitude"].isna() | loc5["longitude"].isna(), "g2", "g1") -loc_data.append(loc5) - -# %% -# ============================================================= -# Dataset 6 -# ============================================================= -idx = 6 -f6 = as_list(map_json[idx]["data_payload_files"])[0] -ddf6 = read_ds_csv(idx, f6) - -x = ddf6.iloc[2:].copy() -x["datetime_UTC"] = pd.to_datetime(x["TIMESTAMP"], format="%m/%d/%y %H:%M", errors="coerce", utc=True) -x["site_id"] = x["site"] -x["interval_min"] = interval_min(x["datetime_UTC"]) - -vwc_cols = [c for c in x.columns if "VWC" in c] -long = x.melt( - id_vars=["datetime_UTC", "site_id", "interval_min"], - value_vars=vwc_cols, - var_name="name", - value_name="volumetric_water_content_m3_m3", -) -long["depth_m"] = ( - pd.to_numeric(long["name"].str.split("_").str[-1].str.replace("cm", "", regex=False), errors="coerce") / 100 -) -long["replicate"] = 1 -long["is_timeseries"] = True -long["water_potential_kPa"] = np.nan -long["gravimetric_water_content_gH2O_gs"] = np.nan - -df6_harmonized = ensure_harmonized_cols(long) -harmonized_data.append(df6_harmonized) -harmonized_ids.append(dsid(idx)) - -loc6 = ( - ddf6.iloc[2:] - .groupby("site", as_index=False) - .first()[["site", "latitude", "longitude"]] - .rename(columns={"site": "site_id"}) -) -loc6["source_dataset_id"] = dsid(idx) -loc6 = add_loc_qc(loc6) -loc_data.append(loc6) - -# %% -# ============================================================= -# Dataset 7 -# ============================================================= -idx = 7 -f7 = as_list(map_json[idx]["data_payload_files"])[0] -m7 = as_list(map_json[idx]["location_metadata_files"])[0] - -ddf7 = read_ds_csv(idx, f7) -mdf7 = read_ds_csv(idx, m7, encoding='latin-1') - -x = ddf7.copy() -x["datetime_UTC"] = parse_local_to_utc(x["date.time"], "%m/%d/%y %H:%M", "America/Denver") -x["site_id"] = "BM" -x["interval_min"] = interval_min(x["datetime_UTC"]) -x["depth_m"] = pd.to_numeric(x["Depth (cm)"], errors="coerce") / 100 -x["replicate"] = 1 -x["is_timeseries"] = True -x["volumetric_water_content_m3_m3"] = pd.to_numeric(x["Volumetric Water Content"], errors="coerce") -x["water_potential_kPa"] = np.nan -x["gravimetric_water_content_gH2O_gs"] = np.nan - -df7_harmonized = ensure_harmonized_cols(x) -harmonized_data.append(df7_harmonized) -harmonized_ids.append(dsid(idx)) - -lat_col = [c for c in mdf7.columns if "Latitude" in c][0] -lon_col = [c for c in mdf7.columns if "Longitude" in c][0] -loc7 = mdf7.iloc[[0]].rename(columns={"Location": "site_id", lat_col: "latitude", lon_col: "longitude"})[ - ["site_id", "latitude", "longitude"] -].copy() -loc7["source_dataset_id"] = dsid(idx) -loc7 = add_loc_qc(loc7) -loc_data.append(loc7) - -# %% -# ============================================================= -# Dataset 8 -# ============================================================= -idx = 8 -f8 = as_list(map_json[idx]["data_payload_files"]) -m8 = as_list(map_json[idx]["location_metadata_files"]) - -ls8 = [read_ds_csv(idx, x) for x in f8] -mdf8 = pd.concat([read_ds_csv(idx, x) for x in m8], ignore_index=True) - -ls8_h = [] -for d in ls8: - cols = [c for c in d.columns if re.search(r"SoilMoisture|SoilMatricPot", c)] - x = d[["DateTime.MST"] + cols].copy() - x["datetime_UTC"] = parse_local_to_utc(x["DateTime.MST"], "%Y-%m-%d %H:%M:%S", "America/Denver") - x["interval_min"] = interval_min(x["datetime_UTC"]) - x = x.drop(columns=["DateTime.MST"]) - x.columns = [re.sub(r"\.m3\.m3|\.kPa", "", c) for c in x.columns] - - value_cols = [c for c in x.columns if re.search(r"SoilMoisture|SoilMatricPot", c)] - long = x.melt(id_vars=["datetime_UTC", "interval_min"], value_vars=value_cols, var_name="name", value_name="value") - m = long["name"].str.extract(r"^(SoilMoisture|SoilMatricPot)_(.*)$") - long["var"] = m[0] - long["depth_m"] = pd.to_numeric(m[1].str.replace("cm", "", regex=False), errors="coerce") / 100 - long["value"] = pd.to_numeric(long["value"], errors="coerce") - long.loc[long["value"] == -9999, "value"] = np.nan - - wide = ( - long.pivot_table( - index=["datetime_UTC", "interval_min", "depth_m"], - columns="var", - values="value", - aggfunc="first", - ) - .reset_index() - .rename(columns={"SoilMoisture": "volumetric_water_content_m3_m3", "SoilMatricPot": "water_potential_kPa"}) - ) - - wide["replicate"] = 1 - wide["site_id"] = "Slate River OBJ-2" - wide["is_timeseries"] = True - wide["gravimetric_water_content_gH2O_gs"] = np.nan - - ls8_h.append(ensure_harmonized_cols(wide)) - -df8_harmonized = pd.concat(ls8_h, ignore_index=True) -harmonized_data.append(df8_harmonized) -harmonized_ids.append(dsid(idx)) - -loc8 = mdf8.iloc[[0]].rename(columns={"SiteName": "site_id", "Latitude": "latitude", "Longitude": "longitude"})[ - ["site_id", "latitude", "longitude"] -].copy() -loc8["source_dataset_id"] = dsid(idx) -loc8 = add_loc_qc(loc8) -loc_data.append(loc8) - -# %% -# ============================================================= -# Dataset 9 -# ============================================================= -idx = 9 -f9 = as_list(map_json[idx]["data_payload_files"])[0] -m9 = as_list(map_json[idx]["location_metadata_files"])[0] - -ddf9 = read_ds_csv(idx, f9) -mdf9 = read_ds_csv(idx, m9) - -x = ddf9.copy() -x["datetime_UTC"] = parse_local_to_utc(x["Collection Date"], "%Y-%m-%d", "America/Denver") -x["site_id"] = x["SampleSiteCode"] -x["interval_min"] = np.nan -x["depth_m"] = 0.2 -x["is_timeseries"] = False -x["water_potential_kPa"] = np.nan - -long = x.melt( - id_vars=["datetime_UTC", "site_id", "interval_min", "depth_m", "is_timeseries", "water_potential_kPa"], - value_vars=["VWC_1", "VWC_2"], - var_name="tmp", - value_name="VWC", -) -long["replicate"] = long["tmp"].str.extract(r"_(\d+)")[0] -long["VWC"] = pd.to_numeric(long["VWC"], errors="coerce") -long.loc[long["VWC"] == -9999.0, "VWC"] = np.nan -long = long[long["VWC"].notna()].copy() -long["volumetric_water_content_m3_m3"] = long["VWC"] -long["gravimetric_water_content_gH2O_gs"] = np.nan - -df9_harmonized = ensure_harmonized_cols(long) -harmonized_data.append(df9_harmonized) -harmonized_ids.append(dsid(idx)) - -loc9_tmp = utm32613_to_latlon(mdf9, "Easting", "Northing") -loc9 = pd.DataFrame( - { - "site_id": loc9_tmp["SampleSiteCode"], - "latitude": loc9_tmp["latitude"], - "longitude": loc9_tmp["longitude"], - } -) -loc9["source_dataset_id"] = dsid(idx) -loc9 = add_loc_qc(loc9) -loc_data.append(loc9) - -# %% -# ============================================================= -# Dataset 10 -# ============================================================= -idx = 10 -f10 = as_list(map_json[idx]["data_payload_files"])[0] -m10 = as_list(map_json[REF_IDX]["location_metadata_files"])[0] - -ddf10 = read_ds_csv(idx, f10) -mdf10 = read_ds_csv(REF_IDX, m10) - -x = ddf10.iloc[1:].copy() -x["datetime_UTC"] = parse_local_to_utc(x["Date"], "%Y-%m-%d", "America/Denver") -x["interval_min"] = interval_min(x["datetime_UTC"]) - -vcols = [c for c in x.columns if re.search(r"vol_water_content", c)] -long = x.melt(id_vars=["datetime_UTC", "interval_min"], value_vars=vcols, var_name="name", value_name="volumetric_water_content_m3_m3") -nm = long["name"].str.extract(r"(.*)\._(.*)") -long["site_id"] = nm[0] -long["depth_m"] = np.select( - [long["site_id"].eq("PLM1"), long["site_id"].eq("PLM2"), long["site_id"].eq("PLM3")], - [0.3, 0.28, 0.2], - default=np.nan, -) -long["is_timeseries"] = True -long["water_potential_kPa"] = np.nan -long["replicate"] = 1 -long["gravimetric_water_content_gH2O_gs"] = np.nan - -df10_harmonized = ensure_harmonized_cols(long) -harmonized_data.append(df10_harmonized) -harmonized_ids.append(dsid(idx)) - -sites = df10_harmonized["site_id"].dropna().astype(str).unique().tolist() -pattern = r"(?:^|)(%s)$" % "|".join([re.escape(s) for s in sites]) if sites else r"$^" -loc10 = mdf10[mdf10["Location_ID"].astype(str).str.contains(pattern, regex=True, na=False)].copy() -loc10["site_id"] = loc10["Location_ID"].str.replace("ER-", "", regex=False) -loc10 = loc10.rename(columns={"Latitude": "latitude", "Longitude": "longitude"})[["site_id", "latitude", "longitude"]] -loc10["source_dataset_id"] = dsid(idx) -loc10 = add_loc_qc(loc10) -loc_data.append(loc10) - -# %% -# ============================================================= -# Dataset 11-14 excluded -# ============================================================= - -# %% -# ============================================================= -# Dataset 15 -# ============================================================= -idx = 15 -f15 = as_list(map_json[idx]["data_payload_files"])[0] -m15 = as_list(map_json[idx]["location_metadata_files"])[0] - -ddf15 = read_ds_csv(idx, f15) -mdf15 = read_ds_csv(idx, m15) - -x = ddf15.iloc[1:].copy() -x["datetime_UTC"] = pd.Timestamp("2019-07-02 06:00:00", tz="UTC") -x["is_timeseries"] = False -x["interval_min"] = np.nan - -vcols = [c for c in x.columns if re.search(r"SM\.VWC", c)] -long = x.melt( - id_vars=["datetime_UTC", "GPS_id", "is_timeseries", "interval_min"], - value_vars=vcols, - var_name="name", - value_name="volumetric_water_content_m3_m3", -) -long["replicate"] = long["name"].str.extract(r"\._(.*)$")[0] -long["volumetric_water_content_m3_m3"] = pd.to_numeric(long["volumetric_water_content_m3_m3"], errors="coerce") / 100 -long["depth_m"] = 0.25 -long["water_potential_kPa"] = np.nan -long["gravimetric_water_content_gH2O_gs"] = np.nan -long["site_id"] = long["GPS_id"] - -df15_harmonized = ensure_harmonized_cols(long) -harmonized_data.append(df15_harmonized) -harmonized_ids.append(dsid(idx)) - -loc15_tmp = utm32613_to_latlon(mdf15, "Easting_m", "Northing_m") -loc15 = pd.DataFrame( - { - "site_id": loc15_tmp["GPS_id"], - "latitude": loc15_tmp["latitude"], - "longitude": loc15_tmp["longitude"], - } -) -loc15["source_dataset_id"] = dsid(idx) -loc15 = add_loc_qc(loc15) -loc_data.append(loc15) - -# %% -# ============================================================= -# Dataset 16 -# ============================================================= -idx = 16 -f16 = as_list(map_json[idx]["data_payload_files"]) -m16 = as_list(map_json[idx]["location_metadata_files"])[0] -s16 = as_list(map_json[idx]["sensor_metadata_files"])[0] - -ls16 = [read_ds_csv(idx, x) for x in f16] -mdf16 = read_ds_csv(idx, m16) -sdf16 = read_ds_csv(idx, s16) - -ls16_h = [] -for i, d in enumerate(ls16): - siten = f16[i] - x = d.copy() - x["datetime_UTC"] = parse_local_to_utc(x["TIMESTAMP_END"], "%Y-%m-%d %H:%M:%S", "Etc/GMT+7") - x["interval_min"] = interval_min(x["datetime_UTC"]) - parts = siten.split("_") - x["SITE_ID"] = parts[6] if len(parts) >= 7 else np.nan - - sw_cols = [c for c in x.columns if re.search(r"SWC|SWP", c)] - long = x.melt(id_vars=["datetime_UTC", "interval_min", "SITE_ID"], value_vars=sw_cols, var_name="VARIABLE", value_name="value") - - meta = sdf16.loc[sdf16["VARIABLE"].astype(str).str.contains("SWC|SWP", regex=True, na=False), ["SITE_ID", "VARIABLE", "HEIGHT"]] - long = long.merge(meta, on=["SITE_ID", "VARIABLE"], how="left") - long["VARIABLE"] = long["VARIABLE"].astype(str).str.split("_").str[0] - long["depth_m"] = pd.to_numeric(long["HEIGHT"], errors="coerce") * -1 - long["replicate"] = 1 - long["is_timeseries"] = True - long["gravimetric_water_content_gH2O_gs"] = np.nan - long = long.drop_duplicates() - - wide = ( - long.pivot_table( - index=["datetime_UTC", "SITE_ID", "depth_m", "replicate", "is_timeseries", "interval_min", "gravimetric_water_content_gH2O_gs"], - columns="VARIABLE", - values="value", - aggfunc="first", - ) - .reset_index() - ) - - wide["SWC"] = pd.to_numeric(wide.get("SWC"), errors="coerce") - wide["SWP"] = pd.to_numeric(wide.get("SWP"), errors="coerce") - wide["volumetric_water_content_m3_m3"] = np.where(wide["SWC"].isin([9999.0, -9999.0]), np.nan, wide["SWC"]) - wide["water_potential_kPa"] = np.where(wide["SWP"].isin([9999.0, -9999.0]), np.nan, wide["SWP"]) - wide = wide.rename(columns={"SITE_ID": "site_id"}) - - ls16_h.append(ensure_harmonized_cols(wide)) - -df16_harmonized = pd.concat(ls16_h, ignore_index=True) -harmonized_data.append(df16_harmonized) -harmonized_ids.append(dsid(idx)) - -loc16 = mdf16.rename(columns={"SITE_ID": "site_id", "LOCATION_LAT": "latitude", "LOCATION_LONG": "longitude"})[ - ["site_id", "latitude", "longitude"] -].copy() -loc16["source_dataset_id"] = dsid(idx) -loc16 = add_loc_qc(loc16) -loc_data.append(loc16) - -# %% -# ============================================================= -# Dataset 17 -# ============================================================= -idx = 17 -f17 = as_list(map_json[idx]["data_payload_files"]) -m17 = as_list(map_json[idx]["location_metadata_files"])[0] - -ls17 = [read_ds_csv(idx, x) for x in f17] -mdf17 = read_ds_csv(idx, m17) - -ls17_h = [] -for i, d in enumerate(ls17): - siten = f17[i] - x = d.iloc[1:].copy() - x["datetime_UTC"] = parse_local_to_utc(x["DateTime"], "%Y-%m-%d %H:%M:%S", "America/Denver") - x["interval_min"] = interval_min(x["datetime_UTC"]) - - site_guess = re.split(r"/|\.", siten) - x["site_id"] = site_guess[2] if len(site_guess) > 2 else np.nan - - mc_cols = [c for c in x.columns if re.search(r"MC", c)] - long = x.melt(id_vars=["datetime_UTC", "interval_min", "site_id"], value_vars=mc_cols, var_name="name", value_name="value") - m = long["name"].str.extract(r"(.*)_(.*)") - long["depth_m"] = pd.to_numeric(m[1].str.replace("m", "", regex=False), errors="coerce") - long["volumetric_water_content_m3_m3"] = pd.to_numeric(long["value"], errors="coerce") - long["replicate"] = 1 - long["is_timeseries"] = True - long["water_potential_kPa"] = np.nan - long["gravimetric_water_content_gH2O_gs"] = np.nan - - ls17_h.append(ensure_harmonized_cols(long)) - -df17_harmonized = pd.concat(ls17_h, ignore_index=True) -harmonized_data.append(df17_harmonized) -harmonized_ids.append(dsid(idx)) - -loc17_tmp = utm32613_to_latlon(mdf17, "Easting", "Northing") -loc17 = pd.DataFrame( - { - "site_id": loc17_tmp["ID"], - "latitude": loc17_tmp["latitude"], - "longitude": loc17_tmp["longitude"], - } -) -loc17["source_dataset_id"] = dsid(idx) -loc17 = add_loc_qc(loc17) -loc17 = loc17[loc17["site_id"].astype(str).str.contains("TMC", na=False)] -loc_data.append(loc17) - -# %% -# ============================================================= -# Dataset 18 -# ============================================================= -idx = 18 -f18 = as_list(map_json[idx]["data_payload_files"])[0] -m18 = as_list(map_json[idx]["location_metadata_files"])[0] - -ddf18 = read_ds_csv(idx, f18, encoding="latin1") -mdf18 = read_ds_csv(idx, m18, encoding="latin1") - -x = ddf18.copy() -x["datetime_UTC"] = parse_local_to_utc(x["Date_Collected"], "%m/%d/%Y", "America/Denver") -x["interval_min"] = np.nan -x["site_id"] = (x["Field_Site"].astype(str) + "_" + x["Plot"].astype(str) + "_" + x["Topographic_Position"].astype(str)).str.replace(r"_$", "", regex=True) -x["depth_m"] = np.select( - [ - x["Depth_Increment"].eq("0-5 cm"), - x["Depth_Increment"].isin(["5-15 cm", "15-May"]), - x["Depth_Increment"].eq("15 cm +"), - ], - [0.025, 0.10, 0.15], - default=np.nan, -) -x["is_timeseries"] = False -x["water_potential_kPa"] = np.nan -x["volumetric_water_content_m3_m3"] = np.nan -x["gravimetric_water_content_gH2O_gs"] = pd.to_numeric(x["Soil Water Content (g H2O per gram soil)"], errors="coerce") -x["replicate"] = x["Replicate"] - -df18_harmonized = ensure_harmonized_cols(x) -harmonized_data.append(df18_harmonized) -harmonized_ids.append(dsid(idx)) - -loc18 = mdf18.copy() -loc18["site_id"] = (loc18["Field_Site"].astype(str) + "_" + loc18["Plot"].astype(str) + "_" + loc18["Topographic_Position"].astype(str)).str.replace(r"_$", "", regex=True) -loc18["latitude"] = pd.to_numeric(loc18["Latitude"].astype(str).str.extract(r"([-+]?\d*\.?\d+)")[0], errors="coerce") -loc18["longitude"] = pd.to_numeric(loc18["Longitude"].astype(str).str.extract(r"([-+]?\d*\.?\d+)")[0], errors="coerce") -loc18["source_dataset_id"] = dsid(idx) -loc18 = loc18[["site_id", "latitude", "longitude", "source_dataset_id"]].drop_duplicates() -loc18 = add_loc_qc(loc18) -loc_data.append(loc18) - -# %% -# ============================================================= -# Dataset 19-22 excluded -# ============================================================= - -# %% -# ============================================================= -# Dataset 23 -# ============================================================= -idx = 23 -m23 = as_list(map_json[REF_IDX]["location_metadata_files"])[0] - -ddf23 = read_ds_csv(idx, "WM_SWC.csv") -sdf23 = read_ds_csv(idx, "sensor_metadata.csv") -pdf23 = read_ds_csv(idx, "plot_metadata.csv") -mdf23 = read_ds_csv(REF_IDX, m23) - -x = ddf23[ddf23["Sensor.Type"] == "SWC"].copy() -x["datetime_UTC"] = parse_local_to_utc(x["Date Time"], "%Y-%m-%d %H:%M:%S", "America/Denver") - -smeta = sdf23[sdf23["Sensor Type"] == "SWC"].copy() -smeta["Sensor.SN"] = pd.to_numeric(smeta["Sensor.SN"], errors="coerce") -smeta = smeta[["Sensor.SN", "Depth.cm"]] - -x["Sensor.SN"] = pd.to_numeric(x["Sensor.SN"], errors="coerce") -x = x.merge(smeta, on="Sensor.SN", how="left") -x = x.merge(pdf23, on="Plot.Location", how="left") -x = x[x["Treatment"] == "control"].copy() - -x["site_id"] = x["Point.Location"].astype(str) + "_" + x["Veg"].astype(str) -x["depth_m"] = pd.to_numeric(x["Depth.cm"], errors="coerce") / 100 -x["is_timeseries"] = True -x["volumetric_water_content_m3_m3"] = pd.to_numeric(x["Measurement"], errors="coerce") -x["water_potential_kPa"] = np.nan -x["gravimetric_water_content_gH2O_gs"] = np.nan - -x["rep_key"] = x["site_id"].astype(str) + "|" + x["Sensor.SN"].astype(str) + "|" + x["depth_m"].astype(str) -x["replicate"] = pd.factorize(x["rep_key"])[0] + 1 - -x = x.sort_values(["site_id", "depth_m", "replicate", "datetime_UTC"]) -x["interval_min"] = x.groupby(["site_id", "depth_m", "replicate"])["datetime_UTC"].diff().dt.total_seconds() / 60.0 -x.loc[x["interval_min"] < 0, "interval_min"] = np.nan - -df23_harmonized = ensure_harmonized_cols(x) -harmonized_data.append(df23_harmonized) -harmonized_ids.append(dsid(idx)) - -loc23 = pd.DataFrame({"site_id": df23_harmonized["site_id"].dropna().unique()}) -loc23["latitude"] = np.nan -loc23["longitude"] = np.nan -loc23 = loc23[loc23["site_id"].isin(df23_harmonized["site_id"].dropna().unique())].copy() -loc23["source_dataset_id"] = dsid(idx) -loc23 = add_loc_qc(loc23).drop_duplicates() -loc_data.append(loc23) - -loc23 = mdf23.rename(columns={"Location_ID": "site_id", "Latitude": "latitude", "Longitude": "longitude"})[ - ["site_id", "latitude", "longitude"] -].copy() -loc23 = loc23[loc23["site_id"].isin(df23_harmonized["site_id"].str.split('_').str[0].dropna().unique())].copy() -loc23["source_dataset_id"] = dsid(idx) -loc23 = add_loc_qc(loc23) -loc_data.append(loc23) - - -# %% -# ============================================================= -# Dataset 24 -# ============================================================= -idx = 24 -f24 = as_list(map_json[idx]["data_payload_files"])[0] -m24 = as_list(map_json[idx]["location_metadata_files"])[0] - -ddf24 = read_ds_csv(idx, f24, sep=";") -mdf24 = read_ds_csv(idx, m24, sep=";") - -x = ddf24.copy() -keep = ["Timestamp"] + [c for c in x.columns if re.search(r"^P[12]_MP_(15|30|60)$", c)] -x = x[keep] -x["datetime_UTC"] = parse_local_to_utc(x["Timestamp"], "%Y-%m-%d", "America/Denver") -x["interval_min"] = interval_min(x["datetime_UTC"]) -x = x.drop(columns=["Timestamp"]) - -mp_cols = [c for c in x.columns if re.search(r"^P[12]_MP_(15|30|60)$", c)] -long = x.melt(id_vars=["datetime_UTC", "interval_min"], value_vars=mp_cols, var_name="name", value_name="MP") -m = long["name"].str.extract(r"^(P[12])_(MP)_(\d+)$") -long["site_id"] = m[0] -long["depth_cm"] = m[2] -long = long[long["MP"].notna()].copy() -long["depth_m"] = pd.to_numeric(long["depth_cm"], errors="coerce") / 100 -long["replicate"] = 1 -long["is_timeseries"] = True -long["volumetric_water_content_m3_m3"] = np.nan -long["gravimetric_water_content_gH2O_gs"] = np.nan -long["water_potential_kPa"] = pd.to_numeric(long["MP"], errors="coerce") -long = long[long["water_potential_kPa"].notna()] - -df24_harmonized = ensure_harmonized_cols(long) -harmonized_data.append(df24_harmonized) -harmonized_ids.append(dsid(idx)) - -loc24 = mdf24.rename(columns={"ID": "site_id", "Latitude": "latitude", "Longitude": "longitude"})[ - ["site_id", "latitude", "longitude"] -].copy() -loc24 = loc24[loc24["site_id"].isin(df24_harmonized["site_id"].dropna().unique())].drop_duplicates() -loc24["source_dataset_id"] = dsid(idx) -loc24 = add_loc_qc(loc24) -loc_data.append(loc24) - -# %% -# ============================================================= -# Dataset 25 -# ============================================================= -idx = 25 -ddf25_conifer = read_ds_csv(idx, "Carbone_conifer.csv") -ddf25_aspen = read_ds_csv(idx, "Carbone_aspen.csv", header=None) -ddf25_aspen.columns = ddf25_conifer.columns - -ddf25 = pd.concat( - [ - ddf25_aspen.assign(site_id="aspen"), - ddf25_conifer.assign(site_id="conifer"), - ], - ignore_index=True, -) - -x = ddf25.copy() -x["datetime_UTC"] = parse_local_to_utc( - x["Year"].astype(str) + "-" + x["Month"].astype(str) + "-" + x["Day"].astype(str) + "-" + x["Hour"].astype(str), - "%Y-%m-%d-%H", - "America/Denver", -) -x = x.sort_values(["site_id", "datetime_UTC"]) -x["interval_min"] = x.groupby("site_id")["datetime_UTC"].diff().dt.total_seconds() / 60.0 - -vwc_cols = [c for c in x.columns if re.search(r"vwc", c, flags=re.IGNORECASE)] -long = x.melt( - id_vars=["datetime_UTC", "site_id", "interval_min"], - value_vars=vwc_cols, - var_name="name", - value_name="volumetric_water_content_m3_m3", -) -long["depth_m"] = pd.to_numeric(long["name"].str.extract(r"(\d+)")[0], errors="coerce") / 100 -long["replicate"] = 1 -long["is_timeseries"] = True -long["water_potential_kPa"] = np.nan -long["gravimetric_water_content_gH2O_gs"] = np.nan -long["volumetric_water_content_m3_m3"] = pd.to_numeric(long["volumetric_water_content_m3_m3"], errors="coerce") -long = long[long["volumetric_water_content_m3_m3"].notna()] - -df25_harmonized = ensure_harmonized_cols(long) -harmonized_data.append(df25_harmonized) -harmonized_ids.append(dsid(idx)) - -loc25 = pd.DataFrame( - { - "site_id": ["aspen", "conifer"], - "latitude": [38.9581589, 38.9581589], - "longitude": [-106.9855848, -106.9855848], - } -).drop_duplicates() -loc25["source_dataset_id"] = dsid(idx) -loc25 = add_loc_qc(loc25) -loc_data.append(loc25) - -# %% -# ============================================================= -# Dataset 26 -# ============================================================= -idx = 26 -f26 = as_list(map_json[idx]["data_payload_files"])[0] -m26 = as_list(map_json[REF_IDX]["location_metadata_files"])[0] - -ddf26 = read_ds_csv(idx, f26, index_col=0) -mdf26 = read_ds_csv(REF_IDX, m26) - -x = ddf26.copy() -x["datetime_UTC"] = parse_local_to_utc(x["Collection date"], "%m/%d/%y", "America/Denver") -x["site_id"] = x["SampleSiteCode"] -x["depth_m"] = (pd.to_numeric(x["Top sample depth_cm"], errors="coerce") + pd.to_numeric(x["Bottom sample depth_cm"], errors="coerce")) / 2 / 100 -x["replicate"] = 1 -x["is_timeseries"] = False -x["interval_min"] = np.nan -x["volumetric_water_content_m3_m3"] = pd.to_numeric(x["water content %vol"], errors="coerce") / 100 -x["gravimetric_water_content_gH2O_gs"] = np.nan -x["water_potential_kPa"] = np.nan - -df26_harmonized = ensure_harmonized_cols(x) -harmonized_data.append(df26_harmonized) -harmonized_ids.append(dsid(idx)) - -loc26 = mdf26.rename(columns={"Location_ID": "site_id", "Latitude": "latitude", "Longitude": "longitude"})[ - ["site_id", "latitude", "longitude"] -].copy() -loc26 = loc26[loc26["site_id"].isin(df26_harmonized["site_id"].dropna().unique())] -loc26["source_dataset_id"] = dsid(idx) -loc26 = add_loc_qc(loc26) -loc_data.append(loc26) - -# %% -# ============================================================= -# Dataset 27 -# ============================================================= -idx = 27 -f27 = as_list(map_json[idx]["data_payload_files"]) -m27 = as_list(map_json[REF_IDX]["location_metadata_files"])[0] - -ls27 = [read_ds_csv(idx, x) for x in f27] -mdf27 = read_ds_csv(REF_IDX, m27) - -ls27_h = [] -for i, d in enumerate(ls27): - x = d.copy() - dt = pd.to_datetime(x["Time"], format="%m/%d/%Y %H:%M", errors="coerce") - dt = dt.dt.tz_localize("Etc/GMT+7", ambiguous="NaT", nonexistent="shift_forward") - x["datetime_UTC"] = dt.dt.tz_convert("UTC") - x["interval_min"] = interval_min(x["datetime_UTC"]) - - cols = [c for c in x.columns if re.search(r"^S[1-4]_wc_\(m3/m3\)$|^S5_wp_\(kPa\)$", c)] - long = x.melt(id_vars=["datetime_UTC", "interval_min"], value_vars=cols, var_name="name", value_name="value") - - long["variable"] = np.select( - [ - long["name"].str.contains("_wc_", regex=False, na=False), - long["name"].str.contains("_wp_", regex=False, na=False), - ], - ["volumetric_water_content_m3_m3", "water_potential_kPa"], - default=np.nan, - ) - - fname = f27[i] - long["depth_m"] = np.select( - [ - long["name"].str.contains(r"^S1_", regex=True, na=False) & pd.Series([bool(re.search(r"ER-PHS4", fname))] * len(long)), - long["name"].str.contains(r"^S1_", regex=True, na=False), - long["name"].str.contains(r"^S2_", regex=True, na=False), - long["name"].str.contains(r"^S3_", regex=True, na=False), - long["name"].str.contains(r"^S4_", regex=True, na=False), - long["name"].str.contains(r"^S5_", regex=True, na=False), - ], - [1.06, 1.15, 0.60, 0.30, 0.10, 0.30], - default=np.nan, - ) - - site_match = re.search(r"ER-PHS[1-4]", fname) - site_id_val = site_match.group(0) if site_match else np.nan - long["site_id"] = site_id_val - long["replicate"] = 1 - long["is_timeseries"] = True - - long = long[long["value"].notna() & long["variable"].notna() & long["depth_m"].notna()].copy() - wide = ( - long.pivot_table( - index=["datetime_UTC", "site_id", "depth_m", "replicate", "is_timeseries", "interval_min"], - columns="variable", - values="value", - aggfunc="first", - ) - .reset_index() - ) - - wide["volumetric_water_content_m3_m3"] = pd.to_numeric(wide.get("volumetric_water_content_m3_m3"), errors="coerce") - wide["water_potential_kPa"] = pd.to_numeric(wide.get("water_potential_kPa"), errors="coerce") - wide["gravimetric_water_content_gH2O_gs"] = np.nan - - ls27_h.append(ensure_harmonized_cols(wide)) - -df27_harmonized = pd.concat(ls27_h, ignore_index=True) -harmonized_data.append(df27_harmonized) -harmonized_ids.append(dsid(idx)) - -loc27 = mdf27.rename(columns={"Location_ID": "site_id", "Latitude": "latitude", "Longitude": "longitude"})[ - ["site_id", "latitude", "longitude"] -].copy() -loc27 = loc27[loc27["site_id"].isin(df27_harmonized["site_id"].dropna().unique())] -loc27["source_dataset_id"] = dsid(idx) -loc27 = add_loc_qc(loc27) -loc_data.append(loc27) - -# %% -# ============================================================= -# Location deduplication and harmonization -# ============================================================= - -# Here we find sites with identical names across datasets or -# sites within a small distance thresholds and collapse them to a -# single uuid - -import uuid -from itertools import combinations - -# Concat location data -loc_df = pd.concat(loc_data, ignore_index=True) -loc_df['latitude'] = pd.to_numeric(loc_df['latitude'], errors='coerce') -loc_df['longitude'] = pd.to_numeric(loc_df['longitude'], errors='coerce') - -# Set thresholds and toggle site_id matching -# Matching thresholds (tune these) -coord_match_meters_strict = 5 # highly likely same footprint - -# If TRUE, identical site_id across datasets is considered evidence -use_cross_dataset_site_id = True - - -# Function to normalize names for comparison -def normalize_name(x) -> str: - if pd.isna(x): - return "" - x = str(x).lower() - x = re.sub(r"[^a-z0-9]+", "", x) - return x.strip() - - -# Function to find safe distance between sites using Haversine formula -def safe_dist_m(lon1: float, lat1: float, lon2: float, lat2: float) -> float: - """Calculate distance in meters between two lat/lon points using Haversine formula""" - if any(pd.isna([lon1, lat1, lon2, lat2])): - return np.inf - - # Haversine formula - R = 6371000 # Earth radius in meters - phi1 = np.radians(lat1) - phi2 = np.radians(lat2) - delta_phi = np.radians(lat2 - lat1) - delta_lambda = np.radians(lon2 - lon1) - - a = np.sin(delta_phi / 2) ** 2 + np.cos(phi1) * np.cos(phi2) * np.sin(delta_lambda / 2) ** 2 - c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a)) - - return R * c - - -# Add site_name column (equivalent to site_id) -loc_df['site_name'] = loc_df['site_id'] - -# Prepare location dataframe -loc_df = loc_df.reset_index(drop=True) -loc_df['row_id'] = loc_df.index -loc_df['site_id'] = loc_df['site_id'].astype(str) -loc_df['source_dataset_id'] = loc_df['source_dataset_id'].astype(str) -loc_df['site_name_norm'] = loc_df['site_name'].apply(normalize_name) -loc_df['has_coords'] = loc_df['latitude'].notna() & loc_df['longitude'].notna() - -n = len(loc_df) -if n == 0: - raise ValueError("No rows in location file.") - - -# Build pairwise links -# Link rows likely to be same location -def is_match_pair(i: int, j: int, df: pd.DataFrame) -> bool: - """Check if two rows should be considered the same location""" - a = df.iloc[i] - b = df.iloc[j] - - # Strong coordinate match - d_m = safe_dist_m(a['longitude'], a['latitude'], b['longitude'], b['latitude']) - if np.isfinite(d_m) and d_m <= coord_match_meters_strict: - return True - - # Same site_id across different datasets - if use_cross_dataset_site_id: - same_site_id = ( - pd.notna(a['site_id']) and pd.notna(b['site_id']) and - str(a['site_id']) != "" and str(b['site_id']) != "" and - a['site_id'] == b['site_id'] - ) - if same_site_id: - return True - - return False - - -# Find all matching pairs using Union-Find algorithm -class UnionFind: - def __init__(self, n): - self.parent = list(range(n)) - self.rank = [0] * n - - def find(self, x): - if self.parent[x] != x: - self.parent[x] = self.find(self.parent[x]) - return self.parent[x] - - def union(self, x, y): - px, py = self.find(x), self.find(y) - if px == py: - return - if self.rank[px] < self.rank[py]: - px, py = py, px - self.parent[py] = px - if self.rank[px] == self.rank[py]: - self.rank[px] += 1 - - -uf = UnionFind(n) - -# Check all pairs and union matching locations -for i, j in combinations(range(n), 2): - if is_match_pair(i, j, loc_df): - uf.union(i, j) - -# Assign component IDs -comp_map = {} -next_comp_id = 0 -for i in range(n): - root = uf.find(i) - if root not in comp_map: - comp_map[root] = next_comp_id - next_comp_id += 1 - -loc_df['location_component_id'] = [comp_map[uf.find(i)] for i in range(n)] - -# Assign stable UUID per component -comp_ids = sorted(loc_df['location_component_id'].unique()) -uuid_map = pd.DataFrame({ - 'location_component_id': comp_ids, - 'harmonized_location_uuid': [str(uuid.uuid4()) for _ in comp_ids] -}) - -loc_df = loc_df.merge(uuid_map, on='location_component_id', how='left') - -# Optional canonical fields per UUID (centroid + representative name) -canon = loc_df.groupby('harmonized_location_uuid').agg( - latitude_harmonized=('latitude', lambda x: np.nan if x.isna().all() else x.mean()), - longitude_harmonized=('longitude', lambda x: np.nan if x.isna().all() else x.mean()), - n_records_in_uuid=('row_id', 'count'), - n_datasets_in_uuid=('source_dataset_id', 'nunique') -).reset_index() - -# Join all together -loc_df = loc_df.merge(canon, on='harmonized_location_uuid', how='left') -loc_df = loc_df.drop(columns=['site_name', 'site_name_norm', 'has_coords', 'row_id', 'location_component_id']) -loc_df = loc_df.sort_values(['source_dataset_id', 'site_id', 'harmonized_location_uuid']).reset_index(drop=True) - -# Quick QA summary -qa = loc_df.groupby('harmonized_location_uuid').size().reset_index(name='n') -qa = qa.sort_values('n', ascending=False) -qa['flag_multi'] = qa['n'] > 1 - -print(f"UUID groups with >1 member: {qa['flag_multi'].sum()} / {len(qa)}") - -# %% -# ============================================================= -# Write -# ============================================================= - -for ds_identifier, df in zip(harmonized_ids, harmonized_data): - out_file = OUT_DIR / f"{ds_identifier}_harmonized.csv" - df.to_csv(out_file, index=False) - -loc_df = pd.concat(loc_data, ignore_index=True) -loc_df.to_csv(OUT_DIR / "location_data_harmonized.csv", index=False) \ No newline at end of file diff --git a/data/gold/expert_code/harmonize_sm/README.md b/data/gold/expert_code/harmonize_sm/README.md new file mode 100644 index 0000000..ce244ba --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/README.md @@ -0,0 +1,51 @@ +# Expert soil-moisture harmonizer (modular) + +This is the expert harmonization "gold" code, refactored from a single +1200-line monolith into a small package: one shared library plus one module +per dataset. + +``` +harmonize_sm/ +├── common.py # Context (mapping JSON + dsid/ds_path/read_ds_csv), +│ # pure helpers, location dedup + CSV writer +├── dataset_01.py # one harmonize(ctx) -> DatasetResult per dataset +├── dataset_02.py +├── ... # indices 1-10, 15-18, 23-27 (the 19 expert datasets) +├── dataset_27.py +├── datasets.py # registry: {index -> harmonize fn} +└── run.py # driver with --holdout +``` + +## Why it's split this way + +Each dataset block in the old monolith only appended to shared accumulators and +the sole cross-block dependency was the reference dataset (index 0, never a +hold-out target). So each block is independent and becomes a self-contained +`harmonize(ctx)` that *returns* its harmonized frame + location frame(s) instead +of mutating globals. The bodies are verbatim copies of the monolith — only the +wiring changed — so the harmonized outputs are identical. + +## Running + +```bash +# Full run (writes per-dataset CSVs + location_data_harmonized.csv) +python run.py + +# Leave-one-cluster-out: just don't run the held-out datasets +python run.py --holdout 1,2,3,6,16,27 # cluster_1 "sg_ph_micromet" +``` + +## Grouped leave-one-cluster-out + +This split is what makes grouped LOO trivial and is why the old cell-splicer +(`ablate_monolith`) is no longer needed: + +- **Executable** (regenerate exemplar outputs without the held-out datasets): + `run.py --holdout ...` — the held-out modules simply aren't run. +- **Code context** (show the agent held-out-free reference code): feed + `common.py` plus the kept `dataset_NN.py` files. + +The harness wraps both via `src/folds/expert_harmonizer.py` +(`run`, `kept_module_paths`, `assemble_source`). Requesting a hold-out that has +no module (the reference dataset 0, or the excluded datasets 11-14 / 19-22) is +rejected, so there is no coupling to `config/cv_folds.yaml`. diff --git a/data/gold/expert_code/harmonize_sm/common.py b/data/gold/expert_code/harmonize_sm/common.py new file mode 100644 index 0000000..92194c6 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/common.py @@ -0,0 +1,337 @@ +"""Shared library for the expert soil-moisture harmonization. + +This is the ``common`` half of the refactored expert script. It holds +everything the per-dataset modules (``dataset_NN.py``) share: + +* :class:`Context` — the loaded mapping JSON plus the dataset-access helpers + (``dsid`` / ``ds_path`` / ``read_ds_csv``) that used to be module-level + globals in the monolith. +* the pure transformation helpers (``parse_local_to_utc`` etc.), unchanged. +* :class:`DatasetResult` — what each ``harmonize(ctx)`` returns. +* :func:`harmonize_locations` — the cross-dataset location dedup footer. +* :func:`write_outputs` — the per-dataset + location CSV writer. + +The bodies of the helpers and the two footers are verbatim copies of the +expert monolith, so the harmonized outputs are identical; only the wiring +(globals -> Context, accumulator lists -> return values) changed. +""" +from __future__ import annotations + +import json +from dataclasses import dataclass, field +from pathlib import Path + +import numpy as np +import pandas as pd +from pyproj import Transformer + +# ============================================================= +# Config / Paths (defaults; overridable via Context.load) +# ============================================================= + +DEFAULT_BASE_DIR = Path(Path.home(), "ess-dive_wfsfa_soil_datasets") +DEFAULT_OUT_DIR = Path("data", "processed", "ess-dive_wfsfa_soil_datasets") + +# map_json[0] is the reference dataset (lookup) +REF_IDX = 0 + + +@dataclass +class DatasetResult: + """What a per-dataset ``harmonize(ctx)`` returns. + + ``locations`` is a list because some datasets (e.g. index 23) contribute + more than one location frame. + """ + + dataset_id: str + harmonized: pd.DataFrame + locations: list = field(default_factory=list) + + +@dataclass +class Context: + """Shared, read-only state passed to every per-dataset harmonizer.""" + + map_json: list[dict] + base_dir: Path = DEFAULT_BASE_DIR + out_dir: Path = DEFAULT_OUT_DIR + ref_idx: int = REF_IDX + + @classmethod + def load( + cls, + base_dir: Path | None = None, + out_dir: Path | None = None, + mapping_path: Path | None = None, + ) -> "Context": + """Load the gold mapping JSON and build a Context. + + Mirrors the monolith's Config cell: the mapping lives next to the + processed outputs, and the output directory is created on load. + """ + base_dir = Path(base_dir) if base_dir is not None else DEFAULT_BASE_DIR + out_dir = Path(out_dir) if out_dir is not None else DEFAULT_OUT_DIR + if mapping_path is None: + mapping_path = out_dir / "sm_data_harmonization_mapping.json" + out_dir.mkdir(parents=True, exist_ok=True) + with Path(mapping_path).open("r", encoding="utf-8") as f: + map_json: list[dict] = json.load(f) + return cls(map_json=map_json, base_dir=base_dir, out_dir=out_dir, ref_idx=REF_IDX) + + # --- dataset-access helpers (were module-level globals in the monolith) --- + + def dsid(self, idx: int) -> str: + return self.map_json[idx]["dataset_identifier"] + + def ds_path(self, idx: int) -> Path: + return self.base_dir / self.dsid(idx) + + def read_ds_csv( + self, idx: int, filename: str, encoding="utf-8", errors="ignore", **kwargs + ) -> pd.DataFrame: + return pd.read_csv(self.ds_path(idx) / filename, encoding=encoding, **kwargs) + + +# ============================================================= +# Helpers (verbatim from the monolith) +# ============================================================= + + +def as_list(x): + if x is None: + return [] + return x if isinstance(x, list) else [x] + + +def parse_local_to_utc(series: pd.Series, fmt: str | None, tz: str) -> pd.Series: + dt = pd.to_datetime(series, format=fmt, errors="coerce") + dt = dt.dt.tz_localize(tz, ambiguous="NaT", nonexistent="shift_forward") + return dt.dt.tz_convert("UTC") + + +def interval_min(s: pd.Series) -> pd.Series: + return s.diff().dt.total_seconds() / 60.0 + + +def utm32613_to_latlon(df: pd.DataFrame, e_col: str, n_col: str) -> pd.DataFrame: + tr = Transformer.from_crs("EPSG:32613", "EPSG:4326", always_xy=True) + e = pd.to_numeric(df[e_col], errors="coerce").values + n = pd.to_numeric(df[n_col], errors="coerce").values + lon, lat = tr.transform(e, n) + out = df.copy() + out["longitude"] = lon + out["latitude"] = lat + return out + + +def ensure_harmonized_cols(df: pd.DataFrame) -> pd.DataFrame: + cols = [ + "datetime_UTC", + "site_id", + "depth_m", + "replicate", + "is_timeseries", + "interval_min", + "volumetric_water_content_m3_m3", + "gravimetric_water_content_gH2O_gs", + "water_potential_kPa", + ] + for c in cols: + if c not in df.columns: + df[c] = np.nan + return df[cols] + + +def add_loc_qc(df: pd.DataFrame) -> pd.DataFrame: + if "qc_flag" not in df.columns: + df["qc_flag"] = np.where(df["latitude"].isna() | df["longitude"].isna(), "g2", np.nan) + return df + + +# ============================================================= +# Location deduplication and harmonization (footer; verbatim body) +# ============================================================= + + +def harmonize_locations(loc_data: list[pd.DataFrame]) -> pd.DataFrame: + """Collapse sites with identical names / near-identical coordinates to one uuid. + + Verbatim port of the monolith's dedup cell. Returns the deduplicated + ``loc_df`` and prints the QA summary. Note: as in the expert script, the + written ``location_data_harmonized.csv`` is the *raw* concat of ``loc_data`` + (see :func:`write_outputs`), not this deduplicated frame. + """ + import re + import uuid + from itertools import combinations + + # Concat location data + loc_df = pd.concat(loc_data, ignore_index=True) + loc_df['latitude'] = pd.to_numeric(loc_df['latitude'], errors='coerce') + loc_df['longitude'] = pd.to_numeric(loc_df['longitude'], errors='coerce') + + # Set thresholds and toggle site_id matching + # Matching thresholds (tune these) + coord_match_meters_strict = 5 # highly likely same footprint + + # If TRUE, identical site_id across datasets is considered evidence + use_cross_dataset_site_id = True + + # Function to normalize names for comparison + def normalize_name(x) -> str: + if pd.isna(x): + return "" + x = str(x).lower() + x = re.sub(r"[^a-z0-9]+", "", x) + return x.strip() + + # Function to find safe distance between sites using Haversine formula + def safe_dist_m(lon1: float, lat1: float, lon2: float, lat2: float) -> float: + """Calculate distance in meters between two lat/lon points using Haversine formula""" + if any(pd.isna([lon1, lat1, lon2, lat2])): + return np.inf + + # Haversine formula + R = 6371000 # Earth radius in meters + phi1 = np.radians(lat1) + phi2 = np.radians(lat2) + delta_phi = np.radians(lat2 - lat1) + delta_lambda = np.radians(lon2 - lon1) + + a = np.sin(delta_phi / 2) ** 2 + np.cos(phi1) * np.cos(phi2) * np.sin(delta_lambda / 2) ** 2 + c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a)) + + return R * c + + # Add site_name column (equivalent to site_id) + loc_df['site_name'] = loc_df['site_id'] + + # Prepare location dataframe + loc_df = loc_df.reset_index(drop=True) + loc_df['row_id'] = loc_df.index + loc_df['site_id'] = loc_df['site_id'].astype(str) + loc_df['source_dataset_id'] = loc_df['source_dataset_id'].astype(str) + loc_df['site_name_norm'] = loc_df['site_name'].apply(normalize_name) + loc_df['has_coords'] = loc_df['latitude'].notna() & loc_df['longitude'].notna() + + n = len(loc_df) + if n == 0: + raise ValueError("No rows in location file.") + + # Build pairwise links + # Link rows likely to be same location + def is_match_pair(i: int, j: int, df: pd.DataFrame) -> bool: + """Check if two rows should be considered the same location""" + a = df.iloc[i] + b = df.iloc[j] + + # Strong coordinate match + d_m = safe_dist_m(a['longitude'], a['latitude'], b['longitude'], b['latitude']) + if np.isfinite(d_m) and d_m <= coord_match_meters_strict: + return True + + # Same site_id across different datasets + if use_cross_dataset_site_id: + same_site_id = ( + pd.notna(a['site_id']) and pd.notna(b['site_id']) and + str(a['site_id']) != "" and str(b['site_id']) != "" and + a['site_id'] == b['site_id'] + ) + if same_site_id: + return True + + return False + + # Find all matching pairs using Union-Find algorithm + class UnionFind: + def __init__(self, n): + self.parent = list(range(n)) + self.rank = [0] * n + + def find(self, x): + if self.parent[x] != x: + self.parent[x] = self.find(self.parent[x]) + return self.parent[x] + + def union(self, x, y): + px, py = self.find(x), self.find(y) + if px == py: + return + if self.rank[px] < self.rank[py]: + px, py = py, px + self.parent[py] = px + if self.rank[px] == self.rank[py]: + self.rank[px] += 1 + + uf = UnionFind(n) + + # Check all pairs and union matching locations + for i, j in combinations(range(n), 2): + if is_match_pair(i, j, loc_df): + uf.union(i, j) + + # Assign component IDs + comp_map = {} + next_comp_id = 0 + for i in range(n): + root = uf.find(i) + if root not in comp_map: + comp_map[root] = next_comp_id + next_comp_id += 1 + + loc_df['location_component_id'] = [comp_map[uf.find(i)] for i in range(n)] + + # Assign stable UUID per component + comp_ids = sorted(loc_df['location_component_id'].unique()) + uuid_map = pd.DataFrame({ + 'location_component_id': comp_ids, + 'harmonized_location_uuid': [str(uuid.uuid4()) for _ in comp_ids] + }) + + loc_df = loc_df.merge(uuid_map, on='location_component_id', how='left') + + # Optional canonical fields per UUID (centroid + representative name) + canon = loc_df.groupby('harmonized_location_uuid').agg( + latitude_harmonized=('latitude', lambda x: np.nan if x.isna().all() else x.mean()), + longitude_harmonized=('longitude', lambda x: np.nan if x.isna().all() else x.mean()), + n_records_in_uuid=('row_id', 'count'), + n_datasets_in_uuid=('source_dataset_id', 'nunique') + ).reset_index() + + # Join all together + loc_df = loc_df.merge(canon, on='harmonized_location_uuid', how='left') + loc_df = loc_df.drop(columns=['site_name', 'site_name_norm', 'has_coords', 'row_id', 'location_component_id']) + loc_df = loc_df.sort_values(['source_dataset_id', 'site_id', 'harmonized_location_uuid']).reset_index(drop=True) + + # Quick QA summary + qa = loc_df.groupby('harmonized_location_uuid').size().reset_index(name='n') + qa = qa.sort_values('n', ascending=False) + qa['flag_multi'] = qa['n'] > 1 + + print(f"UUID groups with >1 member: {qa['flag_multi'].sum()} / {len(qa)}") + + return loc_df + + +# ============================================================= +# Write (footer; verbatim body) +# ============================================================= + + +def write_outputs( + out_dir: Path, + harmonized_ids: list[str], + harmonized_data: list[pd.DataFrame], + loc_data: list[pd.DataFrame], +) -> None: + """Write per-dataset harmonized CSVs and the concatenated location CSV.""" + out_dir = Path(out_dir) + out_dir.mkdir(parents=True, exist_ok=True) + + for ds_identifier, df in zip(harmonized_ids, harmonized_data): + out_file = out_dir / f"{ds_identifier}_harmonized.csv" + df.to_csv(out_file, index=False) + + loc_df = pd.concat(loc_data, ignore_index=True) + loc_df.to_csv(out_dir / "location_data_harmonized.csv", index=False) diff --git a/data/gold/expert_code/harmonize_sm/dataset_01.py b/data/gold/expert_code/harmonize_sm/dataset_01.py new file mode 100644 index 0000000..59500d2 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_01.py @@ -0,0 +1,106 @@ +"""Expert harmonization for dataset index 1. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import re +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + parse_local_to_utc, + interval_min, + utm32613_to_latlon, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 1, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 1 + # ============================================================= + idx = 1 + f1 = as_list(map_json[idx]["data_payload_files"]) + m1 = as_list(map_json[idx]["location_metadata_files"])[0] + + ls1 = [read_ds_csv(idx, x) for x in f1] + mdf1 = read_ds_csv(idx, m1) + + ls1_h = [] + for i, d in enumerate(ls1): + cols = [c for c in d.columns if re.search(r"Water_Content|VWC|Potential", c)] + x = d[["Time"] + cols].copy() + + x["datetime_UTC"] = parse_local_to_utc(x["Time"], "%Y-%m-%d %H:%M:%S", "America/Denver") + x["interval_min"] = interval_min(x["datetime_UTC"]) + x = x.drop(columns=["Time"]) + + long = x.melt(id_vars=["datetime_UTC", "interval_min"], var_name="name", value_name="value") + m = long["name"].str.extract(r"^(.*)_at_(.*)$") + long["var_name"] = m[0] + long["depth"] = m[1] + long = long.dropna(subset=["depth"]) + + long["value"] = pd.to_numeric(long["value"], errors="coerce") + long.loc[long["value"] == -9999, "value"] = np.nan + long = long.dropna(subset=["value"]) + + long["dest_var"] = np.where( + long["var_name"].str.contains(r"Water_Content|VWC", regex=True, na=False), + "volumetric_water_content_m3_m3", + np.where( + long["var_name"].str.contains(r"Potential", regex=True, na=False), + "water_potential_kPa", + np.nan, + ), + ) + long = long.dropna(subset=["dest_var"]) + + long["replicate"] = long.groupby(["datetime_UTC", "depth", "dest_var"]).cumcount() + 1 + + wide = ( + long.pivot_table( + index=["datetime_UTC", "interval_min", "depth", "replicate"], + columns="dest_var", + values="value", + aggfunc="first", + ) + .reset_index() + ) + + wide["depth_m"] = pd.to_numeric(wide["depth"].str.replace("cm", "", regex=False), errors="coerce") / 100 + wide["site_id"] = re.sub(r"\.csv$", "", f1[i]) + wide["is_timeseries"] = True + wide["gravimetric_water_content_gH2O_gs"] = np.nan + + ls1_h.append(ensure_harmonized_cols(wide)) + + df1_harmonized = pd.concat(ls1_h, ignore_index=True) + __harmonized = df1_harmonized + __dataset_id = dsid(idx) + + loc1_tmp = utm32613_to_latlon(mdf1, "Easting", "Northing") + loc1 = pd.DataFrame( + { + "site_id": loc1_tmp["Name"], + "latitude": loc1_tmp["latitude"], + "longitude": loc1_tmp["longitude"], + } + ) + loc1["source_dataset_id"] = dsid(idx) + loc1 = add_loc_qc(loc1) + __locations.append(loc1) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_02.py b/data/gold/expert_code/harmonize_sm/dataset_02.py new file mode 100644 index 0000000..2a5c871 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_02.py @@ -0,0 +1,75 @@ +"""Expert harmonization for dataset index 2. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import re +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + parse_local_to_utc, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 2, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 2 + # ============================================================= + idx = 2 + f2 = as_list(map_json[idx]["data_payload_files"]) + m2 = as_list(map_json[idx]["location_metadata_files"])[0] + + ls2 = [read_ds_csv(idx, x) for x in f2] + mdf2 = read_ds_csv(idx, m2) + + ls2_h = [] + for i, d in enumerate(ls2): + cols = [c for c in d.columns if re.search(r"Moisture", c)] + x = d[["DateTime"] + cols].copy() + x["datetime_UTC"] = parse_local_to_utc(x["DateTime"], "%m/%d/%Y %I:%M:%S %p", "America/Denver") + x["interval_min"] = interval_min(x["datetime_UTC"]) + x = x.drop(columns=["DateTime"]) + + long = x.melt( + id_vars=["datetime_UTC", "interval_min"], + var_name="name", + value_name="volumetric_water_content_m3_m3", + ) + long["depth_m"] = ( + pd.to_numeric(long["name"].str.split("_").str[-1].str.replace("cm", "", regex=False), errors="coerce") / 100 + ) + long["site_id"] = re.sub(r"\.csv$", "", f2[i]) + long["is_timeseries"] = True + long["water_potential_kPa"] = np.nan + long["gravimetric_water_content_gH2O_gs"] = np.nan + long["replicate"] = long.groupby(["datetime_UTC", "depth_m", "name"]).cumcount() + 1 + + ls2_h.append(ensure_harmonized_cols(long)) + + df2_harmonized = pd.concat(ls2_h, ignore_index=True) + __harmonized = df2_harmonized + __dataset_id = dsid(idx) + + loc2 = mdf2.rename(columns={"Name": "site_id", "Lat ": "latitude", "Lon": "longitude"})[ + ["site_id", "latitude", "longitude"] + ].copy() + loc2["source_dataset_id"] = dsid(idx) + loc2 = add_loc_qc(loc2) + __locations.append(loc2) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_03.py b/data/gold/expert_code/harmonize_sm/dataset_03.py new file mode 100644 index 0000000..725d49f --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_03.py @@ -0,0 +1,70 @@ +"""Expert harmonization for dataset index 3. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + parse_local_to_utc, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 3, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 3 + # ============================================================= + idx = 3 + f3 = as_list(map_json[idx]["data_payload_files"])[0] + m3 = as_list(map_json[idx]["location_metadata_files"])[0] + + df3 = read_ds_csv(idx, f3) + mdf3 = read_ds_csv(idx, m3) + + x = df3.copy() + x["datetime_UTC"] = parse_local_to_utc(x["TIMESTAMP"], "%Y-%m-%d %H:%M", "America/Denver") + x["interval_min"] = interval_min(x["datetime_UTC"]) + + mp_cols = [c for c in x.columns if "_MP" in c] + for c in mp_cols: + x[c] = pd.to_numeric(x[c], errors="coerce") + + long = x.melt(id_vars=["datetime_UTC", "interval_min"], value_vars=mp_cols, var_name="name", value_name="water_potential_kPa") + long["depth_m"] = pd.to_numeric( + long["name"].str.replace(r"[._]", " ", regex=True).str.split().str[1].str.replace("cm", "", regex=False), + errors="coerce", + ) / 100 + long["site_id"] = "ER-LLN1" + long["is_timeseries"] = True + long["water_potential_kPa"] = long["water_potential_kPa"].where(~np.isnan(long["water_potential_kPa"]), np.nan) + long["volumetric_water_content_m3_m3"] = np.nan + long["gravimetric_water_content_gH2O_gs"] = np.nan + long["replicate"] = long.groupby(["datetime_UTC", "depth_m"]).cumcount() + 1 + + df3_harmonized = ensure_harmonized_cols(long) + __harmonized = df3_harmonized + __dataset_id = dsid(idx) + + loc3 = mdf3.rename(columns={"Location_ID": "site_id", "Latitude": "latitude", "Longitude": "longitude"})[ + ["site_id", "latitude", "longitude"] + ].copy() + loc3["source_dataset_id"] = dsid(idx) + loc3 = add_loc_qc(loc3) + __locations.append(loc3) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_04.py b/data/gold/expert_code/harmonize_sm/dataset_04.py new file mode 100644 index 0000000..8e9bf54 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_04.py @@ -0,0 +1,66 @@ +"""Expert harmonization for dataset index 4. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 4, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 4 + # ============================================================= + idx = 4 + f4 = as_list(map_json[idx]["data_payload_files"])[0] + m4 = as_list(map_json[idx]["location_metadata_files"])[0] + ddf4 = read_ds_csv(idx, f4) + mdf4 = read_ds_csv(idx, m4) + + x = pd.concat([ddf4.reset_index(drop=True), mdf4.reset_index(drop=True)], axis=1) + dt = pd.to_datetime(x["datetime"], errors="coerce").dt.tz_localize("America/Denver", ambiguous="NaT", nonexistent="shift_forward") + x["datetime_UTC"] = dt.dt.tz_convert("UTC") + x["site_id"] = x["site"] + x["depth_m"] = np.nan + x["replicate"] = 1 + x["is_timeseries"] = True + x["interval_min"] = 24 * 60 + x["volumetric_water_content_m3_m3"] = pd.to_numeric(x["swc"], errors="coerce") + x["water_potential_kPa"] = pd.to_numeric(x["swp"], errors="coerce") + x["gravimetric_water_content_gH2O_gs"] = np.nan + x = x.sort_values(["datetime_UTC", "site_id"]) + x = x[x["site_id"] != "tb"] + + df4_harmonized = ensure_harmonized_cols(x) + __harmonized = df4_harmonized + __dataset_id = dsid(idx) + + loc4 = pd.DataFrame( + { + "site_id": ["ph1", "ph2", "sg5"], + "latitude": [38.92, 38.922583, 38.926250], + "longitude": [-106.95, -106.947288, -106.98], + } + ) + loc4["source_dataset_id"] = dsid(idx) + loc4 = add_loc_qc(loc4) + __locations.append(loc4) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_05.py b/data/gold/expert_code/harmonize_sm/dataset_05.py new file mode 100644 index 0000000..6007d47 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_05.py @@ -0,0 +1,78 @@ +"""Expert harmonization for dataset index 5. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + interval_min, + ensure_harmonized_cols, +) + + +def harmonize(ctx): + """Harmonize dataset 5, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 5 + # ============================================================= + idx = 5 + f5 = as_list(map_json[idx]["data_payload_files"])[0] + m5 = as_list(map_json[REF_IDX]["location_metadata_files"])[0] + + ddf5 = read_ds_csv(idx, f5) + mdf5 = read_ds_csv(REF_IDX, m5) + + x = ddf5.copy() + dt_local = pd.to_datetime(x["Date Time"], errors="coerce").dt.tz_localize("America/Denver", ambiguous="NaT", nonexistent="shift_forward") + x["datetime_UTC"] = dt_local.dt.tz_convert("UTC") + x["site_id"] = x["SFA_Location"] + x["depth_m"] = np.nan + x["replicate"] = 1 + x["is_timeseries"] = True + x["interval_min"] = 60 + x["volumetric_water_content_m3_m3"] = pd.to_numeric(x["Measurement"], errors="coerce") + x["water_potential_kPa"] = np.nan + x["gravimetric_water_content_gH2O_gs"] = np.nan + + df5_harmonized = ensure_harmonized_cols(x) + __harmonized = df5_harmonized + __dataset_id = dsid(idx) + + loc5 = pd.DataFrame({"site_id": df5_harmonized["site_id"].dropna().unique()}) + + loc5a = loc5.merge( + mdf5[["Location_ID", "Latitude", "Longitude"]], + left_on="site_id", + right_on="Location_ID", + how="left", + ).drop(columns=["Location_ID"]) + + loc5b = loc5.merge( + mdf5[["Parent_Location_ID", "Latitude", "Longitude"]], + left_on="site_id", + right_on="Parent_Location_ID", + how="left", + ).drop(columns=["Parent_Location_ID"]) + + loc5m = loc5a.merge(loc5b, on="site_id", how="left", suffixes=(".x", ".y")) + loc5m["latitude"] = loc5m["Latitude.x"].combine_first(loc5m["Latitude.y"]) + loc5m["longitude"] = loc5m["Longitude.x"].combine_first(loc5m["Longitude.y"]) + loc5m["source_dataset_id"] = dsid(idx) + + loc5 = loc5m[["site_id", "latitude", "longitude", "source_dataset_id"]].drop_duplicates(subset=["site_id"]) + loc5["qc_flag"] = np.where(loc5["latitude"].isna() | loc5["longitude"].isna(), "g2", "g1") + __locations.append(loc5) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_06.py b/data/gold/expert_code/harmonize_sm/dataset_06.py new file mode 100644 index 0000000..3b6e3a3 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_06.py @@ -0,0 +1,69 @@ +"""Expert harmonization for dataset index 6. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 6, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 6 + # ============================================================= + idx = 6 + f6 = as_list(map_json[idx]["data_payload_files"])[0] + ddf6 = read_ds_csv(idx, f6) + + x = ddf6.iloc[2:].copy() + x["datetime_UTC"] = pd.to_datetime(x["TIMESTAMP"], format="%m/%d/%y %H:%M", errors="coerce", utc=True) + x["site_id"] = x["site"] + x["interval_min"] = interval_min(x["datetime_UTC"]) + + vwc_cols = [c for c in x.columns if "VWC" in c] + long = x.melt( + id_vars=["datetime_UTC", "site_id", "interval_min"], + value_vars=vwc_cols, + var_name="name", + value_name="volumetric_water_content_m3_m3", + ) + long["depth_m"] = ( + pd.to_numeric(long["name"].str.split("_").str[-1].str.replace("cm", "", regex=False), errors="coerce") / 100 + ) + long["replicate"] = 1 + long["is_timeseries"] = True + long["water_potential_kPa"] = np.nan + long["gravimetric_water_content_gH2O_gs"] = np.nan + + df6_harmonized = ensure_harmonized_cols(long) + __harmonized = df6_harmonized + __dataset_id = dsid(idx) + + loc6 = ( + ddf6.iloc[2:] + .groupby("site", as_index=False) + .first()[["site", "latitude", "longitude"]] + .rename(columns={"site": "site_id"}) + ) + loc6["source_dataset_id"] = dsid(idx) + loc6 = add_loc_qc(loc6) + __locations.append(loc6) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_07.py b/data/gold/expert_code/harmonize_sm/dataset_07.py new file mode 100644 index 0000000..5b49aa5 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_07.py @@ -0,0 +1,63 @@ +"""Expert harmonization for dataset index 7. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + parse_local_to_utc, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 7, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 7 + # ============================================================= + idx = 7 + f7 = as_list(map_json[idx]["data_payload_files"])[0] + m7 = as_list(map_json[idx]["location_metadata_files"])[0] + + ddf7 = read_ds_csv(idx, f7) + mdf7 = read_ds_csv(idx, m7, encoding='latin-1') + + x = ddf7.copy() + x["datetime_UTC"] = parse_local_to_utc(x["date.time"], "%m/%d/%y %H:%M", "America/Denver") + x["site_id"] = "BM" + x["interval_min"] = interval_min(x["datetime_UTC"]) + x["depth_m"] = pd.to_numeric(x["Depth (cm)"], errors="coerce") / 100 + x["replicate"] = 1 + x["is_timeseries"] = True + x["volumetric_water_content_m3_m3"] = pd.to_numeric(x["Volumetric Water Content"], errors="coerce") + x["water_potential_kPa"] = np.nan + x["gravimetric_water_content_gH2O_gs"] = np.nan + + df7_harmonized = ensure_harmonized_cols(x) + __harmonized = df7_harmonized + __dataset_id = dsid(idx) + + lat_col = [c for c in mdf7.columns if "Latitude" in c][0] + lon_col = [c for c in mdf7.columns if "Longitude" in c][0] + loc7 = mdf7.iloc[[0]].rename(columns={"Location": "site_id", lat_col: "latitude", lon_col: "longitude"})[ + ["site_id", "latitude", "longitude"] + ].copy() + loc7["source_dataset_id"] = dsid(idx) + loc7 = add_loc_qc(loc7) + __locations.append(loc7) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_08.py b/data/gold/expert_code/harmonize_sm/dataset_08.py new file mode 100644 index 0000000..954fc68 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_08.py @@ -0,0 +1,86 @@ +"""Expert harmonization for dataset index 8. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import re +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + parse_local_to_utc, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 8, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 8 + # ============================================================= + idx = 8 + f8 = as_list(map_json[idx]["data_payload_files"]) + m8 = as_list(map_json[idx]["location_metadata_files"]) + + ls8 = [read_ds_csv(idx, x) for x in f8] + mdf8 = pd.concat([read_ds_csv(idx, x) for x in m8], ignore_index=True) + + ls8_h = [] + for d in ls8: + cols = [c for c in d.columns if re.search(r"SoilMoisture|SoilMatricPot", c)] + x = d[["DateTime.MST"] + cols].copy() + x["datetime_UTC"] = parse_local_to_utc(x["DateTime.MST"], "%Y-%m-%d %H:%M:%S", "America/Denver") + x["interval_min"] = interval_min(x["datetime_UTC"]) + x = x.drop(columns=["DateTime.MST"]) + x.columns = [re.sub(r"\.m3\.m3|\.kPa", "", c) for c in x.columns] + + value_cols = [c for c in x.columns if re.search(r"SoilMoisture|SoilMatricPot", c)] + long = x.melt(id_vars=["datetime_UTC", "interval_min"], value_vars=value_cols, var_name="name", value_name="value") + m = long["name"].str.extract(r"^(SoilMoisture|SoilMatricPot)_(.*)$") + long["var"] = m[0] + long["depth_m"] = pd.to_numeric(m[1].str.replace("cm", "", regex=False), errors="coerce") / 100 + long["value"] = pd.to_numeric(long["value"], errors="coerce") + long.loc[long["value"] == -9999, "value"] = np.nan + + wide = ( + long.pivot_table( + index=["datetime_UTC", "interval_min", "depth_m"], + columns="var", + values="value", + aggfunc="first", + ) + .reset_index() + .rename(columns={"SoilMoisture": "volumetric_water_content_m3_m3", "SoilMatricPot": "water_potential_kPa"}) + ) + + wide["replicate"] = 1 + wide["site_id"] = "Slate River OBJ-2" + wide["is_timeseries"] = True + wide["gravimetric_water_content_gH2O_gs"] = np.nan + + ls8_h.append(ensure_harmonized_cols(wide)) + + df8_harmonized = pd.concat(ls8_h, ignore_index=True) + __harmonized = df8_harmonized + __dataset_id = dsid(idx) + + loc8 = mdf8.iloc[[0]].rename(columns={"SiteName": "site_id", "Latitude": "latitude", "Longitude": "longitude"})[ + ["site_id", "latitude", "longitude"] + ].copy() + loc8["source_dataset_id"] = dsid(idx) + loc8 = add_loc_qc(loc8) + __locations.append(loc8) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_09.py b/data/gold/expert_code/harmonize_sm/dataset_09.py new file mode 100644 index 0000000..968a203 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_09.py @@ -0,0 +1,77 @@ +"""Expert harmonization for dataset index 9. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + parse_local_to_utc, + interval_min, + utm32613_to_latlon, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 9, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 9 + # ============================================================= + idx = 9 + f9 = as_list(map_json[idx]["data_payload_files"])[0] + m9 = as_list(map_json[idx]["location_metadata_files"])[0] + + ddf9 = read_ds_csv(idx, f9) + mdf9 = read_ds_csv(idx, m9) + + x = ddf9.copy() + x["datetime_UTC"] = parse_local_to_utc(x["Collection Date"], "%Y-%m-%d", "America/Denver") + x["site_id"] = x["SampleSiteCode"] + x["interval_min"] = np.nan + x["depth_m"] = 0.2 + x["is_timeseries"] = False + x["water_potential_kPa"] = np.nan + + long = x.melt( + id_vars=["datetime_UTC", "site_id", "interval_min", "depth_m", "is_timeseries", "water_potential_kPa"], + value_vars=["VWC_1", "VWC_2"], + var_name="tmp", + value_name="VWC", + ) + long["replicate"] = long["tmp"].str.extract(r"_(\d+)")[0] + long["VWC"] = pd.to_numeric(long["VWC"], errors="coerce") + long.loc[long["VWC"] == -9999.0, "VWC"] = np.nan + long = long[long["VWC"].notna()].copy() + long["volumetric_water_content_m3_m3"] = long["VWC"] + long["gravimetric_water_content_gH2O_gs"] = np.nan + + df9_harmonized = ensure_harmonized_cols(long) + __harmonized = df9_harmonized + __dataset_id = dsid(idx) + + loc9_tmp = utm32613_to_latlon(mdf9, "Easting", "Northing") + loc9 = pd.DataFrame( + { + "site_id": loc9_tmp["SampleSiteCode"], + "latitude": loc9_tmp["latitude"], + "longitude": loc9_tmp["longitude"], + } + ) + loc9["source_dataset_id"] = dsid(idx) + loc9 = add_loc_qc(loc9) + __locations.append(loc9) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_10.py b/data/gold/expert_code/harmonize_sm/dataset_10.py new file mode 100644 index 0000000..2da60f7 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_10.py @@ -0,0 +1,70 @@ +"""Expert harmonization for dataset index 10. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import re +import numpy as np + +from common import ( + DatasetResult, + as_list, + parse_local_to_utc, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 10, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 10 + # ============================================================= + idx = 10 + f10 = as_list(map_json[idx]["data_payload_files"])[0] + m10 = as_list(map_json[REF_IDX]["location_metadata_files"])[0] + + ddf10 = read_ds_csv(idx, f10) + mdf10 = read_ds_csv(REF_IDX, m10) + + x = ddf10.iloc[1:].copy() + x["datetime_UTC"] = parse_local_to_utc(x["Date"], "%Y-%m-%d", "America/Denver") + x["interval_min"] = interval_min(x["datetime_UTC"]) + + vcols = [c for c in x.columns if re.search(r"vol_water_content", c)] + long = x.melt(id_vars=["datetime_UTC", "interval_min"], value_vars=vcols, var_name="name", value_name="volumetric_water_content_m3_m3") + nm = long["name"].str.extract(r"(.*)\._(.*)") + long["site_id"] = nm[0] + long["depth_m"] = np.select( + [long["site_id"].eq("PLM1"), long["site_id"].eq("PLM2"), long["site_id"].eq("PLM3")], + [0.3, 0.28, 0.2], + default=np.nan, + ) + long["is_timeseries"] = True + long["water_potential_kPa"] = np.nan + long["replicate"] = 1 + long["gravimetric_water_content_gH2O_gs"] = np.nan + + df10_harmonized = ensure_harmonized_cols(long) + __harmonized = df10_harmonized + __dataset_id = dsid(idx) + + sites = df10_harmonized["site_id"].dropna().astype(str).unique().tolist() + pattern = r"(?:^|)(%s)$" % "|".join([re.escape(s) for s in sites]) if sites else r"$^" + loc10 = mdf10[mdf10["Location_ID"].astype(str).str.contains(pattern, regex=True, na=False)].copy() + loc10["site_id"] = loc10["Location_ID"].str.replace("ER-", "", regex=False) + loc10 = loc10.rename(columns={"Latitude": "latitude", "Longitude": "longitude"})[["site_id", "latitude", "longitude"]] + loc10["source_dataset_id"] = dsid(idx) + loc10 = add_loc_qc(loc10) + __locations.append(loc10) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_15.py b/data/gold/expert_code/harmonize_sm/dataset_15.py new file mode 100644 index 0000000..03ae0a8 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_15.py @@ -0,0 +1,75 @@ +"""Expert harmonization for dataset index 15. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import re +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + interval_min, + utm32613_to_latlon, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 15, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 15 + # ============================================================= + idx = 15 + f15 = as_list(map_json[idx]["data_payload_files"])[0] + m15 = as_list(map_json[idx]["location_metadata_files"])[0] + + ddf15 = read_ds_csv(idx, f15) + mdf15 = read_ds_csv(idx, m15) + + x = ddf15.iloc[1:].copy() + x["datetime_UTC"] = pd.Timestamp("2019-07-02 06:00:00", tz="UTC") + x["is_timeseries"] = False + x["interval_min"] = np.nan + + vcols = [c for c in x.columns if re.search(r"SM\.VWC", c)] + long = x.melt( + id_vars=["datetime_UTC", "GPS_id", "is_timeseries", "interval_min"], + value_vars=vcols, + var_name="name", + value_name="volumetric_water_content_m3_m3", + ) + long["replicate"] = long["name"].str.extract(r"\._(.*)$")[0] + long["volumetric_water_content_m3_m3"] = pd.to_numeric(long["volumetric_water_content_m3_m3"], errors="coerce") / 100 + long["depth_m"] = 0.25 + long["water_potential_kPa"] = np.nan + long["gravimetric_water_content_gH2O_gs"] = np.nan + long["site_id"] = long["GPS_id"] + + df15_harmonized = ensure_harmonized_cols(long) + __harmonized = df15_harmonized + __dataset_id = dsid(idx) + + loc15_tmp = utm32613_to_latlon(mdf15, "Easting_m", "Northing_m") + loc15 = pd.DataFrame( + { + "site_id": loc15_tmp["GPS_id"], + "latitude": loc15_tmp["latitude"], + "longitude": loc15_tmp["longitude"], + } + ) + loc15["source_dataset_id"] = dsid(idx) + loc15 = add_loc_qc(loc15) + __locations.append(loc15) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_16.py b/data/gold/expert_code/harmonize_sm/dataset_16.py new file mode 100644 index 0000000..4a01b9d --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_16.py @@ -0,0 +1,92 @@ +"""Expert harmonization for dataset index 16. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import re +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + parse_local_to_utc, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 16, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 16 + # ============================================================= + idx = 16 + f16 = as_list(map_json[idx]["data_payload_files"]) + m16 = as_list(map_json[idx]["location_metadata_files"])[0] + s16 = as_list(map_json[idx]["sensor_metadata_files"])[0] + + ls16 = [read_ds_csv(idx, x) for x in f16] + mdf16 = read_ds_csv(idx, m16) + sdf16 = read_ds_csv(idx, s16) + + ls16_h = [] + for i, d in enumerate(ls16): + siten = f16[i] + x = d.copy() + x["datetime_UTC"] = parse_local_to_utc(x["TIMESTAMP_END"], "%Y-%m-%d %H:%M:%S", "Etc/GMT+7") + x["interval_min"] = interval_min(x["datetime_UTC"]) + parts = siten.split("_") + x["SITE_ID"] = parts[6] if len(parts) >= 7 else np.nan + + sw_cols = [c for c in x.columns if re.search(r"SWC|SWP", c)] + long = x.melt(id_vars=["datetime_UTC", "interval_min", "SITE_ID"], value_vars=sw_cols, var_name="VARIABLE", value_name="value") + + meta = sdf16.loc[sdf16["VARIABLE"].astype(str).str.contains("SWC|SWP", regex=True, na=False), ["SITE_ID", "VARIABLE", "HEIGHT"]] + long = long.merge(meta, on=["SITE_ID", "VARIABLE"], how="left") + long["VARIABLE"] = long["VARIABLE"].astype(str).str.split("_").str[0] + long["depth_m"] = pd.to_numeric(long["HEIGHT"], errors="coerce") * -1 + long["replicate"] = 1 + long["is_timeseries"] = True + long["gravimetric_water_content_gH2O_gs"] = np.nan + long = long.drop_duplicates() + + wide = ( + long.pivot_table( + index=["datetime_UTC", "SITE_ID", "depth_m", "replicate", "is_timeseries", "interval_min", "gravimetric_water_content_gH2O_gs"], + columns="VARIABLE", + values="value", + aggfunc="first", + ) + .reset_index() + ) + + wide["SWC"] = pd.to_numeric(wide.get("SWC"), errors="coerce") + wide["SWP"] = pd.to_numeric(wide.get("SWP"), errors="coerce") + wide["volumetric_water_content_m3_m3"] = np.where(wide["SWC"].isin([9999.0, -9999.0]), np.nan, wide["SWC"]) + wide["water_potential_kPa"] = np.where(wide["SWP"].isin([9999.0, -9999.0]), np.nan, wide["SWP"]) + wide = wide.rename(columns={"SITE_ID": "site_id"}) + + ls16_h.append(ensure_harmonized_cols(wide)) + + df16_harmonized = pd.concat(ls16_h, ignore_index=True) + __harmonized = df16_harmonized + __dataset_id = dsid(idx) + + loc16 = mdf16.rename(columns={"SITE_ID": "site_id", "LOCATION_LAT": "latitude", "LOCATION_LONG": "longitude"})[ + ["site_id", "latitude", "longitude"] + ].copy() + loc16["source_dataset_id"] = dsid(idx) + loc16 = add_loc_qc(loc16) + __locations.append(loc16) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_17.py b/data/gold/expert_code/harmonize_sm/dataset_17.py new file mode 100644 index 0000000..43e5ab3 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_17.py @@ -0,0 +1,80 @@ +"""Expert harmonization for dataset index 17. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import re +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + parse_local_to_utc, + interval_min, + utm32613_to_latlon, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 17, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 17 + # ============================================================= + idx = 17 + f17 = as_list(map_json[idx]["data_payload_files"]) + m17 = as_list(map_json[idx]["location_metadata_files"])[0] + + ls17 = [read_ds_csv(idx, x) for x in f17] + mdf17 = read_ds_csv(idx, m17) + + ls17_h = [] + for i, d in enumerate(ls17): + siten = f17[i] + x = d.iloc[1:].copy() + x["datetime_UTC"] = parse_local_to_utc(x["DateTime"], "%Y-%m-%d %H:%M:%S", "America/Denver") + x["interval_min"] = interval_min(x["datetime_UTC"]) + + site_guess = re.split(r"/|\.", siten) + x["site_id"] = site_guess[2] if len(site_guess) > 2 else np.nan + + mc_cols = [c for c in x.columns if re.search(r"MC", c)] + long = x.melt(id_vars=["datetime_UTC", "interval_min", "site_id"], value_vars=mc_cols, var_name="name", value_name="value") + m = long["name"].str.extract(r"(.*)_(.*)") + long["depth_m"] = pd.to_numeric(m[1].str.replace("m", "", regex=False), errors="coerce") + long["volumetric_water_content_m3_m3"] = pd.to_numeric(long["value"], errors="coerce") + long["replicate"] = 1 + long["is_timeseries"] = True + long["water_potential_kPa"] = np.nan + long["gravimetric_water_content_gH2O_gs"] = np.nan + + ls17_h.append(ensure_harmonized_cols(long)) + + df17_harmonized = pd.concat(ls17_h, ignore_index=True) + __harmonized = df17_harmonized + __dataset_id = dsid(idx) + + loc17_tmp = utm32613_to_latlon(mdf17, "Easting", "Northing") + loc17 = pd.DataFrame( + { + "site_id": loc17_tmp["ID"], + "latitude": loc17_tmp["latitude"], + "longitude": loc17_tmp["longitude"], + } + ) + loc17["source_dataset_id"] = dsid(idx) + loc17 = add_loc_qc(loc17) + loc17 = loc17[loc17["site_id"].astype(str).str.contains("TMC", na=False)] + __locations.append(loc17) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_18.py b/data/gold/expert_code/harmonize_sm/dataset_18.py new file mode 100644 index 0000000..aee29c6 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_18.py @@ -0,0 +1,71 @@ +"""Expert harmonization for dataset index 18. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + parse_local_to_utc, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 18, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 18 + # ============================================================= + idx = 18 + f18 = as_list(map_json[idx]["data_payload_files"])[0] + m18 = as_list(map_json[idx]["location_metadata_files"])[0] + + ddf18 = read_ds_csv(idx, f18, encoding="latin1") + mdf18 = read_ds_csv(idx, m18, encoding="latin1") + + x = ddf18.copy() + x["datetime_UTC"] = parse_local_to_utc(x["Date_Collected"], "%m/%d/%Y", "America/Denver") + x["interval_min"] = np.nan + x["site_id"] = (x["Field_Site"].astype(str) + "_" + x["Plot"].astype(str) + "_" + x["Topographic_Position"].astype(str)).str.replace(r"_$", "", regex=True) + x["depth_m"] = np.select( + [ + x["Depth_Increment"].eq("0-5 cm"), + x["Depth_Increment"].isin(["5-15 cm", "15-May"]), + x["Depth_Increment"].eq("15 cm +"), + ], + [0.025, 0.10, 0.15], + default=np.nan, + ) + x["is_timeseries"] = False + x["water_potential_kPa"] = np.nan + x["volumetric_water_content_m3_m3"] = np.nan + x["gravimetric_water_content_gH2O_gs"] = pd.to_numeric(x["Soil Water Content (g H2O per gram soil)"], errors="coerce") + x["replicate"] = x["Replicate"] + + df18_harmonized = ensure_harmonized_cols(x) + __harmonized = df18_harmonized + __dataset_id = dsid(idx) + + loc18 = mdf18.copy() + loc18["site_id"] = (loc18["Field_Site"].astype(str) + "_" + loc18["Plot"].astype(str) + "_" + loc18["Topographic_Position"].astype(str)).str.replace(r"_$", "", regex=True) + loc18["latitude"] = pd.to_numeric(loc18["Latitude"].astype(str).str.extract(r"([-+]?\d*\.?\d+)")[0], errors="coerce") + loc18["longitude"] = pd.to_numeric(loc18["Longitude"].astype(str).str.extract(r"([-+]?\d*\.?\d+)")[0], errors="coerce") + loc18["source_dataset_id"] = dsid(idx) + loc18 = loc18[["site_id", "latitude", "longitude", "source_dataset_id"]].drop_duplicates() + loc18 = add_loc_qc(loc18) + __locations.append(loc18) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_23.py b/data/gold/expert_code/harmonize_sm/dataset_23.py new file mode 100644 index 0000000..bf8ac86 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_23.py @@ -0,0 +1,87 @@ +"""Expert harmonization for dataset index 23. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + parse_local_to_utc, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 23, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 23 + # ============================================================= + idx = 23 + m23 = as_list(map_json[REF_IDX]["location_metadata_files"])[0] + + ddf23 = read_ds_csv(idx, "WM_SWC.csv") + sdf23 = read_ds_csv(idx, "sensor_metadata.csv") + pdf23 = read_ds_csv(idx, "plot_metadata.csv") + mdf23 = read_ds_csv(REF_IDX, m23) + + x = ddf23[ddf23["Sensor.Type"] == "SWC"].copy() + x["datetime_UTC"] = parse_local_to_utc(x["Date Time"], "%Y-%m-%d %H:%M:%S", "America/Denver") + + smeta = sdf23[sdf23["Sensor Type"] == "SWC"].copy() + smeta["Sensor.SN"] = pd.to_numeric(smeta["Sensor.SN"], errors="coerce") + smeta = smeta[["Sensor.SN", "Depth.cm"]] + + x["Sensor.SN"] = pd.to_numeric(x["Sensor.SN"], errors="coerce") + x = x.merge(smeta, on="Sensor.SN", how="left") + x = x.merge(pdf23, on="Plot.Location", how="left") + x = x[x["Treatment"] == "control"].copy() + + x["site_id"] = x["Point.Location"].astype(str) + "_" + x["Veg"].astype(str) + x["depth_m"] = pd.to_numeric(x["Depth.cm"], errors="coerce") / 100 + x["is_timeseries"] = True + x["volumetric_water_content_m3_m3"] = pd.to_numeric(x["Measurement"], errors="coerce") + x["water_potential_kPa"] = np.nan + x["gravimetric_water_content_gH2O_gs"] = np.nan + + x["rep_key"] = x["site_id"].astype(str) + "|" + x["Sensor.SN"].astype(str) + "|" + x["depth_m"].astype(str) + x["replicate"] = pd.factorize(x["rep_key"])[0] + 1 + + x = x.sort_values(["site_id", "depth_m", "replicate", "datetime_UTC"]) + x["interval_min"] = x.groupby(["site_id", "depth_m", "replicate"])["datetime_UTC"].diff().dt.total_seconds() / 60.0 + x.loc[x["interval_min"] < 0, "interval_min"] = np.nan + + df23_harmonized = ensure_harmonized_cols(x) + __harmonized = df23_harmonized + __dataset_id = dsid(idx) + + loc23 = pd.DataFrame({"site_id": df23_harmonized["site_id"].dropna().unique()}) + loc23["latitude"] = np.nan + loc23["longitude"] = np.nan + loc23 = loc23[loc23["site_id"].isin(df23_harmonized["site_id"].dropna().unique())].copy() + loc23["source_dataset_id"] = dsid(idx) + loc23 = add_loc_qc(loc23).drop_duplicates() + __locations.append(loc23) + + loc23 = mdf23.rename(columns={"Location_ID": "site_id", "Latitude": "latitude", "Longitude": "longitude"})[ + ["site_id", "latitude", "longitude"] + ].copy() + loc23 = loc23[loc23["site_id"].isin(df23_harmonized["site_id"].str.split('_').str[0].dropna().unique())].copy() + loc23["source_dataset_id"] = dsid(idx) + loc23 = add_loc_qc(loc23) + __locations.append(loc23) + + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_24.py b/data/gold/expert_code/harmonize_sm/dataset_24.py new file mode 100644 index 0000000..db38ef3 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_24.py @@ -0,0 +1,73 @@ +"""Expert harmonization for dataset index 24. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import re +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + parse_local_to_utc, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 24, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 24 + # ============================================================= + idx = 24 + f24 = as_list(map_json[idx]["data_payload_files"])[0] + m24 = as_list(map_json[idx]["location_metadata_files"])[0] + + ddf24 = read_ds_csv(idx, f24, sep=";") + mdf24 = read_ds_csv(idx, m24, sep=";") + + x = ddf24.copy() + keep = ["Timestamp"] + [c for c in x.columns if re.search(r"^P[12]_MP_(15|30|60)$", c)] + x = x[keep] + x["datetime_UTC"] = parse_local_to_utc(x["Timestamp"], "%Y-%m-%d", "America/Denver") + x["interval_min"] = interval_min(x["datetime_UTC"]) + x = x.drop(columns=["Timestamp"]) + + mp_cols = [c for c in x.columns if re.search(r"^P[12]_MP_(15|30|60)$", c)] + long = x.melt(id_vars=["datetime_UTC", "interval_min"], value_vars=mp_cols, var_name="name", value_name="MP") + m = long["name"].str.extract(r"^(P[12])_(MP)_(\d+)$") + long["site_id"] = m[0] + long["depth_cm"] = m[2] + long = long[long["MP"].notna()].copy() + long["depth_m"] = pd.to_numeric(long["depth_cm"], errors="coerce") / 100 + long["replicate"] = 1 + long["is_timeseries"] = True + long["volumetric_water_content_m3_m3"] = np.nan + long["gravimetric_water_content_gH2O_gs"] = np.nan + long["water_potential_kPa"] = pd.to_numeric(long["MP"], errors="coerce") + long = long[long["water_potential_kPa"].notna()] + + df24_harmonized = ensure_harmonized_cols(long) + __harmonized = df24_harmonized + __dataset_id = dsid(idx) + + loc24 = mdf24.rename(columns={"ID": "site_id", "Latitude": "latitude", "Longitude": "longitude"})[ + ["site_id", "latitude", "longitude"] + ].copy() + loc24 = loc24[loc24["site_id"].isin(df24_harmonized["site_id"].dropna().unique())].drop_duplicates() + loc24["source_dataset_id"] = dsid(idx) + loc24 = add_loc_qc(loc24) + __locations.append(loc24) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_25.py b/data/gold/expert_code/harmonize_sm/dataset_25.py new file mode 100644 index 0000000..7a0d4b7 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_25.py @@ -0,0 +1,84 @@ +"""Expert harmonization for dataset index 25. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import re +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + parse_local_to_utc, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 25, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 25 + # ============================================================= + idx = 25 + ddf25_conifer = read_ds_csv(idx, "Carbone_conifer.csv") + ddf25_aspen = read_ds_csv(idx, "Carbone_aspen.csv", header=None) + ddf25_aspen.columns = ddf25_conifer.columns + + ddf25 = pd.concat( + [ + ddf25_aspen.assign(site_id="aspen"), + ddf25_conifer.assign(site_id="conifer"), + ], + ignore_index=True, + ) + + x = ddf25.copy() + x["datetime_UTC"] = parse_local_to_utc( + x["Year"].astype(str) + "-" + x["Month"].astype(str) + "-" + x["Day"].astype(str) + "-" + x["Hour"].astype(str), + "%Y-%m-%d-%H", + "America/Denver", + ) + x = x.sort_values(["site_id", "datetime_UTC"]) + x["interval_min"] = x.groupby("site_id")["datetime_UTC"].diff().dt.total_seconds() / 60.0 + + vwc_cols = [c for c in x.columns if re.search(r"vwc", c, flags=re.IGNORECASE)] + long = x.melt( + id_vars=["datetime_UTC", "site_id", "interval_min"], + value_vars=vwc_cols, + var_name="name", + value_name="volumetric_water_content_m3_m3", + ) + long["depth_m"] = pd.to_numeric(long["name"].str.extract(r"(\d+)")[0], errors="coerce") / 100 + long["replicate"] = 1 + long["is_timeseries"] = True + long["water_potential_kPa"] = np.nan + long["gravimetric_water_content_gH2O_gs"] = np.nan + long["volumetric_water_content_m3_m3"] = pd.to_numeric(long["volumetric_water_content_m3_m3"], errors="coerce") + long = long[long["volumetric_water_content_m3_m3"].notna()] + + df25_harmonized = ensure_harmonized_cols(long) + __harmonized = df25_harmonized + __dataset_id = dsid(idx) + + loc25 = pd.DataFrame( + { + "site_id": ["aspen", "conifer"], + "latitude": [38.9581589, 38.9581589], + "longitude": [-106.9855848, -106.9855848], + } + ).drop_duplicates() + loc25["source_dataset_id"] = dsid(idx) + loc25 = add_loc_qc(loc25) + __locations.append(loc25) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_26.py b/data/gold/expert_code/harmonize_sm/dataset_26.py new file mode 100644 index 0000000..5df20e3 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_26.py @@ -0,0 +1,62 @@ +"""Expert harmonization for dataset index 26. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + parse_local_to_utc, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 26, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 26 + # ============================================================= + idx = 26 + f26 = as_list(map_json[idx]["data_payload_files"])[0] + m26 = as_list(map_json[REF_IDX]["location_metadata_files"])[0] + + ddf26 = read_ds_csv(idx, f26, index_col=0) + mdf26 = read_ds_csv(REF_IDX, m26) + + x = ddf26.copy() + x["datetime_UTC"] = parse_local_to_utc(x["Collection date"], "%m/%d/%y", "America/Denver") + x["site_id"] = x["SampleSiteCode"] + x["depth_m"] = (pd.to_numeric(x["Top sample depth_cm"], errors="coerce") + pd.to_numeric(x["Bottom sample depth_cm"], errors="coerce")) / 2 / 100 + x["replicate"] = 1 + x["is_timeseries"] = False + x["interval_min"] = np.nan + x["volumetric_water_content_m3_m3"] = pd.to_numeric(x["water content %vol"], errors="coerce") / 100 + x["gravimetric_water_content_gH2O_gs"] = np.nan + x["water_potential_kPa"] = np.nan + + df26_harmonized = ensure_harmonized_cols(x) + __harmonized = df26_harmonized + __dataset_id = dsid(idx) + + loc26 = mdf26.rename(columns={"Location_ID": "site_id", "Latitude": "latitude", "Longitude": "longitude"})[ + ["site_id", "latitude", "longitude"] + ].copy() + loc26 = loc26[loc26["site_id"].isin(df26_harmonized["site_id"].dropna().unique())] + loc26["source_dataset_id"] = dsid(idx) + loc26 = add_loc_qc(loc26) + __locations.append(loc26) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/dataset_27.py b/data/gold/expert_code/harmonize_sm/dataset_27.py new file mode 100644 index 0000000..5bae069 --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/dataset_27.py @@ -0,0 +1,108 @@ +"""Expert harmonization for dataset index 27. + +Auto-split from the expert monolith; the body below is verbatim except that +the shared-accumulator writes are captured into the returned DatasetResult. +""" +from __future__ import annotations + +import re +import numpy as np +import pandas as pd + +from common import ( + DatasetResult, + as_list, + interval_min, + ensure_harmonized_cols, + add_loc_qc, +) + + +def harmonize(ctx): + """Harmonize dataset 27, returning its DatasetResult.""" + map_json = ctx.map_json + REF_IDX = ctx.ref_idx + dsid = ctx.dsid + read_ds_csv = ctx.read_ds_csv + __locations = [] + + # ============================================================= + # Dataset 27 + # ============================================================= + idx = 27 + f27 = as_list(map_json[idx]["data_payload_files"]) + m27 = as_list(map_json[REF_IDX]["location_metadata_files"])[0] + + ls27 = [read_ds_csv(idx, x) for x in f27] + mdf27 = read_ds_csv(REF_IDX, m27) + + ls27_h = [] + for i, d in enumerate(ls27): + x = d.copy() + dt = pd.to_datetime(x["Time"], format="%m/%d/%Y %H:%M", errors="coerce") + dt = dt.dt.tz_localize("Etc/GMT+7", ambiguous="NaT", nonexistent="shift_forward") + x["datetime_UTC"] = dt.dt.tz_convert("UTC") + x["interval_min"] = interval_min(x["datetime_UTC"]) + + cols = [c for c in x.columns if re.search(r"^S[1-4]_wc_\(m3/m3\)$|^S5_wp_\(kPa\)$", c)] + long = x.melt(id_vars=["datetime_UTC", "interval_min"], value_vars=cols, var_name="name", value_name="value") + + long["variable"] = np.select( + [ + long["name"].str.contains("_wc_", regex=False, na=False), + long["name"].str.contains("_wp_", regex=False, na=False), + ], + ["volumetric_water_content_m3_m3", "water_potential_kPa"], + default=np.nan, + ) + + fname = f27[i] + long["depth_m"] = np.select( + [ + long["name"].str.contains(r"^S1_", regex=True, na=False) & pd.Series([bool(re.search(r"ER-PHS4", fname))] * len(long)), + long["name"].str.contains(r"^S1_", regex=True, na=False), + long["name"].str.contains(r"^S2_", regex=True, na=False), + long["name"].str.contains(r"^S3_", regex=True, na=False), + long["name"].str.contains(r"^S4_", regex=True, na=False), + long["name"].str.contains(r"^S5_", regex=True, na=False), + ], + [1.06, 1.15, 0.60, 0.30, 0.10, 0.30], + default=np.nan, + ) + + site_match = re.search(r"ER-PHS[1-4]", fname) + site_id_val = site_match.group(0) if site_match else np.nan + long["site_id"] = site_id_val + long["replicate"] = 1 + long["is_timeseries"] = True + + long = long[long["value"].notna() & long["variable"].notna() & long["depth_m"].notna()].copy() + wide = ( + long.pivot_table( + index=["datetime_UTC", "site_id", "depth_m", "replicate", "is_timeseries", "interval_min"], + columns="variable", + values="value", + aggfunc="first", + ) + .reset_index() + ) + + wide["volumetric_water_content_m3_m3"] = pd.to_numeric(wide.get("volumetric_water_content_m3_m3"), errors="coerce") + wide["water_potential_kPa"] = pd.to_numeric(wide.get("water_potential_kPa"), errors="coerce") + wide["gravimetric_water_content_gH2O_gs"] = np.nan + + ls27_h.append(ensure_harmonized_cols(wide)) + + df27_harmonized = pd.concat(ls27_h, ignore_index=True) + __harmonized = df27_harmonized + __dataset_id = dsid(idx) + + loc27 = mdf27.rename(columns={"Location_ID": "site_id", "Latitude": "latitude", "Longitude": "longitude"})[ + ["site_id", "latitude", "longitude"] + ].copy() + loc27 = loc27[loc27["site_id"].isin(df27_harmonized["site_id"].dropna().unique())] + loc27["source_dataset_id"] = dsid(idx) + loc27 = add_loc_qc(loc27) + __locations.append(loc27) + + return DatasetResult(__dataset_id, __harmonized, __locations) diff --git a/data/gold/expert_code/harmonize_sm/datasets.py b/data/gold/expert_code/harmonize_sm/datasets.py new file mode 100644 index 0000000..b7eb24e --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/datasets.py @@ -0,0 +1,20 @@ +"""Registry of per-dataset expert harmonizers. + +Each ``dataset_NN.py`` exposes ``harmonize(ctx) -> DatasetResult``. This maps a +dataset's mapping-JSON index to its harmonizer, so the runner can include or +omit whole datasets just by selecting keys — which is what makes leave-one- +cluster-out a matter of *not importing* a module rather than text-splicing the +monolith. + +The indices here are exactly the 19 datasets the expert harmonizes; the +reference dataset (0) and the globally-excluded datasets (11-14, 19-22) have no +module, so requesting them as a hold-out is rejected upstream. +""" +from __future__ import annotations + +from importlib import import_module + +# Datasets the expert harmonizes, in monolith order. +DATASET_INDICES = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 16, 17, 18, 23, 24, 25, 26, 27] + +DATASETS = {idx: import_module(f"dataset_{idx:02d}").harmonize for idx in DATASET_INDICES} diff --git a/data/gold/expert_code/harmonize_sm/run.py b/data/gold/expert_code/harmonize_sm/run.py new file mode 100644 index 0000000..007842c --- /dev/null +++ b/data/gold/expert_code/harmonize_sm/run.py @@ -0,0 +1,98 @@ +"""Run the expert harmonization over a chosen set of datasets. + +This is the modular replacement for running the monolith. Because each dataset +is an independent module that only *returns* its harmonized + location frames +(no shared-global accumulators), holding a cluster out is just:: + + python run.py --holdout 1,2,3,6,16,27 + +No cell-splicing needed: the held-out modules are simply not run, and the rest +stay runnable (the only cross-block dependency is the reference dataset 0, which +is never a hold-out target). + +Run from anywhere; the script puts its own directory on ``sys.path`` so the +sibling modules (``common``, ``datasets``, ``dataset_NN``) import cleanly. +""" +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parent)) + +from common import Context, harmonize_locations, write_outputs # noqa: E402 +from datasets import DATASET_INDICES, DATASETS # noqa: E402 + + +def run( + holdout: set[int] | None = None, + base_dir: Path | None = None, + out_dir: Path | None = None, + mapping_path: Path | None = None, + write: bool = True, +): + """Harmonize every dataset except those in ``holdout``. + + Returns ``(harmonized_ids, harmonized_data, loc_data)``. Mirrors the + monolith's end-to-end behavior: per-dataset frames are accumulated in + order, the location dedup/QA pass runs, and (when ``write``) the per-dataset + and concatenated-location CSVs are written. + """ + holdout = set(holdout or set()) + unknown = sorted(holdout - set(DATASET_INDICES)) + if unknown: + raise ValueError( + f"hold-out index/indices {unknown} have no dataset module; " + f"available datasets: {DATASET_INDICES}" + ) + + ctx = Context.load(base_dir=base_dir, out_dir=out_dir, mapping_path=mapping_path) + + results = [DATASETS[idx](ctx) for idx in DATASET_INDICES if idx not in holdout] + harmonized_ids = [r.dataset_id for r in results] + harmonized_data = [r.harmonized for r in results] + loc_data = [loc for r in results for loc in r.locations] + + # Dedup/QA pass (as in the expert script, its frame is reported, not written). + harmonize_locations(loc_data) + + if write: + write_outputs(ctx.out_dir, harmonized_ids, harmonized_data, loc_data) + + return harmonized_ids, harmonized_data, loc_data + + +def _parse_holdout(raw: str) -> set[int]: + return {int(tok) for tok in raw.split(",") if tok.strip()} + + +def main(argv: list[str] | None = None) -> None: + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument( + "--holdout", + default="", + help="Comma-separated dataset indices to leave out (e.g. '1,2,3,6,16,27').", + ) + parser.add_argument("--base-dir", type=Path, default=None, help="Override raw dataset base dir.") + parser.add_argument("--out-dir", type=Path, default=None, help="Override output dir.") + parser.add_argument("--mapping", type=Path, default=None, help="Override mapping JSON path.") + parser.add_argument( + "--no-write", action="store_true", help="Run without writing CSVs (smoke test)." + ) + args = parser.parse_args(argv) + + holdout = _parse_holdout(args.holdout) + ids, _, _ = run( + holdout=holdout, + base_dir=args.base_dir, + out_dir=args.out_dir, + mapping_path=args.mapping, + write=not args.no_write, + ) + kept = [i for i in DATASET_INDICES if i not in holdout] + print(f"harmonized {len(ids)} datasets (kept indices {kept}); held out {sorted(holdout)}") + + +if __name__ == "__main__": + main() diff --git a/skills/essdive_sm_curator/README.md b/skills/essdive_sm_curator/README.md index f1f7a0a..9930934 100644 --- a/skills/essdive_sm_curator/README.md +++ b/skills/essdive_sm_curator/README.md @@ -271,7 +271,7 @@ Test the two-skill workflow with a mixed DOI list: - Local metadata cache: `data/raw_cache/ess-dive_meta/` - Context dependencies: - data/gold/sm_data_harmonization_mapping.json # for schema reference and examples - - data/gold/expert_code/harmonize_ess-dive_soilmoisture_data.py # for code pattern reference + - data/gold/expert_code/harmonize_sm/ # modular expert harmonizer (common.py + dataset_NN.py) for code pattern reference - data/gold/harmonized_outputs/*.csv # harmonized datasets - Python libraries: `pandas`, `requests`, `json`, `re` diff --git a/skills/essdive_sm_harmonizer/SKILL.md b/skills/essdive_sm_harmonizer/SKILL.md index f7598cd..1cc6f66 100644 --- a/skills/essdive_sm_harmonizer/SKILL.md +++ b/skills/essdive_sm_harmonizer/SKILL.md @@ -11,7 +11,7 @@ metadata: updated: "2026-06-30" context_dependencies: - data/gold/sm_data_harmonization_mapping.json # for schema reference and examples - - data/gold/expert_code/harmonize_ess-dive_soilmoisture_data.py # for code pattern reference + - data/gold/expert_code/harmonize_sm/ # modular expert harmonizer (common.py + dataset_NN.py) for code pattern reference - data/gold/harmonized_outputs/*.csv # harmonized datasets usage: > Invoke this skill when adding a new ESS-DIVE soil moisture dataset to the diff --git a/src/folds/__init__.py b/src/folds/__init__.py index 1cd2684..8e5b7b8 100644 --- a/src/folds/__init__.py +++ b/src/folds/__init__.py @@ -1 +1 @@ -"""Leave-one-out fold preparation (e.g. ablating the expert monolith).""" +"""Leave-one-cluster-out fold preparation over the modular expert harmonizer.""" diff --git a/src/folds/ablate_monolith.py b/src/folds/ablate_monolith.py deleted file mode 100644 index 1adaefb..0000000 --- a/src/folds/ablate_monolith.py +++ /dev/null @@ -1,183 +0,0 @@ -"""Ablate per-dataset blocks from the expert harmonization monolith. - -The expert script ``data/gold/expert_code/harmonize_ess-dive_soilmoisture_data.py`` -is organized as a sequence of Jupyter-style ``# %%`` cells. Nineteen of those -cells are per-dataset harmonization blocks, each anchored by a single -``idx = N`` line whose ``N`` is the dataset's index in the gold mapping JSON. -Every other cell (config, shared helpers, the location-deduplication footer) -carries no ``idx =`` assignment. - -Grouped leave-one-cluster-out evaluation needs an *executable* copy of this -script with a chosen set of dataset blocks removed, so the held-out cluster is -absent from both the agent's exemplar context and any regenerated exemplar -outputs. The blocks only append to shared accumulators and the sole cross-block -dependency is ``REF_IDX = 0`` (index 0 is never a hold-out target), so dropping -a block's cell leaves the remaining script runnable. - -This module performs that splice by whole cell, keyed on ``idx = N``. Requiring -each requested hold-out to correspond to a real block automatically rejects -``REF_IDX`` (0) and the excluded datasets (11-14, 19-22), which carry no block — -so no coupling to ``cv_folds.yaml`` is needed. -""" -from __future__ import annotations - -import json -import re -from pathlib import Path -from typing import Optional - -import typer - -CELL_RE = re.compile(r"^# %%(\[markdown\])?\s*$") -IDX_RE = re.compile(r"^idx\s*=\s*(\d+)\s*$") - -DEFAULT_INPUT = Path("data/gold/expert_code/harmonize_ess-dive_soilmoisture_data.py") -DEFAULT_MAPPING = Path("data/gold/sm_data_harmonization_mapping.json") - - -def split_cells(source: str) -> list[str]: - """Split monolith source into ``# %%`` cells, losslessly. - - Each returned chunk begins with its ``# %%`` marker line, and concatenating - the chunks reproduces ``source`` exactly. - - >>> cells = split_cells("# %%\\na = 1\\n# %%\\nb = 2\\n") - >>> len(cells) - 2 - >>> "".join(cells) == "# %%\\na = 1\\n# %%\\nb = 2\\n" - True - """ - lines = source.splitlines(keepends=True) - cells: list[str] = [] - current: list[str] = [] - for line in lines: - if CELL_RE.match(line) and current: - cells.append("".join(current)) - current = [line] - else: - current.append(line) - if current: - cells.append("".join(current)) - return cells - - -def cell_index(cell: str) -> Optional[int]: - """Return the dataset index a cell harmonizes, or ``None`` for other cells. - - >>> cell_index("# %%\\nidx = 7\\nf7 = read(idx)\\n") - 7 - >>> cell_index("# %%\\nimport pandas as pd\\n") is None - True - """ - for line in cell.splitlines(): - m = IDX_RE.match(line) - if m: - return int(m.group(1)) - return None - - -def block_indices(source: str) -> list[int]: - """All dataset indices that have a block in ``source`` (sorted ascending).""" - found = (cell_index(c) for c in split_cells(source)) - return sorted(i for i in found if i is not None) - - -def ablate(source: str, holdout: set[int]) -> str: - """Return ``source`` with the held-out dataset blocks removed. - - Only whole ``# %%`` cells whose ``idx = N`` is in ``holdout`` are dropped; - every other cell (config, helpers, footer, other datasets) is preserved - verbatim, so the result stays executable. - - Raises ``ValueError`` if any requested index has no block in ``source`` — - this is what rejects ``REF_IDX`` (0) and the excluded datasets. - """ - present = set(block_indices(source)) - missing = sorted(holdout - present) - if missing: - raise ValueError( - f"no dataset block for index/indices {missing}; " - f"available blocks: {sorted(present)}" - ) - kept = [c for c in split_cells(source) if cell_index(c) not in holdout] - result = "".join(kept) - remaining = set(block_indices(result)) - expected = present - holdout - assert remaining == expected, ( - f"ablation removed the wrong blocks: kept {sorted(remaining)}, " - f"expected {sorted(expected)}" - ) - return result - - -def resolve_holdout(tokens: list[str], mapping_path: Optional[Path]) -> set[int]: - """Resolve hold-out tokens (integer indices or ``dataset_identifier``s). - - Integer-looking tokens are used directly; everything else is looked up as a - ``dataset_identifier`` in the gold mapping JSON. - """ - ident_to_idx: dict[str, int] = {} - if mapping_path is not None and mapping_path.exists(): - mapping = json.loads(mapping_path.read_text()) - ident_to_idx = { - e["dataset_identifier"]: e["index"] - for e in mapping - if e.get("dataset_identifier") is not None - } - out: set[int] = set() - for raw in tokens: - tok = raw.strip() - if not tok: - continue - if tok.isdigit(): - out.add(int(tok)) - elif tok in ident_to_idx: - out.add(ident_to_idx[tok]) - else: - raise ValueError( - f"hold-out token {tok!r} is neither an integer index nor a known " - f"dataset_identifier in {mapping_path}" - ) - return out - - -app = typer.Typer( - add_completion=False, - help="Ablate per-dataset blocks from the expert harmonization monolith.", -) - - -@app.command() -def main( - holdout: str = typer.Option( - ..., - "--holdout", - help="Comma-separated dataset indices or dataset_identifiers to remove.", - ), - input: Path = typer.Option( - DEFAULT_INPUT, "--input", "-i", help="Path to the expert monolith." - ), - output: Optional[Path] = typer.Option( - None, "--output", "-o", help="Write ablated script here (default: stdout)." - ), - mapping: Path = typer.Option( - DEFAULT_MAPPING, "--mapping", help="Gold mapping JSON for identifier resolution." - ), -) -> None: - """Write an executable copy of the monolith with the hold-out blocks removed.""" - source = input.read_text() - holdout_idx = resolve_holdout(holdout.split(","), mapping) - ablated = ablate(source, holdout_idx) - if output is None: - typer.echo(ablated, nl=False) - else: - output.write_text(ablated) - typer.echo( - f"wrote {output} (removed blocks {sorted(holdout_idx)}; " - f"{len(block_indices(ablated))} blocks remain)", - err=True, - ) - - -if __name__ == "__main__": - app() diff --git a/src/folds/expert_harmonizer.py b/src/folds/expert_harmonizer.py new file mode 100644 index 0000000..8069e73 --- /dev/null +++ b/src/folds/expert_harmonizer.py @@ -0,0 +1,212 @@ +"""Leave-one-cluster-out over the *modular* expert harmonizer. + +The expert harmonization used to be one 1200-line monolith, and grouped LOO was +done by text-splicing per-dataset ``# %%`` cells out of it (the old +``ablate_monolith``). The harmonizer is now a small package — a shared +``common`` library plus one ``dataset_NN.py`` per dataset — so holding a cluster +out no longer needs any splicing: + +* **Executable** (regenerate exemplar outputs without the held-out datasets): + :func:`run` calls the package runner, which simply doesn't run the held-out + modules. +* **Code context** (show the agent held-out-free reference code): + :func:`kept_module_paths` returns ``common.py`` plus the kept ``dataset_NN.py`` + files — the patterns to feed as exemplars. + +Requiring each hold-out to be a real dataset module automatically rejects the +reference dataset (0) and the globally-excluded datasets (11-14, 19-22), so +there is no coupling to ``cv_folds.yaml``. +""" +from __future__ import annotations + +import json +from pathlib import Path +from typing import Iterable, Optional + +import typer + +PACKAGE_DIR = Path("data/gold/expert_code/harmonize_sm") +DEFAULT_MAPPING = Path("data/gold/sm_data_harmonization_mapping.json") + +# Datasets the expert harmonizes, in module order. Kept in sync with +# ``harmonize_sm/datasets.py``; :func:`block_indices` validates it against disk. +DATASET_INDICES = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 16, 17, 18, 23, 24, 25, 26, 27] + + +def dataset_module_path(idx: int, package_dir: Path = PACKAGE_DIR) -> Path: + """Path to the per-dataset module for ``idx`` (whether or not it exists).""" + return package_dir / f"dataset_{idx:02d}.py" + + +def block_indices(package_dir: Path = PACKAGE_DIR) -> list[int]: + """Dataset indices that actually have a module on disk (sorted ascending).""" + found = [] + for idx in DATASET_INDICES: + if dataset_module_path(idx, package_dir).exists(): + found.append(idx) + return sorted(found) + + +def _validate_holdout(holdout: set[int], package_dir: Path) -> None: + present = set(block_indices(package_dir)) + missing = sorted(holdout - present) + if missing: + raise ValueError( + f"no dataset module for index/indices {missing}; " + f"available datasets: {sorted(present)}" + ) + + +def kept_indices(holdout: Iterable[int], package_dir: Path = PACKAGE_DIR) -> list[int]: + """Dataset indices that remain after removing ``holdout`` (sorted). + + Raises ``ValueError`` if any hold-out index has no module — this is what + rejects the reference dataset (0) and the excluded datasets (11-14, 19-22). + """ + holdout = set(holdout) + _validate_holdout(holdout, package_dir) + return [i for i in block_indices(package_dir) if i not in holdout] + + +def kept_module_paths(holdout: Iterable[int], package_dir: Path = PACKAGE_DIR) -> list[Path]: + """``common.py`` + the kept ``dataset_NN.py`` files, for exemplar code context. + + This is the modular replacement for the old "ablated monolith" text: the set + of reference code patterns the agent should see when the hold-out cluster is + removed. + """ + paths = [package_dir / "common.py"] + paths += [dataset_module_path(i, package_dir) for i in kept_indices(holdout, package_dir)] + return paths + + +def assemble_source(holdout: Iterable[int], package_dir: Path = PACKAGE_DIR) -> str: + """Concatenate the kept modules into one annotated text block (read-only). + + Convenience for tools that want a single string of the held-out-free expert + code (e.g. to drop into a prompt). It is reference material, not a runnable + script — execution goes through :func:`run`. + """ + chunks = [] + for path in kept_module_paths(holdout, package_dir): + chunks.append(f"# ===== {path.name} =====\n{path.read_text()}") + return "\n\n".join(chunks) + + +def resolve_holdout(tokens: Iterable[str], mapping_path: Optional[Path]) -> set[int]: + """Resolve hold-out tokens (integer indices or ``dataset_identifier``s). + + Integer-looking tokens are used directly; everything else is looked up as a + ``dataset_identifier`` in the gold mapping JSON. + """ + ident_to_idx: dict[str, int] = {} + if mapping_path is not None and Path(mapping_path).exists(): + mapping = json.loads(Path(mapping_path).read_text()) + ident_to_idx = { + e["dataset_identifier"]: e["index"] + for e in mapping + if e.get("dataset_identifier") is not None + } + out: set[int] = set() + for raw in tokens: + tok = raw.strip() + if not tok: + continue + if tok.isdigit(): + out.add(int(tok)) + elif tok in ident_to_idx: + out.add(ident_to_idx[tok]) + else: + raise ValueError( + f"hold-out token {tok!r} is neither an integer index nor a known " + f"dataset_identifier in {mapping_path}" + ) + return out + + +def run( + holdout: Iterable[int], + *, + base_dir: Optional[Path] = None, + out_dir: Optional[Path] = None, + mapping_path: Optional[Path] = None, + write: bool = True, +): + """Execute the modular harmonizer with ``holdout`` removed. + + Delegates to ``harmonize_sm/run.py``; imported lazily so this module stays + importable without the scientific stack. + """ + import importlib.util + import sys + + holdout = set(holdout) + _validate_holdout(holdout, PACKAGE_DIR) + + runner_path = PACKAGE_DIR / "run.py" + sys.path.insert(0, str(PACKAGE_DIR.resolve())) + spec = importlib.util.spec_from_file_location("harmonize_sm_run", runner_path) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module.run( + holdout=holdout, + base_dir=base_dir, + out_dir=out_dir, + mapping_path=mapping_path, + write=write, + ) + + +app = typer.Typer( + add_completion=False, + help="Leave-one-cluster-out over the modular expert harmonizer.", +) + + +@app.command("kept") +def cmd_kept( + holdout: str = typer.Option(..., "--holdout", help="Comma-separated indices or identifiers."), + mapping: Path = typer.Option(DEFAULT_MAPPING, "--mapping", help="Gold mapping JSON."), +) -> None: + """List the module files that remain after removing the hold-out datasets.""" + holdout_idx = resolve_holdout(holdout.split(","), mapping) + for path in kept_module_paths(holdout_idx): + typer.echo(str(path)) + + +@app.command("source") +def cmd_source( + holdout: str = typer.Option(..., "--holdout", help="Comma-separated indices or identifiers."), + output: Optional[Path] = typer.Option(None, "--output", "-o", help="Write here (default stdout)."), + mapping: Path = typer.Option(DEFAULT_MAPPING, "--mapping", help="Gold mapping JSON."), +) -> None: + """Emit the held-out-free expert code (common + kept datasets) as one block.""" + holdout_idx = resolve_holdout(holdout.split(","), mapping) + text = assemble_source(holdout_idx) + if output is None: + typer.echo(text, nl=False) + else: + output.write_text(text) + typer.echo( + f"wrote {output} (held out {sorted(holdout_idx)}; " + f"{len(kept_indices(holdout_idx))} datasets remain)", + err=True, + ) + + +@app.command("run") +def cmd_run( + holdout: str = typer.Option("", "--holdout", help="Comma-separated indices or identifiers."), + mapping: Path = typer.Option(DEFAULT_MAPPING, "--mapping", help="Gold mapping JSON."), + no_write: bool = typer.Option(False, "--no-write", help="Run without writing CSVs."), +) -> None: + """Run the harmonizer with the hold-out datasets removed.""" + holdout_idx = resolve_holdout(holdout.split(","), mapping) if holdout.strip() else set() + ids, _, _ = run(holdout_idx, write=not no_write) + typer.echo( + f"harmonized {len(ids)} datasets; held out {sorted(holdout_idx)}", err=True + ) + + +if __name__ == "__main__": + app() diff --git a/tests/test_ablate_monolith.py b/tests/test_ablate_monolith.py deleted file mode 100644 index 47aed19..0000000 --- a/tests/test_ablate_monolith.py +++ /dev/null @@ -1,137 +0,0 @@ -"""Tests for the monolith block-ablation splicer.""" -from __future__ import annotations - -import ast -from pathlib import Path - -import pytest -from typer.testing import CliRunner - -from src.folds.ablate_monolith import ( - app, - ablate, - block_indices, - cell_index, - resolve_holdout, - split_cells, -) - -# A miniature stand-in with the same shape as the expert script: a header cell, -# two dataset blocks anchored by ``idx = N``, an excluded-comment cell, and a -# footer cell. Keeps the test suite independent of the full corpus. -SYNTHETIC = """\ -# %% -acc = [] # header / shared accumulator - -# %% -# ============================================================= -# Dataset 1 -# ============================================================= -idx = 1 -acc.append(("ds1", idx)) - -# %% -# Dataset 2-3 excluded - -# %% -# ============================================================= -# Dataset 5 -# ============================================================= -idx = 5 -acc.append(("ds5", idx)) - -# %% -# footer -print(acc) -""" - - -def test_split_cells_is_lossless(): - assert "".join(split_cells(SYNTHETIC)) == SYNTHETIC - - -def test_block_indices_finds_only_dataset_cells(): - assert block_indices(SYNTHETIC) == [1, 5] - - -def test_cell_index_none_for_non_dataset_cell(): - assert cell_index("# %%\nimport os\n") is None - - -def test_ablate_removes_only_the_holdout_block(): - out = ablate(SYNTHETIC, {1}) - assert block_indices(out) == [5] - assert "ds1" not in out - assert "ds5" in out - assert "header" in out and "footer" in out # surrounding cells preserved - ast.parse(out) # still valid Python - - -def test_ablate_rejects_index_without_a_block(): - with pytest.raises(ValueError): - ablate(SYNTHETIC, {2}) # 2 is only an excluded comment, no block - - -@pytest.mark.parametrize( - "tokens,expected", - [(["1", "5"], {1, 5}), (["1"], {1}), ([" 5 ", ""], {5})], -) -def test_resolve_holdout_integer_tokens(tokens, expected): - assert resolve_holdout(tokens, None) == expected - - -def test_resolve_holdout_identifier_tokens(tmp_path): - mp = tmp_path / "mapping.json" - mp.write_text('[{"index": 5, "dataset_identifier": "ess-dive_abc"}]') - assert resolve_holdout(["ess-dive_abc"], mp) == {5} - - -def test_resolve_holdout_unknown_token_raises(tmp_path): - mp = tmp_path / "mapping.json" - mp.write_text("[]") - with pytest.raises(ValueError): - resolve_holdout(["not-an-index"], mp) - - -def test_cli_writes_ablated_file(tmp_path): - src = tmp_path / "mono.py" - src.write_text(SYNTHETIC) - out = tmp_path / "ablated.py" - result = CliRunner().invoke( - app, ["--holdout", "1", "--input", str(src), "--output", str(out)] - ) - assert result.exit_code == 0, result.output - text = out.read_text() - assert "ds1" not in text and "ds5" in text - - -# --- Against the real expert monolith, when present --- - -GOLD = Path("data/gold/expert_code/harmonize_ess-dive_soilmoisture_data.py") -EXPERT_TARGETS = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 16, 17, 18, 23, 24, 25, 26, 27] -needs_gold = pytest.mark.skipif(not GOLD.exists(), reason="gold expert code not present") - - -@needs_gold -def test_real_monolith_has_19_dataset_blocks(): - assert block_indices(GOLD.read_text()) == EXPERT_TARGETS - - -@needs_gold -def test_real_ablation_of_cluster_stays_parseable_and_keeps_scaffold(): - src = GOLD.read_text() - holdout = {1, 2, 3, 6, 16, 27} # cluster_1 "sg_ph_micromet" - out = ablate(src, holdout) - assert set(block_indices(out)) == set(EXPERT_TARGETS) - holdout - ast.parse(out) # syntactically valid - assert "REF_IDX = 0" in out # shared reference dataset preserved - assert "Location deduplication" in out # footer preserved - assert "idx = 27" not in out # held-out anchor gone - - -@needs_gold -def test_real_ablation_rejects_ref_idx_and_excluded(): - src = GOLD.read_text() - for bad in ({0}, {11}, {19}): - with pytest.raises(ValueError): - ablate(src, bad) diff --git a/tests/test_expert_harmonizer.py b/tests/test_expert_harmonizer.py new file mode 100644 index 0000000..a412580 --- /dev/null +++ b/tests/test_expert_harmonizer.py @@ -0,0 +1,117 @@ +"""Tests for the modular expert harmonizer's leave-one-cluster-out helpers.""" +from __future__ import annotations + +import ast +from pathlib import Path + +import pytest +from typer.testing import CliRunner + +from src.folds.expert_harmonizer import ( + DATASET_INDICES, + app, + assemble_source, + block_indices, + kept_indices, + kept_module_paths, + resolve_holdout, +) + +PACKAGE = Path("data/gold/expert_code/harmonize_sm") +needs_pkg = pytest.mark.skipif(not PACKAGE.exists(), reason="expert package not present") + +# The expert harmonizes these 19 datasets; 0 (reference) and 11-14, 19-22 +# (globally excluded) have no module. +EXPERT_TARGETS = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 16, 17, 18, 23, 24, 25, 26, 27] +CLUSTER_1 = {1, 2, 3, 6, 16, 27} # "sg_ph_micromet" + + +# --- pure helpers (no scientific stack needed) --- + + +@pytest.mark.parametrize( + "tokens,expected", + [(["1", "5"], {1, 5}), (["1"], {1}), ([" 5 ", ""], {5})], +) +def test_resolve_holdout_integer_tokens(tokens, expected): + assert resolve_holdout(tokens, None) == expected + + +def test_resolve_holdout_identifier_tokens(tmp_path): + mp = tmp_path / "mapping.json" + mp.write_text('[{"index": 5, "dataset_identifier": "ess-dive_abc"}]') + assert resolve_holdout(["ess-dive_abc"], mp) == {5} + + +def test_resolve_holdout_unknown_token_raises(tmp_path): + mp = tmp_path / "mapping.json" + mp.write_text("[]") + with pytest.raises(ValueError): + resolve_holdout(["not-an-index"], mp) + + +# --- against the real modular package --- + + +@needs_pkg +def test_block_indices_match_expert_targets(): + assert block_indices() == EXPERT_TARGETS + assert DATASET_INDICES == EXPERT_TARGETS + + +@needs_pkg +def test_kept_indices_removes_only_the_holdout(): + assert kept_indices(CLUSTER_1) == [i for i in EXPERT_TARGETS if i not in CLUSTER_1] + + +@needs_pkg +@pytest.mark.parametrize("bad", [{0}, {11}, {19}, {1, 12}]) +def test_kept_indices_rejects_ref_and_excluded(bad): + with pytest.raises(ValueError): + kept_indices(bad) + + +@needs_pkg +def test_kept_module_paths_drop_holdout_keep_common_and_rest(): + paths = kept_module_paths(CLUSTER_1) + names = [p.name for p in paths] + assert names[0] == "common.py" # shared library always included + for i in CLUSTER_1: + assert f"dataset_{i:02d}.py" not in names + for i in set(EXPERT_TARGETS) - CLUSTER_1: + assert f"dataset_{i:02d}.py" in names + # every kept module is independently valid Python + for p in paths: + ast.parse(p.read_text()) + + +@needs_pkg +def test_assemble_source_excludes_holdout_code_keeps_rest(): + out = assemble_source(CLUSTER_1) + assert "common.py" in out + assert "dataset_27.py" not in out # held out + assert "idx = 27" not in out + assert "dataset_05.py" in out # kept + assert "idx = 5" in out + + +@needs_pkg +def test_cli_source_writes_holdout_free_code(tmp_path): + out = tmp_path / "kept.py" + result = CliRunner().invoke( + app, ["source", "--holdout", "1,2,3,6,16,27", "--mapping", "/nonexistent", "-o", str(out)] + ) + assert result.exit_code == 0, result.output + text = out.read_text() + assert "dataset_27.py" not in text and "dataset_05.py" in text + + +@needs_pkg +def test_cli_kept_lists_remaining_modules(): + result = CliRunner().invoke( + app, ["kept", "--holdout", "1,2,3,6,16,27", "--mapping", "/nonexistent"] + ) + assert result.exit_code == 0, result.output + assert "common.py" in result.output + assert "dataset_27.py" not in result.output + assert "dataset_18.py" in result.output