hysplit

HYSPLIT Reader

HYSPLITReader

Bases: GriddedReader

Source code in monetio/readers/hysplit.py

@register_reader("hysplit")
class HYSPLITReader(GriddedReader):
    def open_dataset(
        self,
        files: str | list[str],
        drange: list[datetime.datetime] | None = None,
        century: int | None = None,
        verbose: bool = False,
        sample_time_stamp: str = "start",
        check_grid: bool = True,
        lazy: bool = False,
        **kwargs: Any,
    ) -> xr.Dataset:
        """
        Reads HYSPLIT binary concentration (cdump) files.

        Parameters
        ----------
        files : Union[str, List[str]]
            File path(s), URL(s), or glob pattern.
        drange : List[datetime.datetime], optional
            Date range to filter, by default None.
        century : int, optional
            Century to use for 2-digit years (e.g. 2000), by default None.
        verbose : bool, optional
            Whether to print verbose output, by default False.
        sample_time_stamp : str, optional
            Time stamp to use ('start' or 'end'), by default "start".
        check_grid : bool, optional
            Whether to fix grid continuity, by default True.
        lazy : bool, optional
            Whether to use Dask for lazy loading, by default False.
        **kwargs : Any
            Additional arguments passed to the driver.

        Returns
        -------
        xr.Dataset
            The processed HYSPLIT dataset.
        """
        # Set up kwargs for read_method
        read_kwargs = {
            "drange": drange,
            "century": century,
            "verbose": verbose,
            "sample_time_stamp": sample_time_stamp,
            "check_grid": check_grid,
        }

        # If it is a single file, we can use open_dataset_hysplit directly
        # If it is multiple files, we use combine_dataset which we'll modernize
        # but better yet, let's use the XarrayDriver with our custom logic.

        # We override the driver.open call to support our specific multi-file logic
        # while still benefiting from FileUtility and potential future driver features.
        ds = self.driver.open(
            files,
            read_method=open_dataset_hysplit,
            lazy=lazy,
            preprocess=None,  # HYSPLIT handles its own preprocessing
            **read_kwargs,
            **kwargs,
        )

        ds = update_history(ds, "Read HYSPLIT data.")
        return ds

    def harmonize(self, ds):
        return ds

open_dataset(files, drange=None, century=None, verbose=False, sample_time_stamp='start', check_grid=True, lazy=False, **kwargs)

Reads HYSPLIT binary concentration (cdump) files.

Parameters:

Name	Type	Description	Default
files	Union[str, List[str]]	File path(s), URL(s), or glob pattern.	required
drange	List[datetime]	Date range to filter, by default None.	None
century	int	Century to use for 2-digit years (e.g. 2000), by default None.	None
verbose	bool	Whether to print verbose output, by default False.	False
sample_time_stamp	str	Time stamp to use ('start' or 'end'), by default "start".	'start'
check_grid	bool	Whether to fix grid continuity, by default True.	True
lazy	bool	Whether to use Dask for lazy loading, by default False.	False
**kwargs	Any	Additional arguments passed to the driver.	{}

Returns:

Type	Description
Dataset	The processed HYSPLIT dataset.

Source code in monetio/readers/hysplit.py

def open_dataset(
    self,
    files: str | list[str],
    drange: list[datetime.datetime] | None = None,
    century: int | None = None,
    verbose: bool = False,
    sample_time_stamp: str = "start",
    check_grid: bool = True,
    lazy: bool = False,
    **kwargs: Any,
) -> xr.Dataset:
    """
    Reads HYSPLIT binary concentration (cdump) files.

    Parameters
    ----------
    files : Union[str, List[str]]
        File path(s), URL(s), or glob pattern.
    drange : List[datetime.datetime], optional
        Date range to filter, by default None.
    century : int, optional
        Century to use for 2-digit years (e.g. 2000), by default None.
    verbose : bool, optional
        Whether to print verbose output, by default False.
    sample_time_stamp : str, optional
        Time stamp to use ('start' or 'end'), by default "start".
    check_grid : bool, optional
        Whether to fix grid continuity, by default True.
    lazy : bool, optional
        Whether to use Dask for lazy loading, by default False.
    **kwargs : Any
        Additional arguments passed to the driver.

    Returns
    -------
    xr.Dataset
        The processed HYSPLIT dataset.
    """
    # Set up kwargs for read_method
    read_kwargs = {
        "drange": drange,
        "century": century,
        "verbose": verbose,
        "sample_time_stamp": sample_time_stamp,
        "check_grid": check_grid,
    }

    # If it is a single file, we can use open_dataset_hysplit directly
    # If it is multiple files, we use combine_dataset which we'll modernize
    # but better yet, let's use the XarrayDriver with our custom logic.

    # We override the driver.open call to support our specific multi-file logic
    # while still benefiting from FileUtility and potential future driver features.
    ds = self.driver.open(
        files,
        read_method=open_dataset_hysplit,
        lazy=lazy,
        preprocess=None,  # HYSPLIT handles its own preprocessing
        **read_kwargs,
        **kwargs,
    )

    ds = update_history(ds, "Read HYSPLIT data.")
    return ds

