CAMELSH

Overview

CAMELSH is the United States hourly hydrological dataset. Hourly resolution hydrological dataset for US catchments, providing high-temporal-resolution data for detailed hydrological analysis.

Dataset Information

• Region: United States
• Temporal Resolution: Hourly
• Module: hydrodataset.camelsh
• Class: Camelsh

Key Features

Hourly Resolution

Unlike daily CAMELS datasets, CAMELSH provides hourly timeseries data, enabling: - Sub-daily hydrological process analysis - Flash flood and storm event studies - High-frequency streamflow dynamics - Detailed precipitation event analysis

Static Attributes

Static catchment attributes include: - Basin area - Mean precipitation - Topographic characteristics - Land cover information - Soil properties - Climate indices

Dynamic Variables

Hourly timeseries variables available: - Streamflow (hourly) - Precipitation (hourly) - Temperature - Potential evapotranspiration - Solar radiation - And more...

Usage

Basic Usage

from hydrodataset.camelsh import Camelsh
from hydrodataset import resolve_data_path

# Initialize dataset
data_path = resolve_data_path("camelsh")
ds = Camelsh(data_path)

# Get basin IDs
basin_ids = ds.read_object_ids()
print(f"Number of basins: {len(basin_ids)}")

# Check available features
print("Static features:", ds.available_static_features)
print("Dynamic features:", ds.available_dynamic_features)

# Read hourly timeseries data
timeseries = ds.read_ts_xrdataset(
    gage_id_lst=basin_ids[:5],
    t_range=["2015-01-01", "2015-01-31"],  # One month of hourly data
    var_lst=["streamflow", "precipitation"]
)
print(timeseries)

# Read attribute data
attributes = ds.read_attr_xrdataset(
    gage_id_lst=basin_ids[:5],
    var_lst=["area", "p_mean"]
)
print(attributes)

Analyzing Storm Events

# Read hourly data for a specific storm event
storm_data = ds.read_ts_xrdataset(
    gage_id_lst=basin_ids[:3],
    t_range=["2015-06-15 00:00:00", "2015-06-20 23:00:00"],
    var_lst=["streamflow", "precipitation", "temperature_mean"]
)

# Analyze sub-daily patterns
import xarray as xr
hourly_precip = storm_data["precipitation"]
daily_total = hourly_precip.resample(time="1D").sum()
print("Daily precipitation totals:", daily_total)

Reading Specific Variables

# Read with specific time range (note hourly timestamps)
ts_data = ds.read_ts_xrdataset(
    gage_id_lst=basin_ids[:10],
    t_range=["2015-01-01 00:00:00", "2015-12-31 23:00:00"],
    var_lst=["streamflow", "precipitation", "temperature_mean"]
)

# Read basin area
areas = ds.read_area(gage_id_lst=basin_ids[:10])

# Read mean precipitation
mean_precip = ds.read_mean_prcp(gage_id_lst=basin_ids[:10])

Data Considerations

Large Data Volumes

Hourly data results in significantly larger datasets compared to daily data: - 24x more data points per day - Larger cache files - Longer initial cache generation time

Time Range Selection

When working with hourly data:

# Specify full datetime for hourly data
t_range = ["2015-01-01 00:00:00", "2015-01-31 23:00:00"]

# Or use date strings (defaults to 00:00:00)
t_range = ["2015-01-01", "2015-01-31"]

API Reference

hydrodataset.camelsh.Camelsh

Bases: HydroDataset

CAMELSH (CAMELS-Hourly) dataset class extending RainfallRunoff.

This class provides access to the CAMELSH dataset, which contains hourly hydrological and meteorological data for various watersheds.

Attributes:

Name	Type	Description
region		Geographic region identifier
download		Whether to download data automatically
ds_description		Dictionary containing dataset file paths

Source code in hydrodataset/camelsh.py