add_species(dset, species=None)

Sum multiple species into a single DataArray/Dataset.

Parameters:

Name	Type	Description	Default
dset	Dataset	Input HYSPLIT dataset.	required
species	list[str]	List of species to sum. If None, all species in 'Species ID' attribute are used.	None

Returns:

Type	Description
Dataset	Dataset with the summed species.

Source code in monetio/readers/hysplit.py

def add_species(dset: xr.Dataset, species: list[str] = None) -> xr.Dataset:
    """
    Sum multiple species into a single DataArray/Dataset.

    Parameters
    ----------
    dset : xr.Dataset
        Input HYSPLIT dataset.
    species : list[str], optional
        List of species to sum. If None, all species in 'Species ID' attribute are used.

    Returns
    -------
    xr.Dataset
        Dataset with the summed species.
    """
    splist = dset.attrs.get("Species ID", [])
    if not species:
        species = splist

    sflist = [s for s in species if s in dset.data_vars]

    if not sflist:
        return dset

    # Vectorized sum across selected species
    total_par = dset[sflist].to_array(dim="species").sum(dim="species")

    # Re-wrap in Dataset to maintain consistency with other readers
    res = total_par.to_dataset(name="_".join(sflist) if len(sflist) < 3 else "summed_species")

    # Transfer attributes
    res.attrs = dset.attrs.copy()
    res.attrs["Species ID"] = sflist
    return update_history(res, f"Added species sum: {sflist}")

check_grid_continuity(dset)

Check if the grid indices x and y are continuous (step of 1).

Parameters:

Name	Type	Description	Default
dset	Dataset	Input dataset.	required

Returns:

Type	Description
bool	True if grid is continuous.

Source code in monetio/readers/hysplit.py

def check_grid_continuity(dset: xr.Dataset) -> bool:
    """
    Check if the grid indices x and y are continuous (step of 1).

    Parameters
    ----------
    dset : xr.Dataset
        Input dataset.

    Returns
    -------
    bool
        True if grid is continuous.
    """
    # Use diff() for backend-agnostic continuity check
    if "x" in dset.dims and dset.x.size > 1:
        if not (dset.x.diff("x") == 1).all():
            return False
    if "y" in dset.dims and dset.y.size > 1:
        if not (dset.y.diff("y") == 1).all():
            return False
    return True

fix_grid_continuity(dset)

Fix grid continuity by reindexing to a full integer range.

Parameters:

Name	Type	Description	Default
dset	Dataset	Input HYSPLIT dataset.	required

Returns:

Type	Description
Dataset	Dataset with continuous grid.

Source code in monetio/readers/hysplit.py

def fix_grid_continuity(dset: xr.Dataset) -> xr.Dataset:
    """
    Fix grid continuity by reindexing to a full integer range.

    Parameters
    ----------
    dset : xr.Dataset
        Input HYSPLIT dataset.

    Returns
    -------
    xr.Dataset
        Dataset with continuous grid.
    """
    if not dset:
        return dset
    if check_grid_continuity(dset):
        return dset

    # Use min/max to avoid .values where possible (triggers 0-d compute if dask)
    x_min, x_max = int(dset.x.min()), int(dset.x.max())
    y_min, y_max = int(dset.y.min()), int(dset.y.max())

    x_new = np.arange(x_min, x_max + 1)
    y_new = np.arange(y_min, y_max + 1)

    # Reindex to the full range, filling gaps with 0
    dset = dset.reindex(x=x_new, y=y_new, fill_value=0)

    # Update lat/lon coordinates for the new grid
    mgrid = get_latlongrid(dset.attrs, x_new, y_new)
    dset = dset.assign_coords(latitude=(("y", "x"), mgrid[1]), longitude=(("y", "x"), mgrid[0]))

    dset = update_history(dset, "Fixed grid continuity and updated coordinates.")

    return dset

get_latlongrid(attrs, xindx, yindx)

Generate 2D latitude and longitude grids from HYSPLIT attributes and indices.

Parameters:

Name	Type	Description	Default
attrs	dict	HYSPLIT grid attributes.	required
xindx	ndarray	X-indices (1-based).	required
yindx	ndarray	Y-indices (1-based).	required

Returns:

Type	Description
list[ndarray]	[longitude_2d, latitude_2d]

Source code in monetio/readers/hysplit.py

def get_latlongrid(attrs: dict, xindx: np.ndarray, yindx: np.ndarray) -> list[np.ndarray]:
    """
    Generate 2D latitude and longitude grids from HYSPLIT attributes and indices.

    Parameters
    ----------
    attrs : dict
        HYSPLIT grid attributes.
    xindx : np.ndarray
        X-indices (1-based).
    yindx : np.ndarray
        Y-indices (1-based).

    Returns
    -------
    list[np.ndarray]
        [longitude_2d, latitude_2d]
    """
    xindx = np.asanyarray(xindx)
    yindx = np.asanyarray(yindx)
    if (xindx <= 0).any() or (yindx <= 0).any():
        raise ValueError("HYSPLIT grid error: indices must be > 0")

    lat_full, lon_full = getlatlon(attrs)

    # Vectorized indexing
    lon_sub = lon_full[xindx - 1]
    lat_sub = lat_full[yindx - 1]

    # Use xarray broadcasting for lazy 2D grid generation
    lon_2d, lat_2d = xr.broadcast(
        xr.DataArray(lon_sub, dims="x", coords={"x": xindx}),
        xr.DataArray(lat_sub, dims="y", coords={"y": yindx}),
    )

    # Return as numpy-like data to match expected signature
    return [lon_2d.transpose("y", "x").data, lat_2d.transpose("y", "x").data]

get_thickness(xrash)

Calculate layer thicknesses from vertical coordinates backend-agnostic.

Parameters:

Name	Type	Description	Default
xrash	Dataset | DataArray	Dataset containing vertical dimension 'z'.	required

Returns:

Type	Description
DataArray	Thickness of each layer.

Source code in monetio/readers/hysplit.py

def get_thickness(xrash: xr.Dataset | xr.DataArray) -> xr.DataArray:
    """
    Calculate layer thicknesses from vertical coordinates backend-agnostic.

    Parameters
    ----------
    xrash : xr.Dataset | xr.DataArray
        Dataset containing vertical dimension 'z'.

    Returns
    -------
    xr.DataArray
        Thickness of each layer.
    """
    # Vectorized approach: thickness = z[i] - z[i-1], where z[-1] = 0
    # This works for both deposition-inclusive (z[0]=0 -> thickness[0]=0)
    # and above-ground (z[0]>0 -> thickness[0]=z[0]) grids.
    z = xrash.z
    # We use shift(fill_value=0) to avoid xr.concat which can be expensive/tricky with indexes
    z_prev = z.shift(z=1, fill_value=0.0)
    delta = z - z_prev
    return delta.rename("thickness")

getlatlon(attrs)

Generate 1D latitude and longitude arrays from HYSPLIT attributes.

Parameters:

Name	Type	Description	Default
attrs	dict	HYSPLIT grid attributes.	required

Returns:

Type	Description
tuple[ndarray, ndarray]	(latitude, longitude)

Source code in monetio/readers/hysplit.py

def getlatlon(attrs: dict) -> tuple[np.ndarray, np.ndarray]:
    """
    Generate 1D latitude and longitude arrays from HYSPLIT attributes.

    Parameters
    ----------
    attrs : dict
        HYSPLIT grid attributes.

    Returns
    -------
    tuple[np.ndarray, np.ndarray]
        (latitude, longitude)
    """
    lon_tolerance = 0.001
    llcrnr_lat = attrs["llcrnr latitude"]
    llcrnr_lon = attrs["llcrnr longitude"]
    nlat = attrs["Number Lat Points"]
    nlon = attrs["Number Lon Points"]
    dlat = attrs["Latitude Spacing"]
    dlon = attrs["Longitude Spacing"]

    # Vectorized generation
    lat = llcrnr_lat + np.arange(nlat) * dlat
    lon = llcrnr_lon + np.arange(nlon) * dlon

    # Vectorized wrap-around
    lon = np.where(lon >= 180 + lon_tolerance, lon - 360, lon)

    return lat, lon

mass_loading(xrash, delta=None)

Calculate mass loading by vertically integrating concentration lazily.