class Camelsh(HydroDataset):
    """CAMELSH (CAMELS-Hourly) dataset class extending RainfallRunoff.

    This class provides access to the CAMELSH dataset, which contains hourly
    hydrological and meteorological data for various watersheds.

    Attributes:
        region: Geographic region identifier
        download: Whether to download data automatically
        ds_description: Dictionary containing dataset file paths
    """

    def __init__(
        self, uri: str, region: Optional[str] = None, download: bool = False
    ) -> None:
        """Initialize CAMELSH dataset.

        Args:
            uri: Path to the data directory
            region: Geographic region identifier (optional)
            download: Whether to download data automatically (default: False)
        """
        super().__init__(uri)
        self.region = region
        self.download = download

        # aqua_fetch only supports local paths
        if not str(uri).startswith("s3://"):
            self.aqua_fetch = CAMELSH(uri)

    @property
    def _attributes_cache_filename(self):
        return "camelsh_attributes.nc"

    @property
    def _timeseries_cache_filename(self):
        return "camelsh_timeseries.nc"

    @property
    def default_t_range(self):
        return ["1980-01-01", "2024-12-31"]

    def read_object_ids(self) -> np.ndarray:
        """List station IDs from Hourly2/Hourly2/ directory (local or OSS)."""
        uri = str(self.data_source_dir).rstrip("/")
        h2_rel = "CAMELSH/Hourly2/Hourly2"
        if self._is_cloud():
            fs = self._make_s3fs()
            oss_path = f"{uri}/{h2_rel}".removeprefix("s3://")
            fnames = [p.split("/")[-1] for p in fs.ls(oss_path)]
        else:
            h2_dir = os.path.join(uri, *h2_rel.split("/"))
            fnames = os.listdir(h2_dir)
        stations = sorted(
            f.split("_")[0] for f in fnames if f.endswith("_hourly.nc")
        )
        return np.array(stations)

    def cache_attributes_to_zarr(self) -> None:
        """Read CAMELSH attribute CSV files from OSS and write attributes zarr to OSS."""
        import zarr

        fs = self._make_s3fs()
        uri = str(self.data_source_dir).rstrip("/")
        attr_dir = f"{uri}/CAMELSH/attributes/attributes"
        oss_attr_dir = attr_dir.removeprefix("s3://")

        # same rename as AquaFetch CAMELSH.static_map
        static_map = {
            "LAT_GAGE": "lat",
            "LNG_GAGE": "long",
            "ELEV_MEAN_M_BASIN": "elev_catch_m",
            "DRAIN_SQKM": "area_km2",
            "ELEV_SITE_M": "elev_gauge_m",
            "SLOPE_PCT": "slope_percent",
            "PDEN_2000_BLOCK": "pop_density_2000_km2",
            "PDEN_DAY_LANDSCAN_2007": "pop_density_2007_km2",
        }

        csv_paths = [f"s3://{p}" for p in fs.glob(f"{oss_attr_dir}/*.csv")]
        dfs = []
        for p in csv_paths:
            sep = "\t" if "attributes_hydroATLAS.csv" in p else ","
            with fs.open(p.removeprefix("s3://"), "rb") as fh:
                df = pd.read_csv(fh, index_col=0, sep=sep, dtype={0: str})
            df.index = df.index.astype(str)
            dfs.append(df)

        static = pd.concat(dfs, axis=1)
        static = static.loc[:, ~static.columns.duplicated()]
        static = static.rename(columns=static_map)
        static.columns = self._clean_feature_names(list(static.columns))
        static.index.name = "basin"

        ids = static.index.tolist()
        n = len(ids)
        zarr_name = self._attributes_cache_filename.replace(".nc", ".zarr")
        out, opts = self._zarr_path_and_opts(zarr_name)

        root = zarr.open_group(out, mode="w", storage_options=opts, zarr_format=2)
        for col in static.columns:
            vals = (
                static[col].values.astype(str)
                if static[col].dtype == object
                else static[col].values
            )
            arr = root.create_array(col, shape=(n,), chunks=(n,), dtype=vals.dtype)
            arr[:] = vals
            arr.attrs["_ARRAY_DIMENSIONS"] = ["basin"]
        basin_arr = root.create_array("basin", shape=(n,), chunks=(n,), dtype=str)
        basin_arr[:] = ids
        basin_arr.attrs["_ARRAY_DIMENSIONS"] = ["basin"]
        root.attrs["coordinates"] = "basin"
        self._write_zarr_units(root, "static")
        zarr.consolidate_metadata(root.store)
        print(f"Attributes zarr written to: {out}")

    _subclass_static_definitions = {
        # Basic station information
        "area": {"specific_name": "area_km2", "unit": "km^2"},
        "p_mean": {"specific_name": "p_mean", "unit": "mm/day"},
        "p_seasonality": {"specific_name": "p_seasonality", "unit": "none"},
        "frac_snow": {"specific_name": "frac_snow", "unit": "none"},
        "aridity": {"specific_name": "aridity_index", "unit": "none"},
    }
    _dynamic_variable_mapping = {
        StandardVariable.STREAMFLOW: {
            "default_source": "nldas",
            "sources": {"nldas": {"specific_name": "q_cms_obs", "unit": "m^3/s"}},
        },
        StandardVariable.PRECIPITATION: {
            "default_source": "nldas",
            "sources": {
                "nldas": {"specific_name": "pcp_mm", "unit": "mm"},
            },
        },
        StandardVariable.TEMPERATURE_MEAN: {
            "default_source": "nldas",
            "sources": {
                "nldas": {"specific_name": "airtemp_c_mean", "unit": "°C"},
            },
        },
        StandardVariable.LONGWAVE_SOLAR_RADIATION: {
            "default_source": "nldas",
            "sources": {
                "nldas": {"specific_name": "lwdown", "unit": "W/m^2"},
            },
        },
        # Shortwave radiation flux downwards (surface)
        StandardVariable.SOLAR_RADIATION: {
            "default_source": "nldas",
            "sources": {
                "nldas": {"specific_name": "swdown", "unit": "W/m^2"},
            },
        },
        # unit in aquafetch is mm/day.in paper is kg/m^2
        StandardVariable.POTENTIAL_EVAPOTRANSPIRATION: {
            "default_source": "nldas",
            "sources": {"nldas": {"specific_name": "pet_mm", "unit": "kg/m^2"}},
        },
        StandardVariable.SURFACE_PRESSURE: {
            "default_source": "nldas",
            "sources": {
                "nldas": {"specific_name": "psurf", "unit": "Pa"},
            },
        },
        # 10-meter above ground Zonal wind speed(east to west)
        StandardVariable.U_WIND_SPEED: {
            "default_source": "nldas",
            "sources": {
                "nldas": {"specific_name": "wind_e", "unit": "m/s"},
            },
        },
        # 10-meter above ground Meridional wind speed(north to south)
        StandardVariable.V_WIND_SPEED: {
            "default_source": "nldas",
            "sources": {
                "nldas": {"specific_name": "wind_n", "unit": "m/s"},
            },
        },
        StandardVariable.RELATIVE_HUMIDITY: {
            "default_source": "nldas",
            "sources": {
                "nldas": {"specific_name": "qair", "unit": "kg/kg"},
            },
        },
        StandardVariable.WATER_LEVEL: {
            "default_source": "nldas",
            "sources": {
                "nldas": {"specific_name": "water_level", "unit": "m"},
            },
        },
        StandardVariable.CAPE: {
            "default_source": "nldas",
            "sources": {
                "nldas": {"specific_name": "cape", "unit": "J/kg"},
            },
        },
        StandardVariable.CRAINF_FRAC: {
            "default_source": "nldas",
            "sources": {
                "nldas": {"specific_name": "crainf_frac", "unit": "Fraction"},
            },
        },
    }

    def cache_timeseries_to_zarr(self, batch_size: int = 1200) -> None:
        """Read CAMELSH NC files from OSS and write zarr to OSS.

        Reads q from Hourly2/Hourly2/{stn}_hourly.nc and forcing from
        timeseries_nonobs/Data/CAMELSH/timeseries_nonobs/{stn}.nc directly on
        OSS (no local copy needed). Resumable via _progress array.

        Time-sliced write: for each basin batch, the store is filled one time
        chunk at a time, so peak memory is ``batch_size x chunk_t x vars`` rather
        than ``batch_size x nt x vars``. This lets the basin chunk grow large
        (1200) on a memory-limited box; the cost is that each station NC is
        re-downloaded once per time chunk.
        """
        import gc
        import zarr
        import netCDF4 as nc4

        fs = self._make_s3fs()
        uri = str(self.data_source_dir).rstrip("/")
        h2_dir = f"{uri}/CAMELSH/Hourly2/Hourly2"
        nonobs_dir = f"{uri}/CAMELSH/timeseries_nonobs/Data/CAMELSH/timeseries_nonobs"
        ts_dir = f"{uri}/CAMELSH/timeseries/Data/CAMELSH/timeseries"

        # Same as AquaFetch CAMELSH.dyn_map
        dyn_map = {
            "Tair": "airtemp_C_mean",
            "PotEvap": "pet_mm",
            "Rainf": "pcp_mm",
            "streamflow": "q_cms_obs",
        }

        def open_nc_from_oss(s3_path: str) -> xr.Dataset:
            buf = fs.cat(s3_path.removeprefix("s3://"))
            nc = nc4.Dataset("inmemory", memory=buf)
            return xr.open_dataset(xr.backends.NetCDF4DataStore(nc))

        def oss_exists(s3_path: str) -> bool:
            return fs.exists(s3_path.removeprefix("s3://"))

        def rss_gb() -> float:
            # Current resident memory (GiB); Linux-only, -1 elsewhere.
            try:
                with open("/proc/self/statm") as f:
                    pages = int(f.read().split()[1])
                return pages * 4096 / 1024 ** 3
            except Exception:
                return -1.0

        # Canonical cleaned var list from _dynamic_variable_mapping
        cleaned_var_lst = [
            info["sources"][info["default_source"]]["specific_name"]
            for info in self._dynamic_variable_mapping.values()
        ]

        stations = self.read_object_ids().tolist()
        n = len(stations)
        times = pd.date_range(self.default_t_range[0], self.default_t_range[1], freq="h")
        nt = len(times)
        times_ns = times.asi8

        zarr_name = self._timeseries_cache_filename.replace(".nc", ".zarr")
        out, opts = self._zarr_path_and_opts(zarr_name)
        chunk_t = 24 * 365 * 9  # 9-year hourly chunk (~378 MiB with chunk_b=1200)
        chunk_b = min(batch_size, n)

        root = zarr.open_group(out, mode="a", storage_options=opts, zarr_format=2)

        if "basin" not in root:
            print(f"Pre-allocating zarr: {n} stations × {nt} timesteps × {len(cleaned_var_lst)} vars")
            for vn in cleaned_var_lst:
                arr = root.create_array(
                    vn, shape=(n, nt), chunks=(chunk_b, chunk_t),
                    dtype="float32", fill_value=np.nan,
                )
                arr.attrs["_ARRAY_DIMENSIONS"] = ["basin", "time"]
            time_arr = root.create_array("time", shape=(nt,), chunks=(chunk_t,), dtype="int64")
            time_arr[:] = times_ns
            time_arr.attrs["_ARRAY_DIMENSIONS"] = ["time"]
            time_arr.attrs["units"] = "nanoseconds since 1970-01-01"
            time_arr.attrs["calendar"] = "proleptic_gregorian"
            basin_arr = root.create_array("basin", shape=(n,), chunks=(n,), dtype=str)
            basin_arr[:] = stations
            basin_arr.attrs["_ARRAY_DIMENSIONS"] = ["basin"]
            prog = root.create_array("_progress", shape=(n,), chunks=(n,), dtype="int8", fill_value=0)
            prog[:] = 0
            prog.attrs["_ARRAY_DIMENSIONS"] = ["basin"]
            root.attrs["coordinates"] = "basin time"
            print("Pre-allocation done.")
        else:
            print(f"Resuming existing zarr store at {out}")

        progress = root["_progress"]
        n_batches = (n + batch_size - 1) // batch_size
        n_tchunks = (nt + chunk_t - 1) // chunk_t

        for batch_idx in range(0, n, batch_size):
            batch_end = min(batch_idx + batch_size, n)
            batch_num = batch_idx // batch_size + 1
            batch_stations = stations[batch_idx:batch_end]
            nb = len(batch_stations)

            if all(progress[batch_idx:batch_end]):
                print(f"Batch {batch_num}/{n_batches}: already done, skipping")
                continue

            print(f"Batch {batch_num}/{n_batches}: {nb} stations x "
                  f"{n_tchunks} time-chunks (re-read per chunk) ...")

            # Fill and write one time chunk at a time to keep peak memory at
            # nb x chunk_t x vars. Each station NC is re-downloaded per time chunk.
            for tc, t0 in enumerate(range(0, nt, chunk_t), start=1):
                t1 = min(t0 + chunk_t, nt)
                times_slice = times[t0:t1]
                slice_arrays = {
                    vn: np.full((nb, t1 - t0), np.nan, dtype=np.float32)
                    for vn in cleaned_var_lst
                }

                for i, stn in enumerate(batch_stations):
                    ds_q = ds_f = None
                    try:
                        # Read streamflow + water_level from Hourly2
                        ds_q = open_nc_from_oss(f"{h2_dir}/{stn}_hourly.nc")
                        ds_q = ds_q.rename({k: v for k, v in dyn_map.items() if k in ds_q.data_vars})
                        q_clean = {v: self._clean_feature_names([v])[0] for v in list(ds_q.data_vars)}
                        ds_q = ds_q.rename(q_clean)

                        # Read forcing from timeseries_nonobs (fallback: timeseries)
                        forcing_path = f"{nonobs_dir}/{stn}.nc"
                        if not oss_exists(forcing_path):
                            forcing_path = f"{ts_dir}/{stn}.nc"
                        ds_f = open_nc_from_oss(forcing_path)
                        ds_f = ds_f.drop_vars("Streamflow", errors="ignore")
                        if "DateTime" in ds_f.coords or "DateTime" in ds_f.dims:
                            ds_f = ds_f.rename({"DateTime": "time"})
                        ds_f = ds_f.rename({k: v for k, v in dyn_map.items() if k in ds_f.data_vars})
                        f_clean = {v: self._clean_feature_names([v])[0] for v in list(ds_f.data_vars)}
                        ds_f = ds_f.rename(f_clean)

                        for vn in cleaned_var_lst:
                            src = ds_q if vn in ds_q else (ds_f if vn in ds_f else None)
                            if src is not None:
                                da = src[vn].reindex(time=times_slice)
                                slice_arrays[vn][i] = da.values.astype(float)

                    except Exception as e:
                        print(f"  Station {stn} (tchunk {tc}): ERROR — {e}")
                        continue
                    finally:
                        # Always release the in-memory netCDF handles, even on error,
                        # or their C-level buffers accumulate and OOM the process.
                        if ds_q is not None:
                            ds_q.close()
                        if ds_f is not None:
                            ds_f.close()

                for vn in cleaned_var_lst:
                    root[vn][batch_idx:batch_end, t0:t1] = slice_arrays[vn]
                del slice_arrays
                gc.collect()  # reclaim closed netCDF buffers before the next chunk
                print(f"  Batch {batch_num}/{n_batches} tchunk {tc}/{n_tchunks} "
                      f"[{t0}:{t1}] written  (RSS={rss_gb():.1f} GiB)")

            progress[batch_idx:batch_end] = 1
            print(f"Batch {batch_num}/{n_batches}: done")

        self._write_zarr_units(root, "dynamic")
        zarr.consolidate_metadata(root.store)
        print(f"Timeseries zarr written to: {out}")

    def cache_timeseries_xrdataset(self, batch_size=100):
        """
        Cache timeseries data to NetCDF files in batches, each batch saved as a separate file

        Args:
            batch_size: Number of stations to process per batch, default is 100 stations
        """
        if not hasattr(self, "aqua_fetch"):
            raise NotImplementedError("aqua_fetch attribute is required")

        # Build mapping from variable names to units
        unit_lookup = {}
        if hasattr(self, "_dynamic_variable_mapping"):
            for std_name, mapping_info in self._dynamic_variable_mapping.items():
                for source, source_info in mapping_info["sources"].items():
                    unit_lookup[source_info["specific_name"]] = source_info["unit"]

        # Get all station IDs
        gage_id_lst = self.read_object_ids().tolist()
        total_stations = len(gage_id_lst)

        # Get original variable list and clean
        original_var_lst = self.aqua_fetch.dynamic_features
        cleaned_var_lst = self._clean_feature_names(original_var_lst)
        var_name_mapping = dict(zip(original_var_lst, cleaned_var_lst))

        print(
            f"Start batch processing {total_stations} stations, {batch_size} stations per batch"
        )
        print(
            f"Total number of batches: {(total_stations + batch_size - 1)//batch_size}"
        )

        # Ensure cache directory exists
        self.cache_dir.mkdir(parents=True, exist_ok=True)

        # Process stations in batches and save independently
        batch_num = 1
        for batch_idx in range(0, total_stations, batch_size):
            batch_end = min(batch_idx + batch_size, total_stations)
            batch_stations = gage_id_lst[batch_idx:batch_end]

            print(
                f"\nProcessing batch {batch_num}/{(total_stations + batch_size - 1)//batch_size}"
            )
            print(
                f"Station range: {batch_idx} - {batch_end-1} (total {len(batch_stations)} stations)"
            )

            try:
                # Get data for this batch
                batch_data = self.aqua_fetch.fetch_stations_features(
                    stations=batch_stations,
                    dynamic_features=original_var_lst,
                    static_features=None,
                    st=self.default_t_range[0],
                    en=self.default_t_range[1],
                    as_dataframe=False,
                )

                dynamic_data = (
                    batch_data[1] if isinstance(batch_data, tuple) else batch_data
                )

                # Process variables
                new_data_vars = {}
                time_coord = dynamic_data.coords["time"]

                for original_var in tqdm(
                    original_var_lst,
                    desc=f"Processing variables (batch {batch_num})",
                    total=len(original_var_lst),
                ):
                    cleaned_var = var_name_mapping[original_var]
                    var_data = []
                    for station in batch_stations:
                        if station in dynamic_data.data_vars:
                            station_data = dynamic_data[station].sel(
                                dynamic_features=original_var
                            )
                            if "dynamic_features" in station_data.coords:
                                station_data = station_data.drop("dynamic_features")
                            var_data.append(station_data)

                    if var_data:
                        combined = xr.concat(var_data, dim="basin")
                        combined["basin"] = batch_stations
                        combined.attrs["units"] = unit_lookup.get(
                            cleaned_var, "unknown"
                        )
                        new_data_vars[cleaned_var] = combined

                # Create Dataset for this batch
                batch_ds = xr.Dataset(
                    data_vars=new_data_vars,
                    coords={
                        "basin": batch_stations,
                        "time": time_coord,
                    },
                )

                # Save this batch to independent file
                batch_filename = f"batch{batch_num:03d}_camelsh_timeseries.nc"
                batch_filepath = self.cache_dir.joinpath(batch_filename)

                print(f"Saving batch {batch_num} to: {batch_filepath}")
                batch_ds.to_netcdf(batch_filepath)
                print(f"Batch {batch_num} saved successfully")

            except Exception as e:
                print(f"Batch {batch_num} processing failed: {e}")
                import traceback

                traceback.print_exc()
                continue

            batch_num += 1

        print(f"\nAll batches processed! Total {batch_num - 1} batch files saved")

    def read_ts_xrdataset(
        self,
        gage_id_lst: list = None,
        t_range: list = None,
        var_lst: list = None,
        sources: dict = None,
        **kwargs,
    ) -> xr.Dataset:
        """Read timeseries data from batch NC files (local) or zarr on OSS (cloud)."""
        if self._is_cloud():
            # Delegate to base: _load_ts_dataset opens the zarr, base handles
            # variable selection, time slicing, and renaming.
            return super().read_ts_xrdataset(
                gage_id_lst=gage_id_lst,
                t_range=t_range,
                var_lst=var_lst,
                sources=sources,
                **kwargs,
            )

        if (
            not hasattr(self, "_dynamic_variable_mapping")
            or not self._dynamic_variable_mapping
        ):
            raise NotImplementedError(
                "This dataset does not support the standardized variable mapping."
            )

        if var_lst is None:
            var_lst = list(self._dynamic_variable_mapping.keys())

        if t_range is None:
            t_range = self.default_t_range

        target_vars_to_fetch = []
        rename_map = {}

        # Process variable name mapping and data source selection
        for std_name in var_lst:
            if std_name not in self._dynamic_variable_mapping:
                raise ValueError(
                    f"'{std_name}' is not a recognized standard variable for this dataset."
                )

            mapping_info = self._dynamic_variable_mapping[std_name]

            # Determine which data source(s) to use
            is_explicit_source = sources and std_name in sources
            sources_to_use = []
            if is_explicit_source:
                provided_sources = sources[std_name]
                if isinstance(provided_sources, list):
                    sources_to_use.extend(provided_sources)
                else:
                    sources_to_use.append(provided_sources)
            else:
                sources_to_use.append(mapping_info["default_source"])

            # Only need suffix when user explicitly requests multiple data sources
            needs_suffix = is_explicit_source and len(sources_to_use) > 1
            for source in sources_to_use:
                if source not in mapping_info["sources"]:
                    raise ValueError(
                        f"Source '{source}' is not available for variable '{std_name}'."
                    )

                actual_var_name = mapping_info["sources"][source]["specific_name"]
                target_vars_to_fetch.append(actual_var_name)
                output_name = f"{std_name}_{source}" if needs_suffix else std_name
                rename_map[actual_var_name] = output_name

        # Find all batch files
        import glob

        batch_pattern = str(self.cache_dir / "batch*_camelsh_timeseries.nc")
        batch_files = sorted(glob.glob(batch_pattern))

        if not batch_files:
            print("No batch cache files found, starting cache creation...")
            self.cache_timeseries_xrdataset()
            batch_files = sorted(glob.glob(batch_pattern))

            if not batch_files:
                raise FileNotFoundError("Cache creation failed, no batch files found")

        print(f"Found {len(batch_files)} batch files")

        # If no stations specified, read all stations
        if gage_id_lst is None:
            print("No station list specified, will read all stations...")
            gage_id_lst = self.read_object_ids().tolist()

        # Convert station IDs to strings (ensure consistency)
        gage_id_lst = [str(gid) for gid in gage_id_lst]

        # Iterate through batch files to find batches containing required stations
        relevant_datasets = []
        for batch_file in batch_files:
            try:
                # First open only coordinates, don't load data
                ds_batch = xr.open_dataset(batch_file)
                batch_basins = [str(b) for b in ds_batch.basin.values]

                # Check if this batch contains required stations
                common_basins = list(set(gage_id_lst) & set(batch_basins))

                if common_basins:
                    print(
                        f"Batch {os.path.basename(batch_file)}: contains {len(common_basins)} required stations"
                    )

                    # Check if variables exist
                    missing_vars = [
                        v for v in target_vars_to_fetch if v not in ds_batch.data_vars
                    ]
                    if missing_vars:
                        ds_batch.close()
                        raise ValueError(
                            f"Batch {os.path.basename(batch_file)} missing variables: {missing_vars}"
                        )

                    # Select variables and stations
                    ds_subset = ds_batch[target_vars_to_fetch]
                    ds_selected = ds_subset.sel(
                        basin=common_basins, time=slice(t_range[0], t_range[1])
                    )

                    relevant_datasets.append(ds_selected)
                    ds_batch.close()
                else:
                    ds_batch.close()

            except Exception as e:
                print(f"Failed to read batch file {batch_file}: {e}")
                continue

        if not relevant_datasets:
            raise ValueError(
                f"Specified stations not found in any batch files: {gage_id_lst}"
            )

        print(f"Reading data from {len(relevant_datasets)} batches...")

        # Merge data from all relevant batches
        if len(relevant_datasets) == 1:
            final_ds = relevant_datasets[0]
        else:
            final_ds = xr.concat(relevant_datasets, dim="basin")

        # Rename to standard variable names
        final_ds = final_ds.rename(rename_map)

        # Ensure stations are arranged in input order
        if len(gage_id_lst) > 0:
            # Only select actually existing stations
            existing_basins = [b for b in gage_id_lst if b in final_ds.basin.values]
            if existing_basins:
                final_ds = final_ds.sel(basin=existing_basins)

        return final_ds

default_t_range property

__init__(uri, region=None, download=False)

Initialize CAMELSH dataset.

Parameters:

Name	Type	Description	Default
uri	str	Path to the data directory	required
region	Optional[str]	Geographic region identifier (optional)	None
download	bool	Whether to download data automatically (default: False)	False

Source code in hydrodataset/camelsh.py

def __init__(
    self, uri: str, region: Optional[str] = None, download: bool = False
) -> None:
    """Initialize CAMELSH dataset.

    Args:
        uri: Path to the data directory
        region: Geographic region identifier (optional)
        download: Whether to download data automatically (default: False)
    """
    super().__init__(uri)
    self.region = region
    self.download = download

    # aqua_fetch only supports local paths
    if not str(uri).startswith("s3://"):
        self.aqua_fetch = CAMELSH(uri)

read_object_ids()

List station IDs from Hourly2/Hourly2/ directory (local or OSS).

Source code in hydrodataset/camelsh.py

def read_object_ids(self) -> np.ndarray:
    """List station IDs from Hourly2/Hourly2/ directory (local or OSS)."""
    uri = str(self.data_source_dir).rstrip("/")
    h2_rel = "CAMELSH/Hourly2/Hourly2"
    if self._is_cloud():
        fs = self._make_s3fs()
        oss_path = f"{uri}/{h2_rel}".removeprefix("s3://")
        fnames = [p.split("/")[-1] for p in fs.ls(oss_path)]
    else:
        h2_dir = os.path.join(uri, *h2_rel.split("/"))
        fnames = os.listdir(h2_dir)
    stations = sorted(
        f.split("_")[0] for f in fnames if f.endswith("_hourly.nc")
    )
    return np.array(stations)