Parameters:

Name	Type	Description	Default
xrash	DataArray | Dataset	Input data with concentration.	required
delta	DataArray | ndarray	Layer thicknesses. If None, calculated from 'z'.	None

Returns:

Type	Description
DataArray | Dataset	Mass loading (sum of conc * delta).

Source code in monetio/readers/hysplit.py

def mass_loading(
    xrash: xr.DataArray | xr.Dataset, delta: xr.DataArray | np.ndarray | None = None
) -> xr.DataArray | xr.Dataset:
    """
    Calculate mass loading by vertically integrating concentration lazily.

    Parameters
    ----------
    xrash : xr.DataArray | xr.Dataset
        Input data with concentration.
    delta : xr.DataArray | np.ndarray, optional
        Layer thicknesses. If None, calculated from 'z'.

    Returns
    -------
    xr.DataArray | xr.Dataset
        Mass loading (sum of conc * delta).
    """
    # 1. Exclude deposition layer for integration (sum over atmospheric layers only)
    xrash_no_dep = remove_dep(xrash)

    # 2. Get thicknesses
    if delta is None:
        weights = get_thickness(xrash_no_dep)
    else:
        if isinstance(delta, np.ndarray):
            # If provided as numpy, align with the dataset
            # We assume delta matches the original z including dep if lengths match
            if len(delta) == len(xrash.z):
                # We need a way to filter delta without .item()
                # But if it's already numpy, it's eager anyway.
                # However, for consistency, let's use xarray.
                full_weights = xr.DataArray(delta, coords={"z": xrash.z}, dims="z")
                weights = remove_dep(full_weights)
            else:
                weights = xr.DataArray(delta, coords={"z": xrash_no_dep.z}, dims="z")
        else:
            weights = remove_dep(delta)

    # 3. Compute lazy mass loading
    # Mask out non-positive weights to be safe
    weights = weights.where(weights > 0, np.nan)

    ml = (xrash_no_dep * weights).sum(dim="z", skipna=True)

    # 4. Provenance and scientific hygiene
    if isinstance(ml, xr.Dataset):
        ml = update_history(ml, "Calculated mass loading using standardized preprocessing.")
    elif isinstance(ml, xr.DataArray):
        # Ensure name is reasonable
        if hasattr(xrash, "name"):
            ml.name = f"{xrash.name}_mass_loading"
        # Update history if attributes are accessible
        if hasattr(ml, "attrs"):
            if "history" in xrash.attrs:
                ml.attrs["history"] = xrash.attrs["history"]
            ml = update_history(ml, "Calculated mass loading using standardized preprocessing.")

    return ml

remove_dep(xrash)

Mask the deposition layer (z=0) if present backend-agnostic. Keeps the same shape but replaces z=0 values with NaN to remain lazy.

Parameters:

Name	Type	Description	Default
xrash	Dataset | DataArray	Input data.	required

Returns:

Type	Description
Dataset | DataArray	Data with deposition layer masked.

Source code in monetio/readers/hysplit.py

def remove_dep(xrash: xr.Dataset | xr.DataArray) -> xr.Dataset | xr.DataArray:
    """
    Mask the deposition layer (z=0) if present backend-agnostic.
    Keeps the same shape but replaces z=0 values with NaN to remain lazy.

    Parameters
    ----------
    xrash : xr.Dataset | xr.DataArray
        Input data.

    Returns
    -------
    xr.Dataset | xr.DataArray
        Data with deposition layer masked.
    """
    return xrash.where(xrash.z > 0)

thickness_hash(xrash)

Map layer heights to their thicknesses.

Parameters:

Name	Type	Description	Default
xrash	Dataset | DataArray	Dataset containing vertical dimension 'z'.	required

Returns:

Type	Description
dict	Dictionary mapping height to thickness.

Source code in monetio/readers/hysplit.py

def thickness_hash(xrash: xr.Dataset | xr.DataArray) -> dict:
    """
    Map layer heights to their thicknesses.

    Parameters
    ----------
    xrash : xr.Dataset | xr.DataArray
        Dataset containing vertical dimension 'z'.

    Returns
    -------
    dict
        Dictionary mapping height to thickness.
    """
    delta = get_thickness(xrash)
    # Mapping requires discrete values, but we can do this without .values for dask-friendliness
    # by assuming coordinate 'z' is manageable in memory (usually < 100 levels)
    xlevs = xrash.z.data
    dhash = dict(zip(xlevs, delta.data))
    return dhash