read_ts_xrdataset(gage_id_lst=None, t_range=None, var_lst=None, sources=None, **kwargs)

Read timeseries data from batch NC files (local) or zarr on OSS (cloud).

Source code in hydrodataset/camelsh.py

def read_ts_xrdataset(
    self,
    gage_id_lst: list = None,
    t_range: list = None,
    var_lst: list = None,
    sources: dict = None,
    **kwargs,
) -> xr.Dataset:
    """Read timeseries data from batch NC files (local) or zarr on OSS (cloud)."""
    if self._is_cloud():
        # Delegate to base: _load_ts_dataset opens the zarr, base handles
        # variable selection, time slicing, and renaming.
        return super().read_ts_xrdataset(
            gage_id_lst=gage_id_lst,
            t_range=t_range,
            var_lst=var_lst,
            sources=sources,
            **kwargs,
        )

    if (
        not hasattr(self, "_dynamic_variable_mapping")
        or not self._dynamic_variable_mapping
    ):
        raise NotImplementedError(
            "This dataset does not support the standardized variable mapping."
        )

    if var_lst is None:
        var_lst = list(self._dynamic_variable_mapping.keys())

    if t_range is None:
        t_range = self.default_t_range

    target_vars_to_fetch = []
    rename_map = {}

    # Process variable name mapping and data source selection
    for std_name in var_lst:
        if std_name not in self._dynamic_variable_mapping:
            raise ValueError(
                f"'{std_name}' is not a recognized standard variable for this dataset."
            )

        mapping_info = self._dynamic_variable_mapping[std_name]

        # Determine which data source(s) to use
        is_explicit_source = sources and std_name in sources
        sources_to_use = []
        if is_explicit_source:
            provided_sources = sources[std_name]
            if isinstance(provided_sources, list):
                sources_to_use.extend(provided_sources)
            else:
                sources_to_use.append(provided_sources)
        else:
            sources_to_use.append(mapping_info["default_source"])

        # Only need suffix when user explicitly requests multiple data sources
        needs_suffix = is_explicit_source and len(sources_to_use) > 1
        for source in sources_to_use:
            if source not in mapping_info["sources"]:
                raise ValueError(
                    f"Source '{source}' is not available for variable '{std_name}'."
                )

            actual_var_name = mapping_info["sources"][source]["specific_name"]
            target_vars_to_fetch.append(actual_var_name)
            output_name = f"{std_name}_{source}" if needs_suffix else std_name
            rename_map[actual_var_name] = output_name

    # Find all batch files
    import glob

    batch_pattern = str(self.cache_dir / "batch*_camelsh_timeseries.nc")
    batch_files = sorted(glob.glob(batch_pattern))

    if not batch_files:
        print("No batch cache files found, starting cache creation...")
        self.cache_timeseries_xrdataset()
        batch_files = sorted(glob.glob(batch_pattern))

        if not batch_files:
            raise FileNotFoundError("Cache creation failed, no batch files found")

    print(f"Found {len(batch_files)} batch files")

    # If no stations specified, read all stations
    if gage_id_lst is None:
        print("No station list specified, will read all stations...")
        gage_id_lst = self.read_object_ids().tolist()

    # Convert station IDs to strings (ensure consistency)
    gage_id_lst = [str(gid) for gid in gage_id_lst]

    # Iterate through batch files to find batches containing required stations
    relevant_datasets = []
    for batch_file in batch_files:
        try:
            # First open only coordinates, don't load data
            ds_batch = xr.open_dataset(batch_file)
            batch_basins = [str(b) for b in ds_batch.basin.values]

            # Check if this batch contains required stations
            common_basins = list(set(gage_id_lst) & set(batch_basins))

            if common_basins:
                print(
                    f"Batch {os.path.basename(batch_file)}: contains {len(common_basins)} required stations"
                )

                # Check if variables exist
                missing_vars = [
                    v for v in target_vars_to_fetch if v not in ds_batch.data_vars
                ]
                if missing_vars:
                    ds_batch.close()
                    raise ValueError(
                        f"Batch {os.path.basename(batch_file)} missing variables: {missing_vars}"
                    )

                # Select variables and stations
                ds_subset = ds_batch[target_vars_to_fetch]
                ds_selected = ds_subset.sel(
                    basin=common_basins, time=slice(t_range[0], t_range[1])
                )

                relevant_datasets.append(ds_selected)
                ds_batch.close()
            else:
                ds_batch.close()

        except Exception as e:
            print(f"Failed to read batch file {batch_file}: {e}")
            continue

    if not relevant_datasets:
        raise ValueError(
            f"Specified stations not found in any batch files: {gage_id_lst}"
        )

    print(f"Reading data from {len(relevant_datasets)} batches...")

    # Merge data from all relevant batches
    if len(relevant_datasets) == 1:
        final_ds = relevant_datasets[0]
    else:
        final_ds = xr.concat(relevant_datasets, dim="basin")

    # Rename to standard variable names
    final_ds = final_ds.rename(rename_map)

    # Ensure stations are arranged in input order
    if len(gage_id_lst) > 0:
        # Only select actually existing stations
        existing_basins = [b for b in gage_id_lst if b in final_ds.basin.values]
        if existing_basins:
            final_ds = final_ds.sel(basin=existing_basins)

    return final_ds