Explore Module

High-level analysis and visualization tools for climate data exploration.

Overview

The climakitae.explore module provides user-friendly classes and functions for climate data analysis, including:

Warming level analysis — Analyze data by global warming thresholds
Time series visualization — Plot multi-model ensemble trends
Uncertainty quantification — CMIP6 ensemble analysis
Threshold analysis — Climate threshold and exceedance analysis
Typical meteorological year — TMY analysis for engineering applications
Vulnerability assessment — Climate vulnerability analysis tools

Warming Levels

High-level container for the warming levels analysis workflow.

Orchestrates data retrieval, GWL-window slicing, optional anomaly computation, and optional in-memory loading for all available warming levels.

Attributes:

Name	Type	Description
`wl_params`	`WarmingLevelChoose`	Parameter panel that controls variable, resolution, window size, and other query options.
`warming_levels`	`list of float`	Warming level thresholds derived from the GWL lookup table.
`gwl_times`	`DataFrame or None`	Lookup table mapping each simulation/scenario combination to the calendar year it crosses each warming level. Populated by :py:meth:`calculate`.
`catalog_data`	`DataArray or None`	Raw stacked data retrieved from the catalog. Populated by :py:meth:`calculate`.
`sliced_data`	`dict of {str: xarray.DataArray}`	Mapping from warming level string (e.g. `"1.5"`) to the corresponding time-windowed DataArray. Populated by :py:meth:`calculate`.
`gwl_snapshots`	`DataArray`	Time-mean snapshots for every warming level, concatenated along a `warming_level` dimension. Populated by :py:meth:`calculate`.

Examples:

>>> wl = WarmingLevels()
>>> wl.wl_params.variable = "Maximum air temperature at 2m"
>>> wl.calculate()  # blocks until all levels are computed
>>> snapshot_2deg = wl.gwl_snapshots.sel(warming_level="2.0")

Source code in climakitae/explore/warming.py

def __init__(self):
    self.wl_params = WarmingLevelChoose()
    # self.warming_levels = ["0.8", "1.2", "1.5", "2.0", "3.0", "4.0"]
    self.warming_levels = _check_available_warming_levels()
    self.gwl_times = None  # Placeholder for the warming level times
    self.catalog_data = None  # Placeholder for the catalog data

`find_warming_slice(level, gwl_times)`

Slice catalog data to the window around a single warming level.

Groups catalog_data by all_sims, calls :py:func:get_sliced_data for every simulation, then cleans and relabels the result.

Parameters:

Name	Type	Description	Default
`level`	`str`	Warming level threshold to slice around, e.g. `"1.5"`.	required
`gwl_times`	`DataFrame`	GWL lookup table produced by :py:meth:`calculate`. Rows are indexed by `(sim_str, ensemble, scenario)` tuples; columns are warming level strings.	required

Returns:

Type	Description
`DataArray`	DataArray with an expanded `warming_level` dimension set to `[level]`. The `all_sims` coordinate is relabelled from tuples to `"simulation_scenario"` strings.

Source code in climakitae/explore/warming.py

def find_warming_slice(self, level: str, gwl_times: pd.DataFrame) -> xr.DataArray:
    """Slice catalog data to the window around a single warming level.

    Groups ``catalog_data`` by ``all_sims``, calls :py:func:`get_sliced_data`
    for every simulation, then cleans and relabels the result.

    Parameters
    ----------
    level : str
        Warming level threshold to slice around, e.g. ``"1.5"``.
    gwl_times : pandas.DataFrame
        GWL lookup table produced by :py:meth:`calculate`.  Rows are
        indexed by ``(sim_str, ensemble, scenario)`` tuples; columns are
        warming level strings.

    Returns
    -------
    xarray.DataArray
        DataArray with an expanded ``warming_level`` dimension set to
        ``[level]``.  The ``all_sims`` coordinate is relabelled from
        tuples to ``"simulation_scenario"`` strings.
    """
    warming_data = self.catalog_data.groupby("all_sims").map(
        get_sliced_data,
        level=level,
        years=gwl_times,
        months=self.wl_params.months,
        window=self.wl_params.window,
        anom=self.wl_params.anom,
    )
    warming_data = warming_data.expand_dims({"warming_level": [level]})
    warming_data = warming_data.assign_attrs(
        window=self.wl_params.window, months=self.wl_params.months
    )

    # Cleaning data
    warming_data = clean_warm_data(warming_data)

    # Relabeling `all_sims` dimension
    new_warm_data = warming_data.drop_vars("all_sims")
    new_warm_data["all_sims"] = relabel_axis(warming_data["all_sims"])
    return new_warm_data

`calculate()`

Run the full warming levels pipeline and populate instance attributes.

Retrieves catalog data using wl_params, stacks simulations and scenarios, removes invalid combinations, then iterates over every warming level to produce time-windowed slices and time-mean snapshots.

Results are stored in-place:

catalog_data — raw stacked DataArray
sliced_data — dict mapping level → windowed DataArray
gwl_snapshots — concatenated time-mean DataArray

Notes

Progress is reported via a :py:class:tqdm progress bar, one tick per warming level. If wl_params.load_data is True, each slice is loaded into memory before moving to the next level.

Source code in climakitae/explore/warming.py

def calculate(self):
    """Run the full warming levels pipeline and populate instance attributes.

    Retrieves catalog data using ``wl_params``, stacks simulations and
    scenarios, removes invalid combinations, then iterates over every
    warming level to produce time-windowed slices and time-mean snapshots.

    Results are stored in-place:

    - ``catalog_data`` — raw stacked DataArray
    - ``sliced_data``  — dict mapping level → windowed DataArray
    - ``gwl_snapshots`` — concatenated time-mean DataArray

    Notes
    -----
    Progress is reported via a :py:class:`tqdm` progress bar, one tick per
    warming level.  If ``wl_params.load_data`` is ``True``, each slice is
    loaded into memory before moving to the next level.
    """
    # manually reset to all SSPs, in case it was inadvertently changed by
    # temporarily have ['Dynamical','Statistical'] for downscaling_method
    self.wl_params.scenario_ssp = SSPS

    # Postage data and anomalies
    self.catalog_data = self.wl_params.retrieve()
    self.catalog_data = self.catalog_data.stack(all_sims=["simulation", "scenario"])

    # Dropping invalid simulations that come up from stacking scenarios and simulations together
    self.catalog_data = _drop_invalid_sims(self.catalog_data, self.wl_params)

    if self.wl_params.anom == "Yes":
        self.gwl_times = read_csv_file(GWL_1981_2010_FILE, index_col=[0, 1, 2])
    else:
        self.gwl_times = read_csv_file(GWL_1850_1900_FILE, index_col=[0, 1, 2])
    self.gwl_times = self.gwl_times.dropna(how="all")
    self.catalog_data = clean_list(self.catalog_data, self.gwl_times)

    self.sliced_data = {}
    self.gwl_snapshots = {}
    for level in tqdm(
        self.wl_params.warming_levels, desc="Computing each warming level"
    ):
        warm_slice = self.find_warming_slice(level, self.gwl_times)
        if self.wl_params.load_data:
            warm_slice = load(warm_slice, progress_bar=True)

        # Add GWL snapshots
        self.gwl_snapshots[level] = warm_slice.mean(dim="time", skipna=True)

        # Renaming time dimension for warming slice once "time" is all computed on
        freq_strs = {"monthly": "months", "daily": "days", "hourly": "hours"}
        warm_slice = warm_slice.rename(
            {"time": f"{freq_strs[warm_slice.frequency]}_from_center"}
        )
        self.sliced_data[level] = warm_slice

    self.gwl_snapshots = xr.concat(self.gwl_snapshots.values(), dim="warming_level")

Time Series Visualization

Holds the instance of TimeSeriesParameters that is used for the following purposes: 1) to display a panel that previews various time-series transforms (explore), and 2) to save the transform represented by the current state of that preview into a new variable (output_current).

Parameters:

Name	Type	Description	Default
`data`	`DataArray`	Time series array with no spatial coordinates.	required

Attributes:

Name	Type	Description
`choices`	`TimeSeriesParameters`	Param object containing time series data and analysis parameters.

Source code in climakitae/explore/timeseries.py

def __init__(self, data: xr.DataArray):
    if (
        type(data) != xr.core.dataarray.DataArray
    ):  # Data is NOT in the form of xr.DataArray
        raise ValueError(
            "Please pass an xarray DataArray (e.g. as output by DataParameters.retrieve())."
        )
    else:
        raise_error = False
        error_message = ""
        if "lat" in data.coords:  # Data is NOT area averaged
            raise_error = True
            error_message += "Please pass a timeseries (area average)."
        if not any(
            ["Historical + " in v for v in data.scenario.values]
        ):  # Append historical = False
            if raise_error == True:
                error_message += "\n"
            else:
                raise_error = True
            error_message += (
                "Please append the historical period in your data retrieval."
            )
        if raise_error:  # If any errors
            raise ValueError(error_message)

    self.choices = TimeSeriesParameters(data)

`output_current()`

Output the current attributes of the class to a DataArray object. Allows the data to be easily accessed by the user after modifying the attributes directly in the explore panel, for example.

Returns:

Type	Description
`DataArray`

Source code in climakitae/explore/timeseries.py

def output_current(self) -> xr.DataArray:
    """Output the current attributes of the class to a DataArray object.
    Allows the data to be easily accessed by the user after modifying the attributes directly in the explore panel, for example.

    Returns
    -------
    xr.DataArray

    """
    to_output = self.choices.transform_data()
    attrs_to_add = dict(self.choices.param.values())
    to_output = _update_attrs(to_output, attrs_to_add)
    return to_output

Uncertainty Analysis (CMIP Ensemble Optimization)

A class for holding relevant data options for cmip preprocessing

Attributes:

Name	Type	Description
`variable`	`str`	variable name, cf-compliant (or cmip6 variable name)
`area_subset`	`str`	geographic boundary name (states/counties)
`location`	`str`	geographic area name (name of county/state)
`timescale`	`str`	frequency of data
`area_average`	`bool`	average computed across domain

Methods:

Name	Description
`_cmip_clip`	CMIP6-specific subsetting

Source code in climakitae/explore/uncertainty.py

def __init__(
    self,
    variable: str = "tas",  ## set up for temp uncertainty notebook
    area_subset: str = "states",
    location: str = "California",
    timescale: str = "monthly",
    area_average: bool = True,
) -> None:
    self.variable = variable
    self.area_subset = area_subset
    self.location = location
    self.timescale = timescale
    self.area_average = area_average

Typical Meteorological Year

Encapsulate the code needed to generate Typical Meteorological Year (TMY) files.

Uses WRF hourly data to produce TMYs. User provides the start and end years along with location to generate file.

How to set location: The location can either be provided as latitude and longitude coordinates or as the name of a HadISD station in California. Do not set latitude and longitude if using a HadISD station. If latitude and longitude are set along with a custom value for station_name (NOT a HadISD station), the custom station name will be used in file headers where appropriate.

How to set time period: The time period can either be set with a time approach using start_year and end_year or with a warming level approach using warming_level. If the warming level approach is used, a 30-year period is obtained centered around the given warming level and the start and end years are taken for that warming level.

Parameters:

Name	Type	Description	Default
`start_year`	`str`	Initial year of TMY period (time approach)	`UNSET`
`end_year`	`str`	Final year of TMY period (time approach)	`UNSET`
`warming_level`	`float \| int`	Desired warming level (warming level approach)	`UNSET`
`station_name`	`str`	Long name of desired station	`UNSET`
`latitude`	`float \| int(optional)`	Latitude for TMY data if station_name not set	`UNSET`
`longitude`	`float \| int(optional)`	Longitude for TMY data if station_name not set	`UNSET`
`reanalysis`	`bool`	True to use ERA5 reanalysis instead of models	`False`
`verbose`	`bool`	True to increase verbosity	`True`

Attributes:

Name	Type	Description
`start_year`	`str`	Initial year of TMY period
`end_year`	`str`	Final year of TMY period
`warming_level`	`float \| int`	Warming level value
`lat_range`	`tuple`	Pair of latitudes that bracket `latitude`
`lon_range`	`tuple`	Pair of longitudes that bracket `longitude`
`reanalysis`	`bool`	True to use ERA5 reanalysis instead of models
`simulations`	`list[str]`	List of included simulations
`scenario`	`list[str]`	List of scenarios
`vars_and_units`	`dict[str, str]`	Dictionary of all required variables and units
`verbose`	`bool`	True to increase verbosity
`cdf_climatology`	`Dataset`	CDF climatology data
`cdf_monthly`	`Dataset`	CDF monthly data by model
`weighted_fs_sum`	`Dataset`	Weighted F-S statistic results
`top_months`	`DataFrame`	Table of top months by model
`all_vars`	`Dataset`	All loaded variables for TMY
`tmy_data_to_export`	`dict[Dataframe]`	Dictionary of TMY results by simulation
`_skip_last`	`bool`	Internal flag to track last year for warming level approach

Source code in climakitae/explore/typical_meteorological_year.py

def __init__(
    self,
    start_year: int = UNSET,
    end_year: int = UNSET,
    warming_level: float | int = UNSET,
    station_name: str = UNSET,
    latitude: float | int = UNSET,
    longitude: float | int = UNSET,
    reanalysis: bool = False,
    verbose: bool = True,
):

    # Here we go through a few different ways to get the TMY location
    match latitude, longitude, station_name:
        # UNSET will match to object type
        # Case 1: All variables set
        case float() | int(), float() | int(), str():
            if is_HadISD(station_name):
                raise ValueError(
                    "Do not set `latitude` and `longitude` when using a HadISD station for `station_name`. Change `station_name` value if using custom location."
                )
            else:
                print(
                    f"Initializing TMY object for custom location: {latitude} N, {longitude} W with name '{station_name}'."
                )
                self._set_loc_from_lat_lon(latitude, longitude)
                self.stn_name = station_name
        # Case 2: lat/lon provided, no station_name string
        case float() | int(), float() | int(), object():
            print(
                f"Initializing TMY object for custom location: {latitude} N, {longitude} W."
            )
            self._set_loc_from_lat_lon(latitude, longitude)
        # Case 3: station name provided, lat/lon not numeric
        case object(), object(), str():
            print(f"Initializing TMY object for {station_name}.")
            self._set_loc_from_stn_name(station_name)
        # Last case: something else provided
        case _:
            raise ValueError(
                "No valid station name or latitude and longitude provided."
            )
    # Time variables
    match start_year, end_year, warming_level:
        # User set all options - bad
        case float() | int(), float() | int(), float() | int():
            raise ValueError(
                "Variables `start_year` and `end_year` cannot be paired with `warming_level`. Set either `start_year` and `end_year` OR `warming_level."
            )
        # Some other combo - unset variable will be saved as UNSET
        case _:
            self.start_year = start_year
            self.end_year = end_year
            self.warming_level = warming_level
            if isinstance(self.warming_level, int):
                self.warming_level = float(self.warming_level)
    # Whether to drop the last month as a possible match
    self._skip_last = False
    if self.warming_level is not UNSET:
        # True for warming levels because final hours get lost to UTC conversion
        self._skip_last = True
    # Ranges used in get_data to pull smaller dataset without warnings
    self.lat_range = (self.stn_lat - 0.1, self.stn_lat + 0.1)
    self.lon_range = (self.stn_lon - 0.1, self.stn_lon + 0.1)
    # Determine if using models or ERA5 reanalysis
    if reanalysis:
        # Check that valid time period is requested
        if self.warming_level is not UNSET:
            raise ValueError(
                "Cannot use warming levels with ERA5 reanalysis. Please set start_year and end_year parameters."
            )
        if (self.start_year < 1981) or (self.end_year > 2019):
            raise ValueError(
                f"Valid start and end years for ERA5 reanalysis must be between 1981 and 2019. User provided start year {self.start_year} and end year {self.end_year}"
            )
        self.use_era5 = True
        self.simulations = ["WRF_ERA5_reanalysis"]
    else:  # use models
        self.use_era5 = False
        # These 4 simulations have the solar variables needed
        self.simulations = [
            "WRF_EC-Earth3_r1i1p1f1",
            "WRF_MPI-ESM1-2-HR_r3i1p1f1",
            "WRF_TaiESM1_r1i1p1f1",
            "WRF_MIROC6_r1i1p1f1",
        ]
    # Data only available for these scenarios
    self.scenario = ["Historical Climate", "SSP 3-7.0"]
    # Raw catalog variables to fetch (variable_id → display name)
    # These are fetched directly from the catalog in their native units.
    self._raw_vars = {
        "t2": "Air Temperature at 2m",
        "q2": "Water Vapor Mixing Ratio at 2m",
        "psfc": "Surface Pressure",
        "u10": "u10",
        "v10": "v10",
        "swdnb": "Instantaneous downwelling shortwave flux at bottom",
        "swddni": "Shortwave surface downward direct normal irradiance",
        "swddif": "Shortwave surface downward diffuse irradiance",
        "lwdnb": "Instantaneous downwelling longwave flux at bottom",
    }
    # Full set of TMY variables (including derived) with desired units.
    # Used for display name references throughout the rest of the TMY code.
    self.vars_and_units = {
        "Air Temperature at 2m": "degC",
        "Dew point temperature": "degC",
        "Relative humidity": "[0 to 100]",
        "Instantaneous downwelling shortwave flux at bottom": "W/m2",
        "Shortwave surface downward direct normal irradiance": "W/m2",
        "Shortwave surface downward diffuse irradiance": "W/m2",
        "Instantaneous downwelling longwave flux at bottom": "W/m2",
        "Wind speed at 10m": "m s-1",
        "Wind direction at 10m": "degrees",
        "Surface Pressure": "Pa",
        "Water Vapor Mixing Ratio at 2m": "g kg-1",
    }
    self.verbose = verbose
    # These will get set later in analysis
    self.cdf_climatology = UNSET
    self.cdf_monthly = UNSET
    self.weighted_fs_sum = UNSET
    self.top_months = UNSET
    self.all_vars = UNSET
    self.tmy_data_to_export = UNSET

`load_all_variables()`

Load hourly TMY variables and derive daily statistics for CDF/F-S.

Fetches hourly raw variables via ClimateData for the 8760 profile assembly, then derives ALL daily statistics from the hourly data in local time. This matches the original TMY code's approach and avoids two problems with fetching daily catalog variables directly:

UTC vs local time: Catalog daily variables are pre-aggregated over UTC day boundaries, which differ from local-time days by the station's UTC offset (e.g., 8 hours for California).
Non-determinism: With a Dask distributed client active, lazy reductions (.resample().sum()) can produce slightly different floating-point results on each run due to non-deterministic task ordering. Computing hourly data eagerly to numpy before resampling guarantees deterministic daily statistics.

Source code in climakitae/explore/typical_meteorological_year.py

def load_all_variables(self):
    """Load hourly TMY variables and derive daily statistics for CDF/F-S.

    Fetches hourly raw variables via ClimateData for the 8760 profile
    assembly, then derives ALL daily statistics from the hourly data
    in local time.  This matches the original TMY code's approach and
    avoids two problems with fetching daily catalog variables directly:

    1. **UTC vs local time**: Catalog daily variables are pre-aggregated
       over UTC day boundaries, which differ from local-time days by the
       station's UTC offset (e.g., 8 hours for California).
    2. **Non-determinism**: With a Dask distributed client active, lazy
       reductions (``.resample().sum()``) can produce slightly different
       floating-point results on each run due to non-deterministic task
       ordering.  Computing hourly data eagerly to numpy before
       resampling guarantees deterministic daily statistics.
    """
    print("Loading data from catalog via ClimateData.")

    def _fetch_and_clean(variable_id, table_id="1hr"):
        """Fetch a single variable and clean coordinate cruft."""
        self._vprint(f"  Fetching {variable_id} ({table_id})...")
        da = self._fetch_raw_variable(variable_id, table_id=table_id)
        return da.squeeze().drop_vars(
            [
                "lakemask",
                "landmask",
                "x",
                "y",
                "Lambert_Conformal",
                "centered_year",
            ],
            errors="ignore",
        )

    # --- Hourly variables (needed for 8760 profile assembly) ---
    # First hourly fetch must be synchronous for warming level (sets year range)
    hourly_var_ids = list(self._raw_vars.keys())
    if self.warming_level is not UNSET:
        self._vprint(f"  Getting {hourly_var_ids[0]} (sets year range)...")
        # Fetch raw first to capture centered_year before cleaning drops it
        raw_first = self._fetch_raw_variable(hourly_var_ids[0], table_id="1hr")
        if "centered_year" in raw_first.coords:
            self._sim_centered_years = dict(
                zip(
                    raw_first.simulation.values,
                    raw_first.centered_year.values.ravel(),
                )
            )
        first_var = raw_first.squeeze().drop_vars(
            [
                "lakemask",
                "landmask",
                "x",
                "y",
                "Lambert_Conformal",
                "centered_year",
            ],
            errors="ignore",
        )
        first_var.name = self._raw_vars[hourly_var_ids[0]]
        remaining_hourly = hourly_var_ids[1:]
    else:
        first_var = None
        remaining_hourly = hourly_var_ids

    # --- Daily variables (needed for CDF/F-S analysis) ---
    # Derive daily stats from hourly data in local time.
    # This matches the original code's behavior (hourly → local time →
    # daily resample) and avoids using catalog daily variables which are
    # pre-aggregated over UTC day boundaries (different 24-hour window).
    # We also compute hourly data eagerly here so all downstream
    # operations use deterministic numpy math, avoiding non-deterministic
    # floating-point reduction ordering from the dask distributed scheduler.

    # Fetch remaining hourly variables in parallel
    self._vprint("  Loading hourly variables in parallel...")

    def _fetch_hourly(vid):
        da = _fetch_and_clean(vid, "1hr")
        da.name = self._raw_vars[vid]
        return ("hourly", vid, da)

    with ThreadPoolExecutor(max_workers=4) as executor:
        hourly_futures = [
            executor.submit(_fetch_hourly, vid) for vid in remaining_hourly
        ]
        all_results = [f.result() for f in hourly_futures]

    # Separate hourly results
    hourly_list = [r[2] for r in all_results]
    if first_var is not None:
        hourly_list = [first_var] + hourly_list

    # --- Build hourly dataset for 8760 profile ---
    self._vprint("  Merging raw hourly data.")
    raw_ds = xr.merge(hourly_list)

    # Compute hourly derived variables
    self._vprint("  Computing hourly derived variables.")
    t2_degc = raw_ds["Air Temperature at 2m"] - 273.15
    t2_degc.attrs["units"] = "degC"

    q2_gkg = raw_ds["Water Vapor Mixing Ratio at 2m"] * 1000
    q2_gkg.attrs["units"] = "g kg-1"

    psfc_hpa = raw_ds["Surface Pressure"] / 100.0
    psfc_hpa.attrs["units"] = "hPa"

    rh = compute_relative_humidity(
        pressure=psfc_hpa,
        temperature=t2_degc,
        mixing_ratio=q2_gkg,
    )
    rh.name = "Relative humidity"
    rh.attrs["units"] = "[0 to 100]"

    dew_point_k = compute_dewpointtemp(
        temperature=raw_ds["Air Temperature at 2m"],
        rel_hum=rh,
    )
    dew_point = dew_point_k - 273.15
    dew_point.name = "Dew point temperature"
    dew_point.attrs["units"] = "degC"

    wind_speed = compute_wind_mag(u10=raw_ds["u10"], v10=raw_ds["v10"])
    wind_speed.name = "Wind speed at 10m"

    wind_dir = compute_wind_dir(u10=raw_ds["u10"], v10=raw_ds["v10"])
    wind_dir.name = "Wind direction at 10m"

    t2_out = t2_degc.copy()
    t2_out.name = "Air Temperature at 2m"
    t2_out.attrs["units"] = "degC"

    q2_out = q2_gkg.copy()
    q2_out.name = "Water Vapor Mixing Ratio at 2m"

    derived_list = [t2_out, dew_point, rh, wind_speed, wind_dir, q2_out]
    keep_vars = [
        "Instantaneous downwelling shortwave flux at bottom",
        "Shortwave surface downward direct normal irradiance",
        "Shortwave surface downward diffuse irradiance",
        "Instantaneous downwelling longwave flux at bottom",
        "Surface Pressure",
    ]
    kept_from_raw = [raw_ds[v] for v in keep_vars]
    hourly_ds = xr.merge(derived_list + kept_from_raw)
    # In ERA5 reanalysis case with 1 sim, the "simulation" dim
    # got dropped when datasets were squeezed. Adding it back here.
    if "simulation" not in hourly_ds.dims:
        hourly_ds = hourly_ds.expand_dims("simulation")
    self._hourly_data = hourly_ds

    # --- Build daily dataset for CDF/F-S analysis ---
    # Compute hourly data eagerly so all daily resampling uses
    # deterministic numpy math (not affected by dask scheduler ordering).
    # For a single grid cell this is small (~30yr × 8760hr × 4sims).
    self._vprint("  Computing hourly data for daily resampling...")
    with ProgressBar():
        hourly_computed = hourly_ds.compute()

    self._vprint("  Resampling hourly data to daily statistics...")
    # Temperature: daily max, min, mean (in degC, already converted above)
    daily_tmax = hourly_computed["Air Temperature at 2m"].resample(time="1D").max()
    daily_tmax.name = "Daily max air temperature"
    daily_tmax.attrs["frequency"] = "daily"

    daily_tmin = hourly_computed["Air Temperature at 2m"].resample(time="1D").min()
    daily_tmin.name = "Daily min air temperature"
    daily_tmin.attrs["frequency"] = "daily"

    daily_tmean = (
        hourly_computed["Air Temperature at 2m"].resample(time="1D").mean()
    )
    daily_tmean.name = "Daily mean air temperature"
    daily_tmean.attrs["frequency"] = "daily"

    # Dew point: daily max, min, mean
    daily_dp_max = (
        hourly_computed["Dew point temperature"].resample(time="1D").max()
    )
    daily_dp_max.name = "Daily max dewpoint temperature"
    daily_dp_max.attrs["frequency"] = "daily"

    daily_dp_min = (
        hourly_computed["Dew point temperature"].resample(time="1D").min()
    )
    daily_dp_min.name = "Daily min dewpoint temperature"
    daily_dp_min.attrs["frequency"] = "daily"

    daily_dp_mean = (
        hourly_computed["Dew point temperature"].resample(time="1D").mean()
    )
    daily_dp_mean.name = "Daily mean dewpoint temperature"
    daily_dp_mean.attrs["frequency"] = "daily"

    # Wind speed: daily max, mean
    daily_ws_max = hourly_computed["Wind speed at 10m"].resample(time="1D").max()
    daily_ws_max.name = "Daily max wind speed"
    daily_ws_max.attrs["frequency"] = "daily"

    daily_ws_mean = hourly_computed["Wind speed at 10m"].resample(time="1D").mean()
    daily_ws_mean.name = "Daily mean wind speed"
    daily_ws_mean.attrs["frequency"] = "daily"

    # Radiation: daily sums
    ghi_sum = (
        hourly_computed["Instantaneous downwelling shortwave flux at bottom"]
        .resample(time="1D")
        .sum()
    )
    ghi_sum.name = "Global horizontal irradiance"
    ghi_sum.attrs["frequency"] = "daily"

    dni_sum = (
        hourly_computed["Shortwave surface downward direct normal irradiance"]
        .resample(time="1D")
        .sum()
    )
    dni_sum.name = "Direct normal irradiance"
    dni_sum.attrs["frequency"] = "daily"

    daily_arrays = [
        daily_tmax,
        daily_tmin,
        daily_tmean,
        daily_dp_max,
        daily_dp_min,
        daily_dp_mean,
        daily_ws_max,
        daily_ws_mean,
        ghi_sum,
        dni_sum,
    ]

    self.all_vars = xr.merge(daily_arrays)
    self._vprint("  Daily statistics ready.")

`set_cdf_climatology()`

Calculate the long-term climatology for each index for each month so we can establish the baseline pattern.

Source code in climakitae/explore/typical_meteorological_year.py

def set_cdf_climatology(self):
    """Calculate the long-term climatology for each index for each month so
    we can establish the baseline pattern.
    """
    if self.all_vars is UNSET:
        self.load_all_variables()
    self._vprint("Calculating CDF climatology.")
    self.cdf_climatology = get_cdf(self.all_vars)

`set_cdf_monthly()`

Get CDF for each month and variable.

Source code in climakitae/explore/typical_meteorological_year.py

def set_cdf_monthly(self):
    """Get CDF for each month and variable."""
    if self.all_vars is UNSET:
        self.load_all_variables()
    self._vprint("Calculating monthly CDF.")
    self.cdf_monthly = get_cdf_monthly(self.all_vars)
    # Remove the years for the Pinatubo eruption
    self.cdf_monthly = remove_pinatubo_years(self.cdf_monthly)

`set_weighted_statistic()`

Calculate the weighted F-S statistic.

Source code in climakitae/explore/typical_meteorological_year.py

def set_weighted_statistic(self):
    """Calculate the weighted F-S statistic."""
    if self.cdf_climatology is UNSET:
        self.set_cdf_climatology()
    if self.cdf_monthly is UNSET:
        self.set_cdf_monthly()
    self._vprint("Calculating weighted F-S statistic.")
    self.weighted_fs_sum = compute_weighted_fs_sum(
        self.cdf_climatology, self.cdf_monthly
    )

`set_top_months()`

Calculate top months dataframe.

Source code in climakitae/explore/typical_meteorological_year.py

def set_top_months(self):
    """Calculate top months dataframe."""
    # Pass the weighted F-S sum data for simplicity
    if self.weighted_fs_sum is UNSET:
        self.set_weighted_statistic()
    self._vprint("Finding top months (lowest F-S statistic)")
    self.top_months = get_top_months(
        self.weighted_fs_sum, skip_last=self._skip_last
    )

`show_tmy_data_to_export(simulation)`

Show line plots of TMY data for single model.

Parameters:

Name	Type	Description	Default
`simulation`	`str`	Simulation to display.	required

Source code in climakitae/explore/typical_meteorological_year.py

def show_tmy_data_to_export(self, simulation: str):
    """Show line plots of TMY data for single model.

    Parameters
    ----------
    simulation: str
        Simulation to display.

    """
    if self.tmy_data_to_export is UNSET:
        print("No TMY data generated.")
        print("Please run TMY.generate_tmy() to create TMY data for viewing.")
        return

    # WVMR not in final TMY dataframe
    fig_y = list(self.vars_and_units.keys())
    fig_y.remove("Water Vapor Mixing Ratio at 2m")

    self.tmy_data_to_export[simulation].plot(
        x="time",
        y=fig_y,
        title=f"Typical Meteorological Year ({simulation})",
        subplots=True,
        figsize=(10, 8),
        legend=True,
    )

`run_tmy_analysis()`

Generate typical meteorological year data.

Output will be a list of dataframes per simulation. Print statements throughout the function indicate progress.

Notes

Results are saved to the class variable tmy_data_to_export.

Source code in climakitae/explore/typical_meteorological_year.py

def run_tmy_analysis(self):
    """Generate typical meteorological year data.

    Output will be a list of dataframes per simulation.
    Print statements throughout the function indicate progress.

    Notes
    -----
    Results are saved to the class variable `tmy_data_to_export`.
    """
    print("Assembling TMY data to export.")

    self._vprint("  STEP 1: Computing hourly data")

    # Use cached hourly data instead of re-downloading
    if not hasattr(self, "_hourly_data") or self._hourly_data is None:
        # Fallback: load from catalog if run_tmy_analysis called standalone
        self.load_all_variables()
    with ProgressBar():
        all_vars_ds = self._hourly_data.compute()

    # Construct TMY
    self._vprint(
        "\n  STEP 2: Calculating Typical Meteorological Year per model simulation\n  Progress bar shows code looping through each month in the year.\n"
    )

    tmy_data_to_export = self._make_8760_tables(
        all_vars_ds, self.top_months
    )  # Return dict of TMY by simulation

    self._vprint("  Smoothing data at transitions between months.")
    self._vprint("  Dropping water vapor mixing ratio.")
    # Smooth transition hours
    for sim in tmy_data_to_export:
        tmy_data_to_export[sim] = tmy_data_to_export[sim].reset_index()
        tmy_data_to_export[sim] = self._smooth_month_transition_hours(
            tmy_data_to_export[sim]
        )

        # Mixing ratio was only needed for smoothing relative humidity,
        # so it can be dropped now.
        tmy_data_to_export[sim] = tmy_data_to_export[sim].drop(
            columns="Water Vapor Mixing Ratio at 2m"
        )

        # Add metadata columns needed by EPW header writer.
        # The new-core pipeline squeezes these dimensions, so they must be
        # re-attached before export.
        if self.warming_level is not UNSET:
            tmy_data_to_export[sim]["warming_level"] = self.warming_level
        else:
            if self.use_era5:
                tmy_data_to_export[sim]["scenario"] = "historical_reanalysis"
            else:
                tmy_data_to_export[sim]["scenario"] = "historical+ssp370"

        # Set floating point precision for specific vars
        float_cols = [
            "Air Temperature at 2m",
            "Dew point temperature",
            "Relative humidity",
            "Surface Pressure",
            "Instantaneous downwelling longwave flux at bottom",
            "Instantaneous downwelling shortwave flux at bottom",
            "Shortwave surface downward direct normal irradiance",
            "Shortwave surface downward diffuse irradiance",
            "Wind direction at 10m",
        ]

        # precision 1 point for all vars
        tmy_data_to_export[sim][float_cols] = tmy_data_to_export[sim][
            float_cols
        ].round(1)
        # precsion 2 point for wind speed only
        tmy_data_to_export[sim]["Wind speed at 10m"] = tmy_data_to_export[sim][
            "Wind speed at 10m"
        ].round(2)

    self.tmy_data_to_export = tmy_data_to_export
    self._vprint("TMY analysis complete.")

`export_tmy_data(extension='epw')`

Write TMY data to EPW file.

Parameters:

Name	Type	Description	Default
`extension`	`str`	Desired file extension ('tmy','epw', or 'csv')	`'epw'`

Source code in climakitae/explore/typical_meteorological_year.py

def export_tmy_data(self, extension: str = "epw"):
    """Write TMY data to EPW file.

    Parameters
    ----------
    extension: str
        Desired file extension ('tmy','epw', or 'csv')

    """
    print("Exporting TMY to file.")
    for sim, _ in self.tmy_data_to_export.items():
        # Get right year range
        if self.warming_level is UNSET:
            years = (self.start_year, self.end_year)
            clean_sim = sim
        else:
            centered_year = int(self._sim_centered_years[sim])
            year1 = centered_year - 15
            year2 = centered_year + 14
            years = (year1, year2)
            # Append warming level descriptor to simulation name.
            # Legacy sim names no longer contain "_historical+ssp370"
            # after the new-core migration, so we append directly.
            wl_label = match_str_to_wl(self.warming_level)
            clean_sim = f"{sim}_{wl_label}"
        # Attach centered_year so CSV export can include it
        if self.warming_level is not UNSET:
            self.tmy_data_to_export[sim]["centered_year"] = centered_year
        clean_stn_name = (
            self.stn_name.replace(" ", "_").replace("(", "").replace(")", "")
        )
        filename = f"TMY_{clean_stn_name}_{clean_sim}".lower()
        write_tmy_file(
            filename,
            self.tmy_data_to_export[sim],
            years,
            self.stn_name,
            self.stn_code,
            self.stn_lat,
            self.stn_lon,
            self.stn_state,
            file_ext=extension,
        )

`get_candidate_months()`

Run CDF functions to get top candidates.

This function can be used to view the candidate months without running the entire TMY workflow.

Source code in climakitae/explore/typical_meteorological_year.py

def get_candidate_months(self):
    """Run CDF functions to get top candidates.

    This function can be used to view the candidate months
    without running the entire TMY workflow.
    """
    self._vprint("Getting top months for TMY.")
    self.set_cdf_climatology()
    self.set_cdf_monthly()
    self.set_weighted_statistic()
    self.set_top_months()

`generate_tmy()`

Run the whole TMY workflow.

Source code in climakitae/explore/typical_meteorological_year.py

def generate_tmy(self):
    """Run the whole TMY workflow."""
    # This runs the whole workflow at once
    print("Running TMY workflow.")
    self.load_all_variables()
    self.get_candidate_months()
    self.run_tmy_analysis()
    self.export_tmy_data()

Extreme Meteorological Year

Functions for Shock and Persistence Extreme Meteorological Year creation.

`shock_XMY(extreme=UNSET, start_year=UNSET, end_year=UNSET, warming_level=UNSET, station_name=UNSET, latitude=UNSET, longitude=UNSET, reanalysis=False, verbose=True)`

Encapsulate the code needed to generate Typical Meteorological Year (shock XMY) files.

Uses WRF hourly data to produce shock XMYs. User provides the start and end years along with location to generate file.

How to set location: The location can either be provided as latitude and longitude coordinates or as the name of a HadISD station in California. Do not set latitude and longitude if using a HadISD station. If latitude and longitude are set along with a custom value for station_name (NOT a HadISD station), the custom station name will be used in file headers where appropriate.

How to set time period: The time period can either be set with a time approach using start_year and end_year or with a warming level approach using warming_level. If the warming level approach is used, a 30-year period is obtained centered around the given warming level and the start and end years are taken for that warming level.

How to set extreme: The extreme can either be set as 'hot' or 'cold'. The 'hot' extreme corresponds to extreme events like heat waves, while the 'cold' extreme corresponds to extreme evens like cold snaps.

Parameters:

Name	Type	Description	Default
`extreme`	`str`	Type of shock extreme	`UNSET`
`start_year`	`str`	Initial year of XMY period (time approach)	`UNSET`
`end_year`	`str`	Final year of XMY period (time approach)	`UNSET`
`warming_level`	`float \| int`	Desired warming level (warming level approach)	`UNSET`
`station_name`	`str`	Long name of desired station	`UNSET`
`latitude`	`float \| int(optional)`	Latitude for shock XMY data if station_name not set	`UNSET`
`longitude`	`float \| int(optional)`	Longitude for shock XMY data if station_name not set	`UNSET`
`reanalysis`	`bool`	True to use ERA5 reanalysis instead of models	`False`
`verbose`	`bool`	True to increase verbosity	`True`

Attributes:

Name	Type	Description
`extreme`	`str`	extreme type, either "hot" or "cold"
`start_year`	`str`	Initial year of shock XMY period
`end_year`	`str`	Final year of shock XMY period
`warming_level`	`float \| int`	Warming level value
`lat_range`	`tuple`	Pair of latitudes that bracket `latitude`
`lon_range`	`tuple`	Pair of longitudes that bracket `longitude`
`reanalysis`	`bool`	True to use ERA5 reanalysis instead of models
`simulations`	`list[str]`	List of included simulations
`scenario`	`list[str]`	List of scenarios
`vars_and_units`	`dict[str, str]`	Dictionary of all required variables and units
`verbose`	`bool`	True to increase verbosity
`cdf_climatology`	`Dataset`	CDF climatology data
`cdf_monthly`	`Dataset`	CDF monthly data by model
`top_months`	`DataFrame`	Table of top months by model
`all_vars`	`Dataset`	All loaded variables for shock XMY
`air_temp_vars`	`Dataset`	Air temperature variables, for use in finding candidate months
`xmy_data_to_export`	`dict[Dataframe]`	Dictionary of shock XMY results by simulation
`_skip_last`	`bool`	Internal flag to track last year for warming level approach

Source code in climakitae/explore/extreme_meteorological_year.py

def __init__(
    self,
    extreme: str = UNSET,
    start_year: int = UNSET,
    end_year: int = UNSET,
    warming_level: float | int = UNSET,
    station_name: str = UNSET,
    latitude: float | int = UNSET,
    longitude: float | int = UNSET,
    reanalysis: bool = False,
    verbose: bool = True,
):

    # Here we go through a few different ways to get the shock XMY location
    match latitude, longitude, station_name:
        # UNSET will match to object type
        # Case 1: All variables set
        case float() | int(), float() | int(), str():
            if is_HadISD(station_name):
                raise ValueError(
                    "Do not set `latitude` and `longitude` when using a HadISD station for `station_name`. Change `station_name` value if using custom location."
                )
            else:
                print(
                    f"Initializing shock XMY object for custom location: {latitude} N, {longitude} W with name '{station_name}'."
                )
                self._set_loc_from_lat_lon(latitude, longitude)
                self.stn_name = station_name
        # Case 2: lat/lon provided, no station_name string
        case float() | int(), float() | int(), object():
            print(
                f"Initializing shock XMY object for custom location: {latitude} N, {longitude} W."
            )
            self._set_loc_from_lat_lon(latitude, longitude)
        # Case 3: station name provided, lat/lon not numeric
        case object(), object(), str():
            print(f"Initializing shock XMY object for {station_name}.")
            self._set_loc_from_stn_name(station_name)
        # Last case: something else provided
        case _:
            raise ValueError(
                "No valid station name or latitude and longitude provided."
            )
    # Time variables
    match start_year, end_year, warming_level:
        # User set all options - bad
        case float() | int(), float() | int(), float() | int():
            raise ValueError(
                "Variables `start_year` and `end_year` cannot be paired with `warming_level`. Set either `start_year` and `end_year` OR `warming_level."
            )
        # Some other combo - unset variable will be saved as UNSET
        case _:
            self.start_year = start_year
            self.end_year = end_year
            self.warming_level = warming_level
            if isinstance(self.warming_level, int):
                self.warming_level = float(self.warming_level)
    # extreme type
    if extreme in ["hot", "cold"]:
        self.extreme = extreme
        print(f"Generating {extreme} shock XMY.")
    else:
        raise TypeError("Variable `extreme` must be either 'hot' or 'cold'.")
    # Whether to drop the last month as a possible match
    self._skip_last = False
    if self.warming_level is not UNSET:
        # True for warming levels because final hours get lost to UTC conversion
        self._skip_last = True
    # Ranges used in get_data to pull smaller dataset without warnings
    self.lat_range = (self.stn_lat - 0.1, self.stn_lat + 0.1)
    self.lon_range = (self.stn_lon - 0.1, self.stn_lon + 0.1)
    # Determine if using models or ERA5 reanalysis
    if reanalysis:
        # Check that valid time period is requested
        if self.warming_level is not UNSET:
            raise ValueError(
                "Cannot use warming levels with ERA5 reanalysis. Please set start_year and end_year parameters."
            )
        if (self.start_year < 1981) or (self.end_year > 2019):
            raise ValueError(
                f"Valid start and end years for ERA5 reanalysis must be between 1981 and 2019. User provided start year {self.start_year} and end year {self.end_year}"
            )
        self.use_era5 = True
        self.simulations = ["WRF_ERA5_reanalysis"]
    else:  # use models
        self.use_era5 = False
        # These 4 simulations have the solar variables needed
        self.simulations = [
            "WRF_EC-Earth3_r1i1p1f1",
            "WRF_MPI-ESM1-2-HR_r3i1p1f1",
            "WRF_TaiESM1_r1i1p1f1",
            "WRF_MIROC6_r1i1p1f1",
        ]

    # Data only available for these scenarios
    self.scenario = ["Historical Climate", "SSP 3-7.0"]
    # Raw catalog variables to fetch (variable_id → display name)
    # These are fetched directly from the catalog in their native units.
    self._raw_vars = {
        "t2": "Air Temperature at 2m",
        "q2": "Water Vapor Mixing Ratio at 2m",
        "psfc": "Surface Pressure",
        "u10": "u10",
        "v10": "v10",
        "swdnb": "Instantaneous downwelling shortwave flux at bottom",
        "swddni": "Shortwave surface downward direct normal irradiance",
        "swddif": "Shortwave surface downward diffuse irradiance",
        "lwdnb": "Instantaneous downwelling longwave flux at bottom",
    }
    # Full set of shock XMY variables (including derived) with desired units.
    # Used for display name references throughout the rest of the shock XMY code.
    self.vars_and_units = {
        "Air Temperature at 2m": "degC",
        "Dew point temperature": "degC",
        "Relative humidity": "[0 to 100]",
        "Instantaneous downwelling shortwave flux at bottom": "W/m2",
        "Shortwave surface downward direct normal irradiance": "W/m2",
        "Shortwave surface downward diffuse irradiance": "W/m2",
        "Instantaneous downwelling longwave flux at bottom": "W/m2",
        "Wind speed at 10m": "m s-1",
        "Wind direction at 10m": "degrees",
        "Surface Pressure": "Pa",
        "Water Vapor Mixing Ratio at 2m": "g kg-1",
    }
    self.verbose = verbose
    # These will get set later in analysis
    self.cdf_climatology = UNSET
    self.cdf_monthly = UNSET
    self.top_months = UNSET
    self.all_vars = UNSET
    self.air_temp_vars = UNSET
    self.xmy_data_to_export = UNSET

`load_all_variables()`

Load hourly shock XMY variables and derive daily statistics.

Fetches hourly raw variables via ClimateData for the 8760 profile assembly, then derives ALL daily statistics from the hourly data in local time. This approach avoids two problems with fetching daily catalog variables directly:

UTC vs local time: Catalog daily variables are pre-aggregated over UTC day boundaries, which differ from local-time days by the station's UTC offset (e.g., 8 hours for California).
Non-determinism: With a Dask distributed client active, lazy reductions (.resample().sum()) can produce slightly different floating-point results on each run due to non-deterministic task ordering. Computing hourly data eagerly to numpy before resampling guarantees deterministic daily statistics.

Source code in climakitae/explore/extreme_meteorological_year.py

def load_all_variables(self):
    """Load hourly shock XMY variables and derive daily statistics.

    Fetches hourly raw variables via ClimateData for the 8760 profile
    assembly, then derives ALL daily statistics from the hourly data
    in local time.  This approach avoids two problems with fetching
    daily catalog variables directly:

    1. **UTC vs local time**: Catalog daily variables are pre-aggregated
       over UTC day boundaries, which differ from local-time days by the
       station's UTC offset (e.g., 8 hours for California).
    2. **Non-determinism**: With a Dask distributed client active, lazy
       reductions (``.resample().sum()``) can produce slightly different
       floating-point results on each run due to non-deterministic task
       ordering.  Computing hourly data eagerly to numpy before
       resampling guarantees deterministic daily statistics.
    """
    print("Loading data from catalog via ClimateData.")

    def _fetch_and_clean(variable_id, table_id="1hr"):
        """Fetch a single variable and clean coordinate cruft."""
        self._vprint(f"  Fetching {variable_id} ({table_id})...")
        da = self._fetch_raw_variable(variable_id, table_id=table_id)
        return da.squeeze().drop_vars(
            [
                "lakemask",
                "landmask",
                "x",
                "y",
                "Lambert_Conformal",
                "centered_year",
            ],
            errors="ignore",
        )

    # --- Hourly variables (needed for 8760 profile assembly) ---
    # First hourly fetch must be synchronous for warming level (sets year range)
    hourly_var_ids = list(self._raw_vars.keys())
    if self.warming_level is not UNSET:
        self._vprint(f"  Getting {hourly_var_ids[0]} (sets year range)...")
        # Fetch raw first to capture centered_year before cleaning drops it
        raw_first = self._fetch_raw_variable(hourly_var_ids[0], table_id="1hr")
        if "centered_year" in raw_first.coords:
            self._sim_centered_years = dict(
                zip(
                    raw_first.simulation.values,
                    raw_first.centered_year.values.ravel(),
                )
            )
        first_var = raw_first.squeeze().drop_vars(
            [
                "lakemask",
                "landmask",
                "x",
                "y",
                "Lambert_Conformal",
                "centered_year",
            ],
            errors="ignore",
        )
        first_var.name = self._raw_vars[hourly_var_ids[0]]
        remaining_hourly = hourly_var_ids[1:]
    else:
        first_var = None
        remaining_hourly = hourly_var_ids

    # --- Daily variables resampling ---
    # Derive daily stats from hourly data in local time.
    # This matches the original code's behavior (hourly → local time →
    # daily resample) and avoids using catalog daily variables which are
    # pre-aggregated over UTC day boundaries (different 24-hour window).
    # We also compute hourly data eagerly here so all downstream
    # operations use deterministic numpy math, avoiding non-deterministic
    # floating-point reduction ordering from the dask distributed scheduler.

    # Fetch remaining hourly variables in parallel
    self._vprint("  Loading hourly variables in parallel...")

    def _fetch_hourly(vid):
        da = _fetch_and_clean(vid, "1hr")
        da.name = self._raw_vars[vid]
        return ("hourly", vid, da)

    with ThreadPoolExecutor(max_workers=4) as executor:
        hourly_futures = [
            executor.submit(_fetch_hourly, vid) for vid in remaining_hourly
        ]
        all_results = [f.result() for f in hourly_futures]

    # Separate hourly results
    hourly_list = [r[2] for r in all_results]
    if first_var is not None:
        hourly_list = [first_var] + hourly_list

    # --- Build hourly dataset for 8760 profile ---
    self._vprint("  Merging raw hourly data.")
    raw_ds = xr.merge(hourly_list)

    # Compute hourly derived variables
    self._vprint("  Computing hourly derived variables.")
    t2_degc = raw_ds["Air Temperature at 2m"] - 273.15
    t2_degc.attrs["units"] = "degC"

    q2_gkg = raw_ds["Water Vapor Mixing Ratio at 2m"] * 1000
    q2_gkg.attrs["units"] = "g kg-1"

    psfc_hpa = raw_ds["Surface Pressure"] / 100.0
    psfc_hpa.attrs["units"] = "hPa"

    rh = compute_relative_humidity(
        pressure=psfc_hpa,
        temperature=t2_degc,
        mixing_ratio=q2_gkg,
    )
    rh.name = "Relative humidity"
    rh.attrs["units"] = "[0 to 100]"

    dew_point_k = compute_dewpointtemp(
        temperature=raw_ds["Air Temperature at 2m"],
        rel_hum=rh,
    )
    dew_point = dew_point_k - 273.15
    dew_point.name = "Dew point temperature"
    dew_point.attrs["units"] = "degC"

    wind_speed = compute_wind_mag(u10=raw_ds["u10"], v10=raw_ds["v10"])
    wind_speed.name = "Wind speed at 10m"

    wind_dir = compute_wind_dir(u10=raw_ds["u10"], v10=raw_ds["v10"])
    wind_dir.name = "Wind direction at 10m"

    t2_out = t2_degc.copy()
    t2_out.name = "Air Temperature at 2m"
    t2_out.attrs["units"] = "degC"

    q2_out = q2_gkg.copy()
    q2_out.name = "Water Vapor Mixing Ratio at 2m"

    derived_list = [t2_out, dew_point, rh, wind_speed, wind_dir, q2_out]
    keep_vars = [
        "Instantaneous downwelling shortwave flux at bottom",
        "Shortwave surface downward direct normal irradiance",
        "Shortwave surface downward diffuse irradiance",
        "Instantaneous downwelling longwave flux at bottom",
        "Surface Pressure",
    ]
    kept_from_raw = [raw_ds[v] for v in keep_vars]
    hourly_ds = xr.merge(derived_list + kept_from_raw)
    # In ERA5 reanalysis case with 1 sim, the "simulation" dim
    # got dropped when datasets were squeezed. Adding it back here.
    if "simulation" not in hourly_ds.dims:
        hourly_ds = hourly_ds.expand_dims("simulation")
    self._hourly_data = hourly_ds

    # --- Build daily dataset ---

    # Compute hourly data eagerly so all daily resampling uses
    # deterministic numpy math (not affected by dask scheduler ordering).
    # For a single grid cell this is small (~30yr × 8760hr × 4sims).
    self._vprint("  Computing hourly data for daily resampling...")
    with ProgressBar():
        hourly_computed = hourly_ds.compute()

    self._vprint("  Resampling hourly data to daily statistics...")
    # Temperature: daily max, min, mean (in degC, already converted above)
    daily_tmax = hourly_computed["Air Temperature at 2m"].resample(time="1D").max()
    daily_tmax.name = "Daily max air temperature"
    daily_tmax.attrs["frequency"] = "daily"

    daily_tmin = hourly_computed["Air Temperature at 2m"].resample(time="1D").min()
    daily_tmin.name = "Daily min air temperature"
    daily_tmin.attrs["frequency"] = "daily"

    daily_tmean = (
        hourly_computed["Air Temperature at 2m"].resample(time="1D").mean()
    )
    daily_tmean.name = "Daily mean air temperature"
    daily_tmean.attrs["frequency"] = "daily"

    # Dew point: daily max, min, mean
    daily_dp_max = (
        hourly_computed["Dew point temperature"].resample(time="1D").max()
    )
    daily_dp_max.name = "Daily max dewpoint temperature"
    daily_dp_max.attrs["frequency"] = "daily"

    daily_dp_min = (
        hourly_computed["Dew point temperature"].resample(time="1D").min()
    )
    daily_dp_min.name = "Daily min dewpoint temperature"
    daily_dp_min.attrs["frequency"] = "daily"

    daily_dp_mean = (
        hourly_computed["Dew point temperature"].resample(time="1D").mean()
    )
    daily_dp_mean.name = "Daily mean dewpoint temperature"
    daily_dp_mean.attrs["frequency"] = "daily"

    # Wind speed: daily max, mean
    daily_ws_max = hourly_computed["Wind speed at 10m"].resample(time="1D").max()
    daily_ws_max.name = "Daily max wind speed"
    daily_ws_max.attrs["frequency"] = "daily"

    daily_ws_mean = hourly_computed["Wind speed at 10m"].resample(time="1D").mean()
    daily_ws_mean.name = "Daily mean wind speed"
    daily_ws_mean.attrs["frequency"] = "daily"

    # Radiation: daily sums
    ghi_sum = (
        hourly_computed["Instantaneous downwelling shortwave flux at bottom"]
        .resample(time="1D")
        .sum()
    )
    ghi_sum.name = "Global horizontal irradiance"
    ghi_sum.attrs["frequency"] = "daily"

    dni_sum = (
        hourly_computed["Shortwave surface downward direct normal irradiance"]
        .resample(time="1D")
        .sum()
    )
    dni_sum.name = "Direct normal irradiance"
    dni_sum.attrs["frequency"] = "daily"

    daily_arrays = [
        daily_tmax,
        daily_tmin,
        daily_tmean,
        daily_dp_max,
        daily_dp_min,
        daily_dp_mean,
        daily_ws_max,
        daily_ws_mean,
        ghi_sum,
        dni_sum,
    ]

    air_temp_var_arrays = [daily_tmax, daily_tmin]

    self.all_vars = xr.merge(daily_arrays)
    self._vprint("  Daily statistics ready.")
    # subet all loaded variables to only air temperature, which is the only variable used to determine the candidate month
    self.air_temp_vars = xr.merge(air_temp_var_arrays)
    self._vprint(
        "   Air temperature variables stored for use in finding candidate months."
    )

`set_cdf_climatology()`

Calculate the long-term climatology for each index for each month so we can establish the baseline pattern.

Source code in climakitae/explore/extreme_meteorological_year.py

def set_cdf_climatology(self):
    """Calculate the long-term climatology for each index for each month so
    we can establish the baseline pattern.
    """
    if self.air_temp_vars is UNSET:
        self.load_all_variables()
    self._vprint("Calculating CDF climatology using only air temp variables.")
    self.cdf_climatology = get_cdf(self.air_temp_vars)

`set_cdf_monthly()`

Get CDF for each month and variable.

Source code in climakitae/explore/extreme_meteorological_year.py

def set_cdf_monthly(self):
    """Get CDF for each month and variable."""
    if self.air_temp_vars is UNSET:
        self.load_all_variables()
    self._vprint("Calculating monthly CDF using only air temp variables.")
    self.cdf_monthly = get_cdf_monthly(self.air_temp_vars)
    # Remove the years for the Pinatubo eruption
    self.cdf_monthly = remove_pinatubo_years(self.cdf_monthly)

`set_top_months()`

Calculate top months dataframe.

Source code in climakitae/explore/extreme_meteorological_year.py

def set_top_months(self):
    """Calculate top months dataframe."""

    if self.cdf_climatology is UNSET:
        self.set_cdf_climatology()

    if self.cdf_monthly is UNSET:
        self.set_cdf_monthly()

    self._vprint(
        "Finding top months (greatest deviation from climatological CDF median)."
    )
    self.top_months = generate_candidate_months(
        self.cdf_monthly, self.cdf_climatology, self.extreme, self._skip_last
    )

`show_xmy_data_to_export(simulation)`

Show line plots of shock XMY data for single model.

Parameters:

Name	Type	Description	Default
`simulation`	`str`	Simulation to display.	required

Source code in climakitae/explore/extreme_meteorological_year.py

def show_xmy_data_to_export(self, simulation: str):
    """Show line plots of shock XMY data for single model.

    Parameters
    ----------
    simulation: str
        Simulation to display.

    """
    if self.xmy_data_to_export is UNSET:
        print("No XMY data generated.")
        print("Please run xmy.generate_xmy() to create XMY data for viewing.")
        return

    # WVMR not in final XMY dataframe
    fig_y = list(self.vars_and_units.keys())
    fig_y.remove("Water Vapor Mixing Ratio at 2m")

    self.xmy_data_to_export[simulation].plot(
        x="time",
        y=fig_y,
        title=f"Shock Extreme Meteorological Year ({simulation})",
        subplots=True,
        figsize=(10, 8),
        legend=True,
    )

`run_xmy_analysis()`

Generate shock extreme meteorological year data.

Output will be a list of dataframes per simulation. Print statements throughout the function indicate progress.

Notes

Results are saved to the class variable xmy_data_to_export.

Source code in climakitae/explore/extreme_meteorological_year.py

def run_xmy_analysis(self):
    """Generate shock extreme meteorological year data.

    Output will be a list of dataframes per simulation.
    Print statements throughout the function indicate progress.

    Notes
    -----
    Results are saved to the class variable `xmy_data_to_export`.
    """
    print("Assembling XMY data to export.")

    self._vprint("  STEP 1: Computing hourly data")

    # Use cached hourly data instead of re-downloading
    if not hasattr(self, "_hourly_data") or self._hourly_data is None:
        # Fallback: load from catalog if run_xmy_analysis called standalone
        self.load_all_variables()
    with ProgressBar():
        all_vars_ds = self._hourly_data.compute()

    # Construct XMY
    self._vprint(
        f"\n  STEP 2: Calculating {self.extreme} Shock Extreme Meteorological Year per model simulation\n  Progress bar shows code looping through each month in the year.\n"
    )

    xmy_data_to_export = self._make_8760_tables(
        all_vars_ds, self.top_months
    )  # Return dict of shock XMY by simulation

    self._vprint("  Smoothing data at transitions between months.")
    self._vprint("  Dropping water vapor mixing ratio.")
    # Smooth transition hours
    for sim in xmy_data_to_export:
        xmy_data_to_export[sim] = xmy_data_to_export[sim].reset_index()
        xmy_data_to_export[sim] = self._smooth_month_transition_hours(
            xmy_data_to_export[sim]
        )

        # Mixing ratio was only needed for smoothing relative humidity,
        # so it can be dropped now.
        xmy_data_to_export[sim] = xmy_data_to_export[sim].drop(
            columns="Water Vapor Mixing Ratio at 2m"
        )

        # Add metadata columns needed by EPW header writer.
        # The new-core pipeline squeezes these dimensions, so they must be
        # re-attached before export.
        if self.warming_level is not UNSET:
            xmy_data_to_export[sim]["warming_level"] = self.warming_level
        else:
            if self.use_era5:
                xmy_data_to_export[sim]["scenario"] = "historical_reanalysis"
            else:
                xmy_data_to_export[sim]["scenario"] = "historical+ssp370"

        # Set floating point precision for specific vars
        float_cols = [
            "Air Temperature at 2m",
            "Dew point temperature",
            "Relative humidity",
            "Surface Pressure",
            "Instantaneous downwelling longwave flux at bottom",
            "Instantaneous downwelling shortwave flux at bottom",
            "Shortwave surface downward direct normal irradiance",
            "Shortwave surface downward diffuse irradiance",
            "Wind direction at 10m",
        ]

        # precision 1 point for all vars
        xmy_data_to_export[sim][float_cols] = xmy_data_to_export[sim][
            float_cols
        ].round(1)
        # precsion 2 point for wind speed only
        xmy_data_to_export[sim]["Wind speed at 10m"] = xmy_data_to_export[sim][
            "Wind speed at 10m"
        ].round(2)

    self.xmy_data_to_export = xmy_data_to_export
    self._vprint("shock XMY analysis complete.")

`export_xmy_data(extension='epw')`

Write XMY data to EPW file.

Parameters:

Name	Type	Description	Default
`extension`	`str`	Desired file extension ('tmy','epw', or 'csv')	`'epw'`

Source code in climakitae/explore/extreme_meteorological_year.py

def export_xmy_data(self, extension: str = "epw"):
    """Write XMY data to EPW file.

    Parameters
    ----------
    extension: str
        Desired file extension ('tmy','epw', or 'csv')

    """
    print("Exporting shock XMY to file.")
    for sim, _ in self.xmy_data_to_export.items():
        # Get right year range
        if self.warming_level is UNSET:
            years = (self.start_year, self.end_year)
            clean_sim = sim
        else:
            centered_year = int(self._sim_centered_years[sim])
            year1 = centered_year - 15
            year2 = centered_year + 14
            years = (year1, year2)
            # Append warming level descriptor to simulation name.
            # Legacy sim names no longer contain "_historical+ssp370"
            # after the new-core migration, so we append directly.
            wl_label = match_str_to_wl(self.warming_level)
            clean_sim = f"{sim}_{wl_label}"
        # Attach centered_year so CSV export can include it
        if self.warming_level is not UNSET:
            self.xmy_data_to_export[sim]["centered_year"] = centered_year
        clean_stn_name = (
            self.stn_name.replace(" ", "_").replace("(", "").replace(")", "")
        )
        filename = f"{self.extreme}_shock_xmy_{clean_stn_name}_{clean_sim}".lower()
        write_tmy_file(
            filename,
            self.xmy_data_to_export[sim],
            years,
            self.stn_name,
            self.stn_code,
            self.stn_lat,
            self.stn_lon,
            self.stn_state,
            file_ext=extension,
            xmy_header={"extreme_type": self.extreme, "xmy_type": "SHOCK"},
        )

`get_candidate_months()`

Run CDF functions to get top candidates.

This function can be used to view the candidate months without running the entire shock XMY workflow.

Source code in climakitae/explore/extreme_meteorological_year.py

def get_candidate_months(self):
    """Run CDF functions to get top candidates.

    This function can be used to view the candidate months
    without running the entire shock XMY workflow.
    """
    self._vprint(f"Getting top months for {self.extreme} shock XMY.")
    self.set_cdf_climatology()
    self.set_cdf_monthly()
    self.set_top_months()

`generate_xmy()`

Run the whole XMY workflow.

Source code in climakitae/explore/extreme_meteorological_year.py

def generate_xmy(self):
    """Run the whole XMY workflow."""
    # This runs the whole workflow at once
    print("Running shock XMY workflow.")
    self.load_all_variables()
    self.get_candidate_months()
    self.run_xmy_analysis()
    self.export_xmy_data()

`persistence_XMY(q=UNSET, start_year=UNSET, end_year=UNSET, warming_level=UNSET, station_name=UNSET, latitude=UNSET, longitude=UNSET, reanalysis=False, verbose=True)`

Encapsulate the code needed to generate Persistence Extreme Meteorological Year (persistence XMY) files.

Uses WRF hourly data to produce persistence XMYs. User provides the start and end years along with location to generate file.

How to set location: The location can either be provided as latitude and longitude coordinates or as the name of a HadISD station in California. Do not set latitude and longitude if using a HadISD station. If latitude and longitude are set along with a custom value for station_name (NOT a HadISD station), the custom station name will be used in file headers where appropriate.

How to set time period: The time period can either be set with a time approach using start_year and end_year or with a warming level approach using warming_level. If the warming level approach is used, a 30-year period is obtained centered around the given warming level and the start and end years are taken for that warming level.

How to set extreme: The extreme is defined by the quantile value, q, and can be a value betwee 0 and 1. For example, q=0.10 will produce a cold extreme, while a q=0.90 will produce a hot extreme.

Parameters:

Name	Type	Description	Default
`q`	`float`	Extreme quantile	`UNSET`
`start_year`	`str`	Initial year of XMY period (time approach)	`UNSET`
`end_year`	`str`	Final year of XMY period (time approach)	`UNSET`
`warming_level`	`float \| int`	Desired warming level (warming level approach)	`UNSET`
`station_name`	`str`	Long name of desired station	`UNSET`
`latitude`	`float \| int(optional)`	Latitude for persistence XMY data if station_name not set	`UNSET`
`longitude`	`float \| int(optional)`	Longitude for persistence XMY data if station_name not set	`UNSET`
`verbose`	`bool`	True to increase verbosity	`True`

Attributes:

Name	Type	Description
`q`	`float`	extreme quantile, ranging from 0 to 1
`start_year`	`str`	Initial year of persistence XMY period
`end_year`	`str`	Final year of persistence XMY period
`warming_level`	`float \| int`	Warming level value
`lat_range`	`tuple`	Pair of latitudes that bracket `latitude`
`lon_range`	`tuple`	Pair of longitudes that bracket `longitude`
`reanalysis`	`bool`	True to use ERA5 reanalysis
`simulations`	`list[str]`	List of included simulations
`scenario`	`list[str]`	List of scenarios
`vars_and_units`	`dict[str, str]`	Dictionary of all required variables and units
`verbose`	`bool`	True to increase verbosity
`top_hours`	`DataFrame`	Table of top months by model
`all_vars`	`Dataset`	All loaded variables for persistence XMY
`air_temp`	`Dataset`	Air temperature, for use in finding candidate hours
`xmy_data_to_export`	`dict[Dataframe]`	Dictionary of persistence XMY results by simulation
`_skip_last`	`bool`	Internal flag to track last year for warming level approach

Source code in climakitae/explore/extreme_meteorological_year.py

def __init__(
    self,
    q: float = UNSET,
    start_year: int = UNSET,
    end_year: int = UNSET,
    warming_level: float | int = UNSET,
    station_name: str = UNSET,
    latitude: float | int = UNSET,
    longitude: float | int = UNSET,
    reanalysis: bool = False,
    verbose: bool = True,
):

    # Here we go through a few different ways to get the persistence XMY location
    match latitude, longitude, station_name:
        # UNSET will match to object type
        # Case 1: All variables set
        case float() | int(), float() | int(), str():
            if is_HadISD(station_name):
                raise ValueError(
                    "Do not set `latitude` and `longitude` when using a HadISD station for `station_name`. Change `station_name` value if using custom location."
                )
            else:
                print(
                    f"Initializing persistence XMY object for custom location: {latitude} N, {longitude} W with name '{station_name}'."
                )
                self._set_loc_from_lat_lon(latitude, longitude)
                self.stn_name = station_name
        # Case 2: lat/lon provided, no station_name string
        case float() | int(), float() | int(), object():
            print(
                f"Initializing persistence XMY object for custom location: {latitude} N, {longitude} W."
            )
            self._set_loc_from_lat_lon(latitude, longitude)
        # Case 3: station name provided, lat/lon not numeric
        case object(), object(), str():
            print(f"Initializing persistence XMY object for {station_name}.")
            self._set_loc_from_stn_name(station_name)
        # Last case: something else provided
        case _:
            raise ValueError(
                "No valid station name or latitude and longitude provided."
            )
    # Time variables
    match start_year, end_year, warming_level:
        # User set all options - bad
        case float() | int(), float() | int(), float() | int():
            raise ValueError(
                "Variables `start_year` and `end_year` cannot be paired with `warming_level`. Set either `start_year` and `end_year` OR `warming_level."
            )
        # Some other combo - unset variable will be saved as UNSET
        case _:
            self.start_year = start_year
            self.end_year = end_year
            self.warming_level = warming_level
            if isinstance(self.warming_level, int):
                self.warming_level = float(self.warming_level)
    # percentile value
    if 0 <= q <= 1:
        self.q = q
        p = q * 100
        p = int(p)
        print(f"Generating p{p} persistence XMY.")
    else:
        raise TypeError("Variable `q` must be a float between 0 and 1 (e.g. 0.5).")
    # Whether to drop the last month as a possible match
    self._skip_last = False
    if self.warming_level is not UNSET:
        # True for warming levels because final hours get lost to UTC conversion
        self._skip_last = True
    # Ranges used in get_data to pull smaller dataset without warnings
    self.lat_range = (self.stn_lat - 0.1, self.stn_lat + 0.1)
    self.lon_range = (self.stn_lon - 0.1, self.stn_lon + 0.1)
    # Determine if using models or ERA5 reanalysis
    if reanalysis:
        # Check that valid time period is requested
        if self.warming_level is not UNSET:
            raise ValueError(
                "Cannot use warming levels with ERA5 reanalysis. Please set start_year and end_year parameters."
            )
        if (self.start_year < 1981) or (self.end_year > 2019):
            raise ValueError(
                f"Valid start and end years for ERA5 reanalysis must be between 1981 and 2019. User provided start year {self.start_year} and end year {self.end_year}"
            )
        self.use_era5 = True
        self.simulations = ["WRF_ERA5_reanalysis"]
    else:  # use models
        self.use_era5 = False
        # These 4 simulations have the solar variables needed
        self.simulations = [
            "WRF_EC-Earth3_r1i1p1f1",
            "WRF_MPI-ESM1-2-HR_r3i1p1f1",
            "WRF_TaiESM1_r1i1p1f1",
            "WRF_MIROC6_r1i1p1f1",
        ]

    # Data only available for these scenarios
    self.scenario = ["Historical Climate", "SSP 3-7.0"]
    # Raw catalog variables to fetch (variable_id → display name)
    # These are fetched directly from the catalog in their native units.
    self._raw_vars = {
        "t2": "Air Temperature at 2m",
        "q2": "Water Vapor Mixing Ratio at 2m",
        "psfc": "Surface Pressure",
        "u10": "u10",
        "v10": "v10",
        "swdnb": "Instantaneous downwelling shortwave flux at bottom",
        "swddni": "Shortwave surface downward direct normal irradiance",
        "swddif": "Shortwave surface downward diffuse irradiance",
        "lwdnb": "Instantaneous downwelling longwave flux at bottom",
    }
    # Full set of persistence XMY variables (including derived) with desired units.
    # Used for display name references throughout the rest of the persistence XMY code.
    self.vars_and_units = {
        "Air Temperature at 2m": "degC",
        "Dew point temperature": "degC",
        "Relative humidity": "[0 to 100]",
        "Instantaneous downwelling shortwave flux at bottom": "W/m2",
        "Shortwave surface downward direct normal irradiance": "W/m2",
        "Shortwave surface downward diffuse irradiance": "W/m2",
        "Instantaneous downwelling longwave flux at bottom": "W/m2",
        "Wind speed at 10m": "m s-1",
        "Wind direction at 10m": "degrees",
        "Surface Pressure": "Pa",
        "Water Vapor Mixing Ratio at 2m": "g kg-1",
    }
    self.verbose = verbose
    # These will get set later in analysis
    self.top_hours = UNSET
    self.all_vars = UNSET
    self.air_temp_var = UNSET
    self.xmy_data_to_export = UNSET

`load_all_variables()`

Load hourly persistence XMY variables and derive daily statistics.

Fetches hourly raw variables via ClimateData for the 8760 profile assembly, then derives ALL daily statistics from the hourly data in local time. This approach avoids two problems with fetching daily catalog variables directly:

UTC vs local time: Catalog daily variables are pre-aggregated over UTC day boundaries, which differ from local-time days by the station's UTC offset (e.g., 8 hours for California).
Non-determinism: With a Dask distributed client active, lazy reductions (.resample().sum()) can produce slightly different floating-point results on each run due to non-deterministic task ordering. Computing hourly data eagerly to numpy before resampling guarantees deterministic daily statistics.

Source code in climakitae/explore/extreme_meteorological_year.py

def load_all_variables(self):
    """Load hourly persistence XMY variables and derive daily statistics.

    Fetches hourly raw variables via ClimateData for the 8760 profile
    assembly, then derives ALL daily statistics from the hourly data
    in local time.  This approach avoids two problems with fetching
    daily catalog variables directly:

    1. **UTC vs local time**: Catalog daily variables are pre-aggregated
       over UTC day boundaries, which differ from local-time days by the
       station's UTC offset (e.g., 8 hours for California).
    2. **Non-determinism**: With a Dask distributed client active, lazy
       reductions (``.resample().sum()``) can produce slightly different
       floating-point results on each run due to non-deterministic task
       ordering.  Computing hourly data eagerly to numpy before
       resampling guarantees deterministic daily statistics.
    """
    print("Loading data from catalog via ClimateData.")

    def _fetch_and_clean(variable_id, table_id="1hr"):
        """Fetch a single variable and clean coordinate cruft."""
        self._vprint(f"  Fetching {variable_id} ({table_id})...")
        da = self._fetch_raw_variable(variable_id, table_id=table_id)
        return da.squeeze().drop_vars(
            [
                "lakemask",
                "landmask",
                "x",
                "y",
                "Lambert_Conformal",
                "centered_year",
            ],
            errors="ignore",
        )

    # --- Hourly variables (needed for 8760 profile assembly) ---
    # First hourly fetch must be synchronous for warming level (sets year range)
    hourly_var_ids = list(self._raw_vars.keys())
    if self.warming_level is not UNSET:
        self._vprint(f"  Getting {hourly_var_ids[0]} (sets year range)...")
        # Fetch raw first to capture centered_year before cleaning drops it
        raw_first = self._fetch_raw_variable(hourly_var_ids[0], table_id="1hr")
        if "centered_year" in raw_first.coords:
            self._sim_centered_years = dict(
                zip(
                    raw_first.simulation.values,
                    raw_first.centered_year.values.ravel(),
                )
            )
        first_var = raw_first.squeeze().drop_vars(
            [
                "lakemask",
                "landmask",
                "x",
                "y",
                "Lambert_Conformal",
                "centered_year",
            ],
            errors="ignore",
        )
        first_var.name = self._raw_vars[hourly_var_ids[0]]
        remaining_hourly = hourly_var_ids[1:]
    else:
        first_var = None
        remaining_hourly = hourly_var_ids

    # --- Daily variables resampling ---
    # Derive daily stats from hourly data in local time.
    # This matches the original code's behavior (hourly → local time →
    # daily resample) and avoids using catalog daily variables which are
    # pre-aggregated over UTC day boundaries (different 24-hour window).
    # We also compute hourly data eagerly here so all downstream
    # operations use deterministic numpy math, avoiding non-deterministic
    # floating-point reduction ordering from the dask distributed scheduler.

    # Fetch remaining hourly variables in parallel
    self._vprint("  Loading hourly variables in parallel...")

    def _fetch_hourly(vid):
        da = _fetch_and_clean(vid, "1hr")
        da.name = self._raw_vars[vid]
        return ("hourly", vid, da)

    with ThreadPoolExecutor(max_workers=4) as executor:
        hourly_futures = [
            executor.submit(_fetch_hourly, vid) for vid in remaining_hourly
        ]
        all_results = [f.result() for f in hourly_futures]

    # Separate hourly results
    hourly_list = [r[2] for r in all_results]
    if first_var is not None:
        hourly_list = [first_var] + hourly_list

    # --- Build hourly dataset for 8760 profile ---
    self._vprint("  Merging raw hourly data.")
    raw_ds = xr.merge(hourly_list)

    # Compute hourly derived variables
    self._vprint("  Computing hourly derived variables.")
    t2_degc = raw_ds["Air Temperature at 2m"] - 273.15
    t2_degc.attrs["units"] = "degC"

    q2_gkg = raw_ds["Water Vapor Mixing Ratio at 2m"] * 1000
    q2_gkg.attrs["units"] = "g kg-1"

    psfc_hpa = raw_ds["Surface Pressure"] / 100.0
    psfc_hpa.attrs["units"] = "hPa"

    rh = compute_relative_humidity(
        pressure=psfc_hpa,
        temperature=t2_degc,
        mixing_ratio=q2_gkg,
    )
    rh.name = "Relative humidity"
    rh.attrs["units"] = "[0 to 100]"

    dew_point_k = compute_dewpointtemp(
        temperature=raw_ds["Air Temperature at 2m"],
        rel_hum=rh,
    )
    dew_point = dew_point_k - 273.15
    dew_point.name = "Dew point temperature"
    dew_point.attrs["units"] = "degC"

    wind_speed = compute_wind_mag(u10=raw_ds["u10"], v10=raw_ds["v10"])
    wind_speed.name = "Wind speed at 10m"

    wind_dir = compute_wind_dir(u10=raw_ds["u10"], v10=raw_ds["v10"])
    wind_dir.name = "Wind direction at 10m"

    t2_out = t2_degc.copy()
    t2_out.name = "Air Temperature at 2m"
    t2_out.attrs["units"] = "degC"

    self.air_temp_var = t2_out
    self._vprint("   Air temperature stored for use in finding candidate hours.")

    q2_out = q2_gkg.copy()
    q2_out.name = "Water Vapor Mixing Ratio at 2m"

    derived_list = [t2_out, dew_point, rh, wind_speed, wind_dir, q2_out]
    keep_vars = [
        "Instantaneous downwelling shortwave flux at bottom",
        "Shortwave surface downward direct normal irradiance",
        "Shortwave surface downward diffuse irradiance",
        "Instantaneous downwelling longwave flux at bottom",
        "Surface Pressure",
    ]
    kept_from_raw = [raw_ds[v] for v in keep_vars]
    hourly_ds = xr.merge(derived_list + kept_from_raw)
    # In ERA5 reanalysis case with 1 sim, the "simulation" dim
    # got dropped when datasets were squeezed. Adding it back here.
    if "simulation" not in hourly_ds.dims:
        hourly_ds = hourly_ds.expand_dims("simulation")
    self._hourly_data = hourly_ds

    # --- Build daily dataset ---

    # Compute hourly data eagerly so all daily resampling uses
    # deterministic numpy math (not affected by dask scheduler ordering).
    # For a single grid cell this is small (~30yr × 8760hr × 4sims).
    self._vprint("  Computing hourly data for daily resampling...")
    with ProgressBar():
        hourly_computed = hourly_ds.compute()

    self._vprint("  Resampling hourly data to daily statistics...")
    # Temperature: daily max, min, mean (in degC, already converted above)
    daily_tmax = hourly_computed["Air Temperature at 2m"].resample(time="1D").max()
    daily_tmax.name = "Daily max air temperature"
    daily_tmax.attrs["frequency"] = "daily"

    daily_tmin = hourly_computed["Air Temperature at 2m"].resample(time="1D").min()
    daily_tmin.name = "Daily min air temperature"
    daily_tmin.attrs["frequency"] = "daily"

    daily_tmean = (
        hourly_computed["Air Temperature at 2m"].resample(time="1D").mean()
    )
    daily_tmean.name = "Daily mean air temperature"
    daily_tmean.attrs["frequency"] = "daily"

    # Dew point: daily max, min, mean
    daily_dp_max = (
        hourly_computed["Dew point temperature"].resample(time="1D").max()
    )
    daily_dp_max.name = "Daily max dewpoint temperature"
    daily_dp_max.attrs["frequency"] = "daily"

    daily_dp_min = (
        hourly_computed["Dew point temperature"].resample(time="1D").min()
    )
    daily_dp_min.name = "Daily min dewpoint temperature"
    daily_dp_min.attrs["frequency"] = "daily"

    daily_dp_mean = (
        hourly_computed["Dew point temperature"].resample(time="1D").mean()
    )
    daily_dp_mean.name = "Daily mean dewpoint temperature"
    daily_dp_mean.attrs["frequency"] = "daily"

    # Wind speed: daily max, mean
    daily_ws_max = hourly_computed["Wind speed at 10m"].resample(time="1D").max()
    daily_ws_max.name = "Daily max wind speed"
    daily_ws_max.attrs["frequency"] = "daily"

    daily_ws_mean = hourly_computed["Wind speed at 10m"].resample(time="1D").mean()
    daily_ws_mean.name = "Daily mean wind speed"
    daily_ws_mean.attrs["frequency"] = "daily"

    # Radiation: daily sums
    ghi_sum = (
        hourly_computed["Instantaneous downwelling shortwave flux at bottom"]
        .resample(time="1D")
        .sum()
    )
    ghi_sum.name = "Global horizontal irradiance"
    ghi_sum.attrs["frequency"] = "daily"

    dni_sum = (
        hourly_computed["Shortwave surface downward direct normal irradiance"]
        .resample(time="1D")
        .sum()
    )
    dni_sum.name = "Direct normal irradiance"
    dni_sum.attrs["frequency"] = "daily"

    daily_arrays = [
        daily_tmax,
        daily_tmin,
        daily_tmean,
        daily_dp_max,
        daily_dp_min,
        daily_dp_mean,
        daily_ws_max,
        daily_ws_mean,
        ghi_sum,
        dni_sum,
    ]

    self.all_vars = xr.merge(daily_arrays)
    self._vprint("  Daily statistics ready.")

`set_top_hours()`

Calculate top months dataframe.

Source code in climakitae/explore/extreme_meteorological_year.py

def set_top_hours(self):
    """Calculate top months dataframe."""

    self._vprint("Finding top hours.")
    self.top_hours = persistence_get_top_hours(
        self.air_temp_var, self.q, self._skip_last
    )
    if self.use_era5:
        # Simulation column gets set to None in this case - correct it
        self.top_hours = self.top_hours.assign(simulation=self.simulations[0])

`show_xmy_data_to_export(simulation)`

Show line plots of persistence XMY data for single model.

Parameters:

Name	Type	Description	Default
`simulation`	`str`	Simulation to display. Use simulation names as defined in self.simuations	required

Source code in climakitae/explore/extreme_meteorological_year.py

def show_xmy_data_to_export(self, simulation: str):
    """Show line plots of persistence XMY data for single model.

    Parameters
    ----------
    simulation: str
        Simulation to display.
        Use simulation names as defined in self.simuations

    """
    if self.xmy_data_to_export is UNSET:
        print("No XMY data generated.")
        print("Please run xmy.generate_xmy() to create XMY data for viewing.")
        return

    # WVMR not in final XMY dataframe
    fig_y = list(self.vars_and_units.keys())
    fig_y.remove("Water Vapor Mixing Ratio at 2m")

    self.xmy_data_to_export[simulation].plot(
        x="time",
        y=fig_y,
        title=f"Shock Extreme Meteorological Year ({simulation})",
        subplots=True,
        figsize=(10, 8),
        legend=True,
    )

`run_xmy_analysis()`

Generate persistence extreme meteorological year data.

Output will be a list of dataframes per simulation. Print statements throughout the function indicate progress.

Notes

Results are saved to the class variable xmy_data_to_export.

Source code in climakitae/explore/extreme_meteorological_year.py

def run_xmy_analysis(self):
    """Generate persistence extreme meteorological year data.

    Output will be a list of dataframes per simulation.
    Print statements throughout the function indicate progress.

    Notes
    -----
    Results are saved to the class variable `xmy_data_to_export`.
    """
    print("Assembling XMY data to export.")

    self._vprint("  STEP 1: Computing hourly data")

    # Use cached hourly data instead of re-downloading
    if not hasattr(self, "_hourly_data") or self._hourly_data is None:
        # Fallback: load from catalog if run_xmy_analysis called standalone
        self.load_all_variables()
    with ProgressBar():
        all_vars_ds = self._hourly_data.compute()

    # When skip_last=True (WL approach), persistence_get_top_hours padded
    # air_temp_var to N_ceil*8760 stamps so the partial trailing year's
    # Jan–Nov hours are still selectable. top_hours["year"] can therefore
    # contain the partial year's dummy label. Pad all_vars_ds to the same
    # N_ceil*8760 length (NaN-fill) so _make_8760_tables sees the same
    # year set and year_idx never goes out of bounds.
    if self._skip_last:
        _hours_per_year = 8760
        _n_total = all_vars_ds.sizes["time"]
        _n_full_years = -(-_n_total // _hours_per_year)  # ceiling
        _missing = _n_full_years * _hours_per_year - _n_total
        if _missing > 0:
            _times = pd.DatetimeIndex(all_vars_ds.time.values)
            _pad_start = _times[-1] + pd.Timedelta(hours=1)
            _pad_times = pd.date_range(start=_pad_start, periods=_missing, freq="h")
            _full_times = _times.append(_pad_times)
            all_vars_ds = all_vars_ds.reindex(time=_full_times)

    # Construct XMY
    p = self.q * 100
    p = int(p)
    self._vprint(
        f"\n  STEP 2: Calculating p{p} Persistence Extreme Meteorological Year per model simulation\n  Progress bar shows code looping through each month in the year.\n"
    )

    xmy_data_to_export = self._make_8760_tables(
        all_vars_ds, self.top_hours
    )  # Return dict of persistence XMY by simulation

    self._vprint("  Smoothing data at transitions between months.")
    self._vprint("  Dropping water vapor mixing ratio.")
    # Smooth transition hours
    for sim in xmy_data_to_export:
        xmy_data_to_export[sim] = xmy_data_to_export[sim].reset_index()
        xmy_data_to_export[sim] = self._smooth_month_transition_hours(
            xmy_data_to_export[sim]
        )

        # Mixing ratio was only needed for smoothing relative humidity,
        # so it can be dropped now.
        xmy_data_to_export[sim] = xmy_data_to_export[sim].drop(
            columns="Water Vapor Mixing Ratio at 2m"
        )

        # Add metadata columns needed by EPW header writer.
        # The new-core pipeline squeezes these dimensions, so they must be
        # re-attached before export.
        if self.warming_level is not UNSET:
            xmy_data_to_export[sim]["warming_level"] = self.warming_level
        else:
            if self.use_era5:
                xmy_data_to_export[sim]["scenario"] = "historical_reanalysis"
            else:
                xmy_data_to_export[sim]["scenario"] = "historical+ssp370"

    self.xmy_data_to_export = xmy_data_to_export
    self._vprint("persistence XMY analysis complete.")

`export_xmy_data(extension='epw')`

Write XMY data to EPW file.

Parameters:

Name	Type	Description	Default
`extension`	`str`	Desired file extension ('tmy','epw', or 'csv')	`'epw'`

Source code in climakitae/explore/extreme_meteorological_year.py

def export_xmy_data(self, extension: str = "epw"):
    """Write XMY data to EPW file.

    Parameters
    ----------
    extension: str
        Desired file extension ('tmy','epw', or 'csv')

    """
    print("Exporting persistence XMY to file.")
    for sim, _ in self.xmy_data_to_export.items():
        # Get right year range
        if self.warming_level is UNSET:
            years = (self.start_year, self.end_year)
            clean_sim = sim
        else:
            centered_year = int(self._sim_centered_years[sim])
            year1 = centered_year - 15
            year2 = centered_year + 14
            years = (year1, year2)
            # Append warming level descriptor to simulation name.
            # Legacy sim names no longer contain "_historical+ssp370"
            # after the new-core migration, so we append directly.
            wl_label = match_str_to_wl(self.warming_level)
            clean_sim = f"{sim}_{wl_label}"
        # Attach centered_year so CSV export can include it
        if self.warming_level is not UNSET:
            self.xmy_data_to_export[sim]["centered_year"] = centered_year
        clean_stn_name = (
            self.stn_name.replace(" ", "_").replace("(", "").replace(")", "")
        )
        p = self.q * 100
        p = int(p)
        filename = "persistence_xmy_p{0}_{1}_{2}".format(
            p,
            clean_stn_name,
            clean_sim,
        ).lower()
        write_tmy_file(
            filename,
            self.xmy_data_to_export[sim],
            years,
            self.stn_name,
            self.stn_code,
            self.stn_lat,
            self.stn_lon,
            self.stn_state,
            file_ext=extension,
            xmy_header={"extreme_type": p, "xmy_type": "PERSISTENCE"},
        )

`get_candidate_hours()`

Get top candidate hours.

This function can be used to view the candidate hours without running the entire persistence XMY workflow.

Source code in climakitae/explore/extreme_meteorological_year.py

def get_candidate_hours(self):
    """Get top candidate hours.

    This function can be used to view the candidate hours
    without running the entire persistence XMY workflow.
    """
    p = self.q * 100
    p = int(p)
    self._vprint(f"Getting top hours for p{p} persistence XMY.")
    self.set_top_hours()

`generate_xmy()`

Run the whole XMY workflow.

Source code in climakitae/explore/extreme_meteorological_year.py

def generate_xmy(self):
    """Run the whole XMY workflow."""
    # This runs the whole workflow at once
    print("Running persistence XMY workflow.")
    self.load_all_variables()
    self.get_candidate_hours()
    self.run_xmy_analysis()
    self.export_xmy_data()

`find_hot_cold_extreme_from_median(sub_month, sub_clim, target=0.5, extreme='cold')`

Identifies hottest or coldest year based on deviation from the median (CDF=0.5) temperature value across years.

Parameters:

Name	Type	Description	Default
`sub_month`	`DataArray`	dims: (year, data, bin_number)	required
`sub_clim`	`DataArray`	climatology with same structure	required
`target`	`float`	CDF level (default 0.5)	`0.5`
`extreme`	`str`	"cold" -> pick minimum deviation "hot" -> pick maximum deviation	`'cold'`

Returns:

Name	Type	Description
`results`	`list`	median values per year
`worst_year`	`scalar`	identified extreme year

Source code in climakitae/explore/extreme_meteorological_year.py

def find_hot_cold_extreme_from_median(
    sub_month: xr.DataArray,
    sub_clim: xr.DataArray,
    target: float = 0.5,
    extreme: str = "cold",  # "cold" or "hot"
) -> tuple[np.ndarray, np.ndarray, np.integer]:
    """
    Identifies hottest or coldest year based on deviation from the median (CDF=0.5)
    temperature value across years.

    Parameters
    ----------
    sub_month : xr.DataArray
        dims: (year, data, bin_number)
    sub_clim : xr.DataArray
        climatology with same structure
    target : float
        CDF level (default 0.5)
    extreme : str
        "cold" -> pick minimum deviation
        "hot"  -> pick maximum deviation

    Returns
    -------
    results : list
        median values per year
    worst_year : scalar
        identified extreme year
    """

    # -----------------------------
    # Step 1: climatology median (p50 reference)
    # -----------------------------
    clim_prob = sub_clim.sel(data="probability")
    clim_bins = sub_clim.sel(data="bins")

    clim_idx = np.abs(clim_prob - target).argmin(dim="bin_number")
    clim_05 = clim_bins.isel(bin_number=clim_idx).values

    # Step 2: per-year median extraction
    # -----------------------------
    yr_prob = sub_month.sel(data="probability")
    yr_bins = sub_month.sel(data="bins")

    idx_list = np.abs(yr_prob - target).argmin(dim="bin_number").values.tolist()

    results = []

    years = sub_month["year"].values

    for i, yr in enumerate(years):
        idx = idx_list[i]

        val = yr_bins.sel(year=yr).isel(bin_number=idx).values
        results.append(val)

    results = np.array(results)

    # -----------------------------
    # Step 3: compute anomalies vs climatology
    # -----------------------------
    anomaly = results - clim_05

    # -----------------------------
    # Step 4: pick extreme year
    # -----------------------------
    if extreme == "cold":
        worst_idx = np.argmin(anomaly)  # most negative deviation
    elif extreme == "hot":
        worst_idx = np.argmax(anomaly)  # most positive deviation
    else:
        raise ValueError(f"Extreme must be 'cold' or 'hot', Received {extreme}.")

    worst_year = years[worst_idx]

    return results, anomaly, worst_year

`generate_candidate_months(cdf_monthly, cdf_climatology, extreme='cold', skip_last=False)`

Run find_hot_cold_extreme_from_median() over entire input dataset. Generate a dataframe of selected years per month and simulation for input to XMY generation function.

Parameters:

Name	Type	Description	Default
`cdf_monthly`	`DataArray`	Monthly CDF to compare against climatological CDF	required
`cdf_climatology`	`DataArray`	CDF representing climatology	required
`extreme`	`str`	The type of shock XMY. "cold" -> pick minimum deviation "hot" -> pick maximum deviation	`'cold'`
`skip_last`	`bool`	True to exclude the final month, e.g. if data missing after time conversion	`False`

Returns:

Name	Type	Description
`top_df`	`dataframe of selected years per simulation and month`

Source code in climakitae/explore/extreme_meteorological_year.py

def generate_candidate_months(
    cdf_monthly: xr.DataArray,
    cdf_climatology: xr.DataArray,
    extreme: str = "cold",
    skip_last: bool = False,  # "cold" or "hot"
) -> pd.DataFrame:
    """
    Run find_hot_cold_extreme_from_median() over entire input dataset.
    Generate a dataframe of selected years per month and simulation
        for input to XMY generation function.

    Parameters
    ----------
    cdf_monthly : xr.DataArray
        Monthly CDF to compare against climatological CDF
    cdf_climatology : xr.DataArray
        CDF representing climatology
    extreme : str
        The type of shock XMY.
        "cold" -> pick minimum deviation
        "hot"  -> pick maximum deviation
    skip_last: bool
        True to exclude the final month, e.g. if data missing after time conversion

    Returns
    -------
    top_df: dataframe of selected years per simulation and month
    """

    # select priority variable based on shock type
    if extreme == "cold":
        var = "Daily min air temperature"
    elif extreme == "hot":
        var = "Daily max air temperature"

    # subset CDFs for priority variable
    subset_clim = cdf_climatology[var]
    subset_month = cdf_monthly[var]

    all_months = cdf_monthly.month.values
    last_month = int(all_months[-1])

    results = []

    # iterate over all simulations and months
    for sim in cdf_monthly.simulation.values:
        clim_sim = subset_clim.sel(simulation=sim)
        month_sim = subset_month.sel(simulation=sim)
        for mon in cdf_monthly.month.values:
            clim_mon = clim_sim.sel(month=mon)
            month_mon = month_sim.sel(month=mon)

            if skip_last and int(mon) == last_month:
                # Drop the final year entirely so it cannot be selected
                month_mon = month_mon.sel(year=month_mon.year.values[:-1])

            _, _, worst_year = find_hot_cold_extreme_from_median(
                month_mon, clim_mon, extreme=extreme
            )
            row = {"month": mon, "simulation": sim, "year": worst_year}
            results.append(row)
    top_df = pd.DataFrame(results)

    return top_df

`persistence_get_top_hours(data, q=0.5, skip_last=False)`

Selects a representative year for each hour using the Cal-Adapt Standard Year methodology, using hourly air temperature.

Computes the specified quantile for each hour of the year across all years, then selects the actual data value closest to that quantile (not interpolated), and returns the year that value was taken from for each of the 8760 hours (one year), for each simulation.

Constructs a table of years or per hour per simulation for use in persistence XMY generation.

Parameters:

Name	Type	Description	Default
`data`	`DataArray`	Hourly air temperature with time and simulation dimensions.	required
`q`	`float`	Quantile value for selecting representative values (0.0 to 1.0). Default is 0.5 (median).	`0.5`
`skip_last`	`bool`	True to exclude the final month, e.g. if data missing after time conversion	`False`

Returns:

Type	Description
`DataFrame`	Standard year table for each warming level and simulation, with days of year as the index and hour of day as the columns. Multi-index columns include Hour, Warming_Level, and Simulation dimensions.

Source code in climakitae/explore/extreme_meteorological_year.py

def persistence_get_top_hours(
    data: xr.DataArray, q: float = 0.5, skip_last: bool = False
) -> pd.DataFrame:
    """
    Selects a representative year for each hour using the Cal-Adapt Standard Year
    methodology, using hourly air temperature.

    Computes the specified quantile for each hour of the year across all years,
    then selects the actual data value closest to that quantile (not interpolated),
    and returns the year that value was taken from for each of the 8760 hours (one year),
    for each simulation.

    Constructs a table of years or per hour per simulation for use in persistence XMY generation.

    Parameters
    ----------
    data : xr.DataArray
        Hourly air temperature with time and simulation dimensions.

    q : float, optional
        Quantile value for selecting representative values (0.0 to 1.0).
        Default is 0.5 (median).

    skip_last: bool
        True to exclude the final month, e.g. if data missing after time conversion

    Returns
    -------
    pd.DataFrame
        Standard year table for each warming level and simulation,
        with days of year as the index and hour of day as the columns.
        Multi-index columns include Hour, Warming_Level, and Simulation dimensions.

    """

    # Check for simulation dimension
    has_simulation = "simulation" in data.dims
    if has_simulation:
        simulations = data.simulation.values
    else:
        simulations = [None]

    if skip_last:
        # WL+local-time leaves the trailing year |offset_hours| short of 8760
        # (e.g. KLAX/PST: last year ends Dec 31 15:00). Pad the tail back to a
        # clean N*8760 cycle, then NaN-mask the trailing month so it can't be
        # selected. NaN (not inf) — nanquantile/nanargmin both skip NaN.
        times = pd.DatetimeIndex(data.time.values)
        n_full_years = -(-len(times) // 8760)  # ceiling division
        missing = n_full_years * 8760 - len(times)

        da_work = data.copy(deep=True)
        if missing > 0:
            pad_start = times[-1] + pd.Timedelta(hours=1)
            pad_times = pd.date_range(start=pad_start, periods=missing, freq="h")
            pad_shape = list(data.shape)
            time_axis = data.dims.index("time")
            pad_shape[time_axis] = len(pad_times)
            pad_coords = {d: data.coords[d] for d in data.dims if d != "time"}
            pad_coords["time"] = pad_times
            pad = xr.DataArray(
                np.full(pad_shape, np.nan), dims=data.dims, coords=pad_coords
            )
            da_work = xr.concat([da_work, pad], dim="time")

        # NaN-mask the entire trailing local-time month
        last_t = pd.DatetimeIndex(da_work.time.values)[-1]
        last_year, last_month = int(last_t.year), int(last_t.month)
        mask = (da_work.time.dt.year == last_year) & (
            da_work.time.dt.month == last_month
        )
        da_work = da_work.where(~mask)
    else:
        da_work = data

    # Get all available time data
    hours_per_year = 8760
    total_hours = len(da_work.time)
    n_years = total_hours // hours_per_year

    print(f"      📊 Processing {total_hours:,} hours ({n_years} years) of data")
    print(f"      🎯 Computing {q*100:.0f}th percentile for each hour of year")

    # Create hour-of-year coordinate for all data (cycling through 1-8760)
    hour_of_year_all = np.tile(np.arange(1, hours_per_year + 1), n_years)[:total_hours]
    da_work = da_work.assign_coords(hour_of_year=("time", hour_of_year_all))

    # Initialize storage for profiles
    df_list = []

    for sim_idx, sim in enumerate(simulations):
        # Select data for this warming level and simulation combination
        if has_simulation:
            subset_data = da_work.isel(simulation=sim_idx)
        else:
            subset_data = da_work

        # Vectorized quantile computation using numpy
        # Reshape raw values into (n_years, hours_per_year) then compute
        # the quantile across years for each hour-of-year position
        values = subset_data.values
        n_total = len(values)

        usable = (n_total // hours_per_year) * hours_per_year
        year_hour_matrix = values[:usable].reshape(-1, hours_per_year)

        # Compute quantile targets for each of the 8760 hour positions
        quantile_targets = np.nanquantile(year_hour_matrix, q, axis=0)  # shape: (8760,)

        # For each hour position, find the actual year whose value is
        # closest to the quantile (avoids interpolation)
        diffs = np.abs(
            year_hour_matrix - quantile_targets[np.newaxis, :]
        )  # (n_years, 8760)

        closest_year_idx = np.nanargmin(diffs, axis=0)  # (8760,)

        # Calendar year for each hour-of-year position
        usable_times = subset_data.time.values[:usable]
        row_years = pd.DatetimeIndex(usable_times).year.values.reshape(
            -1, hours_per_year
        )[:, 0]
        sel_year = row_years[closest_year_idx]  # (8760,) — calendar year, not row index

        # Store the selected years in a pandas dataframe
        df_i = pd.DataFrame(
            {
                "hour": np.arange(1, hours_per_year + 1),
                "simulation": sim,
                "year": sel_year.astype(int),
            }
        )
        df_list.append(df_i)

    # Concatenate list together for all simulations
    top_df = pd.concat(df_list).reset_index(drop=True)

    return top_df

Standard Year Profile

Calculates the Standard Year Climate Profiles using a warming level approach and designated quantiles. The historical baseline for relative profile computation is a warming level of 1.2 C. User specified warming level will be calculated relative to this baseline unless the "no_delta" option is set to True, in which case the raw profile(s) for the requested warming level(s) will be returned.

`export_profile_to_csv(profile, **kwargs)`

Export profile to csv file with a descriptive file name.

Each warming level is saved in a separate file.

Parameters:

Name	Type	Description	Default
`profile`	`DataFrame`	Standard year profile with MultiIndex columns	required
`**kwargs`	`dict`	Keyword arguments for data selection. Allowed keys: variable : str Name of variable used in profile q : float Percentile used in profile global_warming_levels : list[float] List of global warming levels in profile warming_level_winow: int in range (5,25), optional Years around Global Warming Level (+/-) (e.g. 15 means a 30yr window) location : str, tuple(float \| int), list[str \| float \| int] Station name, cached_area, or coordinates cached_area : str, optional Name of cached area used in profile no_delta : bool, default False, optional True if no_delta=True when generating profile approach (Optional) : str, "Warming Level" or "Time" The climate profile approach to use, either - "Time" (which is actually a warming level approach, but centered on an input) or - "Warming Level" (default) centered_year (Optional) : int in range [1980,2099] For approach="Time", the year for which to find a corresponding warming level time_profile_scenario (Optional) : str, default "SSP 3-7.0" SSP scenario from ["SSP 3-7.0", "SSP 2-4.5","SSP 5-8.5"] bias_adjusted_models (optional) : bool, default False, if True only return bias-adjusted WRF models	`{}`

Source code in climakitae/explore/standard_year_profile.py

def export_profile_to_csv(profile: pd.DataFrame, **kwargs: Any) -> None:
    """
    Export profile to csv file with a descriptive file name.

    Each warming level is saved in a separate file.

    Parameters
    ----------
    profile: pd.DataFrame
        Standard year profile with MultiIndex columns

    **kwargs : dict
        Keyword arguments for data selection. Allowed keys:
            variable : str
                Name of variable used in profile
            q : float
                Percentile used in profile
            global_warming_levels : list[float]
                List of global warming levels in profile
            warming_level_winow: int in range (5,25), optional
                Years around Global Warming Level (+/-) (e.g. 15 means a 30yr window)
            location : str, tuple(float | int), list[str | float | int]
                Station name, cached_area, or coordinates
            cached_area : str, optional
                Name of cached area used in profile
            no_delta : bool, default False, optional
                True if no_delta=True when generating profile
            approach (Optional) : str, "Warming Level" or "Time"
                The climate profile approach to use, either
                    - "Time" (which is actually a warming level approach, but centered on an input) or
                    - "Warming Level" (default)
            centered_year (Optional) : int in range [1980,2099]
                For approach="Time", the year for which to find a corresponding warming level
            time_profile_scenario (Optional) : str, default "SSP 3-7.0"
                SSP scenario from ["SSP 3-7.0", "SSP 2-4.5","SSP 5-8.5"]
            bias_adjusted_models (optional) : bool, default False, if True only return bias-adjusted WRF models
    """

    # Get required parameter values
    variable = kwargs.get("variable")
    q = kwargs.get("q")

    # Get warming_level, no_delta, warming_level_window, approach, centered_year inputs, and scenario
    no_delta = kwargs.get("no_delta", False)
    warming_level_window = kwargs.get("warming_level_window", None)
    approach = kwargs.get("approach", None)
    centered_year = kwargs.get("centered_year", None)
    global_warming_levels = kwargs.get("warming_level", None)
    scenario = kwargs.get("time_profile_scenario", None)
    ba_models = kwargs.get("bias_adjusted_models", False)

    # Get variable id string to use in file name
    variable_descriptions = read_csv_file(VARIABLE_DESCRIPTIONS_CSV_PATH)
    var_id = variable_descriptions[
        (variable_descriptions["display_name"] == variable)
        & (variable_descriptions["timescale"] == "hourly")
    ]["variable_id"].item()

    # Get location string based on combination of location variables
    kwargs = _handle_location_params(**kwargs)
    func_list = [_check_cached_area, _check_lat_lon, _check_stations]
    location_str = ""
    for func in func_list:
        location_str = func(location_str, **kwargs)

    # Check profile MultiIndex to pull out data by Global Warming Level
    match profile.keys().nlevels:
        case 1:  # Single WL
            # If 'warming_level' provided, fetch the value within the input list
            if global_warming_levels:
                gwl = global_warming_levels[0]
            else:
                # For time-based approach, do not include gwl in the filename - set it to None
                if approach == "Time":
                    gwl = global_warming_levels
                # Otherwise, use the default
                else:
                    gwl = 1.2

            filename = _get_clean_standardyr_filename(
                var_id,
                q,
                location_str,
                gwl,
                warming_level_window,
                no_delta,
                approach,
                centered_year,
                scenario,
                ba_models,
            )
            profile.to_csv(filename, index=False)
        case 2:  # Multiple WL (WL included in MultiIndex)
            for gwl in global_warming_levels:  # Single file per WL
                filename = _get_clean_standardyr_filename(
                    var_id,
                    q,
                    location_str,
                    gwl,
                    warming_level_window,
                    no_delta,
                    approach,
                    centered_year,
                    scenario,
                    ba_models,
                )
                profile.xs(f"WL_{gwl}", level="Warming_Level", axis=1).to_csv(
                    filename, index=False
                )
        case _:
            raise ValueError(
                f"Profile MultiIndex should have one or two levels. Found {profile.keys().nlevels} levels."
            )

`retrieve_profile_data(**kwargs)`

Backend function for retrieving data needed for computing climate profiles.

Reads in the full hourly data for the 8760 analysis, including all warming levels.

Parameters:

Name	Type	Description	Default
`**kwargs`	`dict`	Keyword arguments for data selection. Allowed keys: - variable (Optional) : str, default "Air Temperature at 2m" - resolution (Optional) : str, default "3 km" - approach (Optional) : str, "Warming Level" or "Time" - centered (Optional) : int - warming_levels (Optional) : List[float], default [1.2] - warming_level_window (Optional): int in range [5,25] - cached_area (Optional) : str or List[str] - latitude (Optional) : float or tuple - longitude (Optional) : float or tuple - stations (Optional) : list[str], default None - units (Optional) : str, default "degF" - no_delta (optional) : bool, default False, if True, do not retrieve historical data, return raw future profile - bias_adjusted_models (optional) : bool, default False, if True only return bias-adjusted WRF models	`{}`

Returns:

Type	Description
`Tuple[Dataset, Dataset]`	(historic_data, future_data, get_data_params) - Historical data at 1.2°C warming, and future data at specified warming levels.

Raises:

Type	Description
`ValueError`	If invalid parameter keys are provided.

Example

historic_data, future_data = retrieve_profile_data( ... variable="Air Temperature at 2m", ... resolution="45 km", ... time_profile_scenario="SSP 2-4.5", ... warming_level=[1.5, 2.0, 3.0], ... units="degF" ... )

historic_data, future_data = retrieve_profile_data( ... warming_level=[2.0] ... )

Notes

Historical data is always retrieved for warming level = 1.2°C. Future data uses user-specified warming levels or defaults.

The function prioritizes location parameters in the following order: 1. cached_area 2. latitude/longitude 3. stations Each parameter will override the lower-priority ones if provided. So if cached_area is given, lat/lon and stations are ignored. If lat/lon are given, stations are ignored. If stations are given, they are used only if neither cached_area nor lat/lon are provided.

If no location parameters are provided, a warning is issued about retrieving the entire CA dataset.

Source code in climakitae/explore/standard_year_profile.py

def retrieve_profile_data(**kwargs: Any) -> Tuple[xr.DataArray, xr.DataArray]:
    """
    Backend function for retrieving data needed for computing climate profiles.

    Reads in the full hourly data for the 8760 analysis, including all warming levels.

    Parameters
    ----------
    **kwargs : dict
        Keyword arguments for data selection. Allowed keys:
        - variable (Optional) : str, default "Air Temperature at 2m"
        - resolution (Optional) : str, default "3 km"
        - approach (Optional) : str, "Warming Level" or "Time"
        - centered (Optional) : int
        - warming_levels (Optional) : List[float], default [1.2]
        - warming_level_window (Optional): int in range [5,25]
        - cached_area (Optional) : str or List[str]
        - latitude (Optional) : float or tuple
        - longitude (Optional) : float or tuple
        - stations (Optional) : list[str], default None
        - units (Optional) : str, default "degF"
        - no_delta (optional) : bool, default False, if True, do not retrieve historical data, return raw future profile
        - bias_adjusted_models (optional) : bool, default False, if True only return bias-adjusted WRF models

    Returns
    -------
    Tuple[xr.Dataset, xr.Dataset]
        (historic_data, future_data, get_data_params) - Historical data at 1.2°C warming,
        and future data at specified warming levels.
    Raises
    ------
    ValueError
        If invalid parameter keys are provided.

    Example
    -------
    >>> historic_data, future_data = retrieve_profile_data(
    ...     variable="Air Temperature at 2m",
    ...     resolution="45 km",
    ...     time_profile_scenario="SSP 2-4.5",
    ...     warming_level=[1.5, 2.0, 3.0],
    ...     units="degF"
    ... )

    >>> historic_data, future_data = retrieve_profile_data(
    ...     warming_level=[2.0]
    ... )

    Notes
    -----
    Historical data is always retrieved for warming level = 1.2°C.
    Future data uses user-specified warming levels or defaults.

    The function prioritizes location parameters in the following order:
    1. cached_area
    2. latitude/longitude
    3. stations
    Each parameter will override the lower-priority ones if provided. So if cached_area
    is given, lat/lon and stations are ignored. If lat/lon are given, stations are
    ignored. If stations are given, they are used only if neither cached_area nor lat/lon
    are provided.

    If no location parameters are provided, a warning is issued about retrieving the
    entire CA dataset.
    """
    no_delta = kwargs.pop("no_delta", False)
    # Define allowed inputs with types and defaults
    # Compute units default separately to avoid runtime evaluation in dictionary
    units_default = (
        "degF"  # Default to degF if user hasn't specified both variable and units
        if kwargs.get("variable", None) is None and kwargs.get("units", None) is None
        else None  # otherwise default to None and let get_data decide
    )

    ALLOWED_INPUTS = {
        "variable": (str, "Air Temperature at 2m"),
        "resolution": (str, "3 km"),
        "approach": (str, "Warming Level"),
        "centered_year": (int, None),
        "time_profile_scenario": (str, "SSP 3-7.0"),
        "warming_level": (list, [1.2]),
        "warming_level_window": (int, None),
        "cached_area": ((str, list), None),
        "latitude": ((float, tuple), None),
        "longitude": ((float, tuple), None),
        "stations": (list, None),
        "units": (str, units_default),
        "bias_adjusted_models": (bool, False),
    }

    # if the user does not enter warming level the analysis is a moot point
    # because the historical data is always at 1.2C
    REQUIRED_INPUTS = []
    for req in REQUIRED_INPUTS:
        if req not in kwargs:
            raise ValueError(f"Missing required input: '{req}'")

    # Validate input keys
    invalid_keys = set(kwargs.keys()) - set(ALLOWED_INPUTS.keys())
    if invalid_keys:
        raise ValueError(
            f"Invalid input(s): {list(invalid_keys)}. "
            f"Allowed inputs are: {list(ALLOWED_INPUTS.keys())}"
        )

    # Validate input types
    for key, value in kwargs.items():
        expected_type, _ = ALLOWED_INPUTS[key]
        # Handle union types (tuples of types)
        if isinstance(expected_type, tuple):
            if not isinstance(value, expected_type):
                type_names = [t.__name__ for t in expected_type]
                raise TypeError(
                    f"Parameter '{key}' must be of type {' or '.join(type_names)}, "
                    f"got {type(value).__name__}"
                )
        else:
            if not isinstance(value, expected_type):
                raise TypeError(
                    f"Parameter '{key}' must be of type {expected_type.__name__}, "
                    f"got {type(value).__name__}"
                )
        # check that warming_level_window is between 5 and 25
        if key == "warming_level_window":
            if value not in range(5, 26):
                raise ValueError(
                    f"Parameter '{key}' must be an integer between 5 and 25, "
                    f"got {value}"
                )
        # check that time_profile_scenario is within ["SSP 3-7.0", "SSP 2-4.5","SSP 5-8.5"]
        if key == "time_profile_scenario":
            if value not in ["SSP 3-7.0", "SSP 2-4.5", "SSP 5-8.5"]:
                raise ValueError(
                    f"Parameter '{key}' must be 'SSP 3-7.0', 'SSP 2-4.5', or 'SSP 5-8.5', "
                    f"received {value}."
                )

    # Validate and update approach parameters
    kwargs = _handle_approach_params(**kwargs)

    # Validate location parameters
    # the bahavior will be to use cached_area if provided
    # otherwise use lat/lon if provided
    # otherwise use stations if provided
    location_params = ["cached_area", "latitude", "longitude", "stations"]
    provided_location_params = [
        key for key in location_params if kwargs.get(key) is not None
    ]

    if "cached_area" in provided_location_params:
        # If cached_area is provided, unset lat/lon and stations
        if "latitude" in kwargs or "longitude" in kwargs:
            kwargs.pop("latitude", None)
            kwargs.pop("longitude", None)
            print(
                "   ⚠️  Note: Using cached_area, ignoring provided latitude/longitude"
            )
        if "stations" in kwargs:
            kwargs.pop("stations", None)
            print("   ⚠️  Note: Using cached_area, ignoring provided stations")
    elif (
        "latitude" in provided_location_params
        or "longitude" in provided_location_params
    ):
        # If lat/lon provided, unset stations
        if "stations" in kwargs:
            kwargs.pop("stations", None)
            print("   ⚠️  Note: Using latitude/longitude, ignoring provided stations")
    elif "stations" in provided_location_params:
        # Stations provided - convert to lat/lon with buffer
        stations = kwargs.pop("stations")
        buffer = _get_buffer_from_resolution(kwargs.get("resolution", "3 km"))
        print(
            f"   📍 Converting {len(stations)} station(s) to lat/lon coordinates with ±{buffer:.2f}° buffer"
        )
        try:
            lat_bounds, lon_bounds = _convert_stations_to_lat_lon(
                stations, buffer=buffer
            )
            kwargs["latitude"] = lat_bounds
            kwargs["longitude"] = lon_bounds
            print(f"      Latitude range: {lat_bounds[0]:.4f} to {lat_bounds[1]:.4f}")
            print(f"      Longitude range: {lon_bounds[0]:.4f} to {lon_bounds[1]:.4f}")
        except ValueError as e:
            raise ValueError(f"Error converting stations to coordinates: {e}")
    else:
        # No location parameters provided - warn about entire CA dataset
        print(
            "   ⚠️  WARNING: No location parameter provided (cached_area, latitude/longitude, or stations)"
        )
        print(
            "      The entire California dataset will be retrieved, which may be very large and slow."
        )
        print(
            "      Consider specifying a cached_area, lat/lon bounds, or specific stations for better performance."
        )

    # Set default parameters for data retrieval
    # Note: if stations were provided, they've been converted to lat/lon above
    get_data_params = {
        "variable": kwargs.get("variable", "Air Temperature at 2m"),
        "resolution": kwargs.get("resolution", "3 km"),
        "downscaling_method": "Dynamical",  # must be WRF, cannot be LOCA
        "timescale": "hourly",  # must be hourly for 8760 analysis
        "area_average": "Yes",
        "units": kwargs.get(
            "units",
            (
                "degF"  # Default to degF if user hasn't specified both variable and units
                if kwargs.get("variable", None) is None
                and kwargs.get("units", None) is None
                else None  # otherwise default to None and let get_data decide
            ),
        ),
        "approach": "Warming Level",
        "warming_level": [1.2],  # Historic global warming level
        "warming_level_window": kwargs.get(
            "warming_level_window", 15
        ),  # Use user input warming level window, if provided. Otherwise, default to 15.
        "cached_area": kwargs.get("cached_area", None),
        "latitude": kwargs.get("latitude", None),
        "longitude": kwargs.get("longitude", None),
    }

    warming_level = kwargs.get("warming_level", [1.2])
    print(
        f"Step 2: Data is being retrieved for warming level {warming_level} across all available simulations and scenarios."
    )

    historic_data = None
    if not no_delta:
        # Retrieve historical data at 1.2°C warming level
        historic_data = get_data(**get_data_params)

    # Update with any user-provided parameters for future data retrieval
    get_data_params.update(kwargs)
    future_data = get_data(**get_data_params)

    # Filter for only bias-adjusted WRF models, if user indicates this
    ba_models = kwargs.get("bias_adjusted_models", False)
    if ba_models:
        print("Filtering data for bias-adjusted models.")
        future_data = _filter_ba_models(future_data)
        if historic_data is not None:
            historic_data = _filter_ba_models(historic_data)
    else:
        None

    print(
        f"Step 3: The climate profile will now be produced at warming level {warming_level}."
    )

    return historic_data, future_data

`get_climate_profile(**kwargs)`

High-level function to compute standard year climate profiles using warming level data.

This function retrieves climate data and computes standard year profiles using the 8760 analysis approach. It combines data retrieval and profile computation in a single call.

Parameters:

Name	Type	Description	Default
`**kwargs`	`dict`	Keyword arguments for data selection. Allowed keys: - variable (Optional) : str, default "Air Temperature at 2m" - resolution (Optional) : str, default "3 km" - approach (Optional) : str, "Warming Level" or "Time", default "Warming Level" - centered_year (Optional) : int - time_profile_scenario (Optional) : str, default "SSP 3-7.0" - warming_level (Required) : List[float], default [1.2] - warming_level_window (Optional): int in range [5,25], default 15 - location (Optional) : str, List[str\|float\|int], or Tuple[str\|float\|int] - units (Optional) : str, default "degF" - q (Optional) : float \| list[float], default 0.5, quantile for profile calculation - no_delta (optional) : bool, default False, if True, do not apply baseline subtraction, return raw future profile - bias_adjusted_models (optional) : bool, default False, if True only return bias-adjusted WRF models	`{}`

Returns:

Type	Description
`DataFrame`	Standard year table for each warming level, with days of year as the index and hour of day as the columns. If multiple warming levels exist, they will be included as additional column levels. Units and metadata are preserved in the DataFrame's attrs dictionary.

Examples:

>>> profile = get_climate_profile(
...     variable="Air Temperature at 2m",
...     warming_level=[1.5, 2.0, 3.0],
...     units="degF"
... )

>>> profile = get_climate_profile(warming_level=[2.0])

Source code in climakitae/explore/standard_year_profile.py

def get_climate_profile(**kwargs: Dict[str, Any]) -> pd.DataFrame:
    """
    High-level function to compute standard year climate profiles using warming level data.

    This function retrieves climate data and computes standard year
    profiles using the 8760 analysis approach. It combines data retrieval
    and profile computation in a single call.

    Parameters
    ----------
    **kwargs : dict
        Keyword arguments for data selection. Allowed keys:
        - variable (Optional) : str, default "Air Temperature at 2m"
        - resolution (Optional) : str, default "3 km"
        - approach (Optional) : str, "Warming Level" or "Time", default "Warming Level"
        - centered_year (Optional) : int
        - time_profile_scenario (Optional) : str, default "SSP 3-7.0"
        - warming_level (Required) : List[float], default [1.2]
        - warming_level_window (Optional): int in range [5,25], default 15
        - location (Optional) : str, List[str|float|int], or Tuple[str|float|int]
        - units (Optional) : str, default "degF"
        - q (Optional) : float | list[float], default 0.5, quantile for profile calculation
        - no_delta (optional) : bool, default False, if True, do not apply baseline subtraction, return raw future profile
        - bias_adjusted_models (optional) : bool, default False, if True only return bias-adjusted WRF models

    Returns
    -------
    pd.DataFrame
        Standard year table for each warming level, with days of year as
        the index and hour of day as the columns. If multiple warming levels exist,
        they will be included as additional column levels. Units and metadata are
        preserved in the DataFrame's attrs dictionary.

    Examples
    --------
    >>> profile = get_climate_profile(
    ...     variable="Air Temperature at 2m",
    ...     warming_level=[1.5, 2.0, 3.0],
    ...     units="degF"
    ... )

    >>> profile = get_climate_profile(warming_level=[2.0])
    """
    # Extract parameters for compute_profile
    q = kwargs.pop("q", 0.5)
    no_delta = kwargs.get("no_delta", False)

    # Parse the location parameter
    kwargs = _handle_location_params(**kwargs)

    # Retrieve the climate data
    print("📊 Retrieving climate data...")
    with tqdm(
        total=2 if not no_delta else 1, desc="Data retrieval", unit="dataset"
    ) as pbar:
        historic_data, future_data = retrieve_profile_data(**kwargs)
        pbar.update(2)

    # Notify users of default values being used in the absence of input parameters
    # relevant for warming_level_window and warming_level
    defaults = {
        "warming_level": [1.2],
        "warming_level_window": 15,
        "approach": "Warming Level",
        "variable": "Air Temperature at 2m",
        "q": 0.5,
        "resolution": "3 km",
        "units": "degF",
    }
    for key, default_val in defaults.items():
        if key in kwargs:
            # skip this key
            continue
        if key == "warming_level":
            # if approach=Time, then default warming level is not used
            if kwargs.get("approach") == "Time":
                continue
            else:
                print(f"Using default '{key}': {default_val}")
                kwargs[key] = default_val
        else:
            if key == "q" and q != 0.5:
                continue
            else:
                print(f"Using default '{key}': {default_val}")
                kwargs[key] = default_val

    # catch invalid selections that return None
    if future_data is None and historic_data is None:
        raise ValueError(
            "No data returned for either historical or future datasets.\nPlease review your data selection parameters."
        )

    # Call compute_profile with the processed data
    # Compute profiles for both historical and future data
    def _fetch_primary_data_variable(ds: xr.Dataset) -> xr.DataArray:
        """
        Helper to extract the primary data variable from a Dataset.

        Parameters
        ----------
        ds : xr.Dataset
            The xarray Dataset from which to extract the primary data variable.

        Returns
        -------
        xr.DataArray
            The primary data variable as a DataArray.
        """
        if not isinstance(ds, xr.Dataset):
            return ds
        var_name = list(ds.data_vars.keys())[0]
        return ds[var_name]

    future_profile_data = _fetch_primary_data_variable(future_data)
    historic_profile_data = _fetch_primary_data_variable(historic_data)

    # Compute profiles for both datasets
    print("⚙️  Computing climate profiles...")

    future_profile = compute_profile(future_profile_data, q=q)
    if no_delta:
        historic_profile = None
    else:
        historic_profile = compute_profile(historic_profile_data, q=q)

    if no_delta:
        print("   ✓ No baseline subtraction requested, returning raw future profile")
        return future_profile

    # Compute the difference profile
    difference_profile = _compute_difference_profile(future_profile, historic_profile)

    print(
        f"✅ Climate profile computation complete! Final shape: {difference_profile.shape}"
    )
    print(
        f"   (Days: {difference_profile.shape[0]}, Hours/Columns: {difference_profile.shape[1]})"
    )

    return difference_profile

`compute_profile(data, q=0.5)`

Calculates the standard year climate profile for warming level data using 8760 analysis.

This function handles global warming levels approach using time_delta coordinate. Processes all 30 years of warming level data centered around the year a warming level is reached, computes the specified quantile for each hour of the year across all years, then selects the actual data value closest to that quantile (not interpolated), and returns a characteristic profile of 8760 hours (one year) for each warming level and simulation combination.

Parameters:

Name	Type	Description	Default
`data`	`DataArray`	Hourly base-line subtracted data for one variable with warming_level, time_delta, and simulation dimensions. Expected to contain ~30 years (262,800 hours) of data for each warming level and simulation.	required
`q`	`float`	Quantile value for selecting representative values (0.0 to 1.0). Default is 0.5 (median).	`0.5`

Returns:

Type	Description
`DataFrame`	Standard year table for each warming level and simulation, with days of year as the index and hour of day as the columns. Multi-index columns include Hour, Warming_Level, and Simulation dimensions.

Source code in climakitae/explore/standard_year_profile.py

def compute_profile(data: xr.DataArray, q=0.5) -> pd.DataFrame:
    """
    Calculates the standard year climate profile for warming level data using 8760
    analysis.

    This function handles global warming levels approach using time_delta coordinate.
    Processes all 30 years of warming level data centered around the year a warming level
    is reached, computes the specified quantile for each hour of the year across all years,
    then selects the actual data value closest to that quantile (not interpolated),
    and returns a characteristic profile of 8760 hours (one year) for each warming level
    and simulation combination.

    Parameters
    ----------
    data : xr.DataArray
        Hourly base-line subtracted data for one variable with warming_level,
        time_delta, and simulation dimensions. Expected to contain ~30 years
        (262,800 hours) of data for each warming level and simulation.

    q : float, optional
        Quantile value for selecting representative values (0.0 to 1.0).
        Default is 0.5 (median).

    Returns
    -------
    pd.DataFrame
        Standard year table for each warming level and simulation,
        with days of year as the index and hour of day as the columns.
        Multi-index columns include Hour, Warming_Level, and Simulation dimensions.

    """
    # Check for simulation dimension
    has_simulation = "simulation" in data.dims
    if has_simulation:
        n_simulations = len(data.simulation)
        simulations = data.simulation.values
    else:
        n_simulations = 1
        simulations = [None]

    # Get all available time_delta data (all 30 years)
    hours_per_year = 8760
    total_hours = len(data.time_delta)
    n_years = total_hours // hours_per_year

    print(f"      📊 Processing {total_hours:,} hours ({n_years} years) of data")
    print(f"      🎯 Computing {q*100:.0f}th percentile for each hour of year")

    # Create hour-of-year coordinate for all data (cycling through 1-8760)
    hour_of_year_all = np.tile(np.arange(1, hours_per_year + 1), n_years)[:total_hours]
    data = data.assign_coords(hour_of_year=("time_delta", hour_of_year_all))

    warming_levels = data.warming_level.values

    # Create helper function to extract meaningful simulation labels
    def _get_simulation_label(sim: str | int | None, sim_idx: int) -> str:
        """Extract meaningful simulation label from simulation identifier."""
        if sim is None:
            return f"Sim_{sim_idx+1}"

        sim_str = str(sim)
        if "WRF_" in sim_str:
            # Parse simulation name format: WRF_GCM_params_scenario
            # Example: WRF_CESM2_r11i1p1f1_historical+ssp245
            parts = sim_str.split("_")
            if len(parts) >= 4:
                gcm = parts[1]  # e.g., CESM2, CNRM-ESM2-1
                params = parts[2]  # e.g., r11i1p1f1
                scenario = parts[3]  # e.g., historical+ssp245

                # Extract SSP from scenario (e.g., ssp245 from historical+ssp245)
                if "ssp" in scenario:
                    ssp_part = scenario.split("ssp")[-1]  # Get part after 'ssp'
                    ssp = f"ssp{ssp_part}"
                else:
                    ssp = "hist"  # fallback for historical-only

                return f"{gcm}-{params}-{ssp}"
            elif len(parts) >= 2:
                # Fallback for shorter format
                return f"{parts[1]}-{sim_idx+1}"
            else:
                return f"Sim_{sim_idx+1}"
        else:
            # Ensure uniqueness by adding index for non-WRF format
            base_name = sim_str.split("_")[0] if "_" in sim_str else sim_str
            return f"{base_name}-{sim_idx+1}"

    # Process all data using quantile computation across years
    print(
        f"      ⚙️ Computing quantiles for {len(warming_levels)} warming level(s) and {n_simulations} simulation(s)"
    )

    # Eagerly compute all dask data at once (one round-trip to scheduler)
    if hasattr(data.data, "chunks"):
        print("      📥 Loading data into memory...")
        from dask.diagnostics import ProgressBar

        with ProgressBar():
            data = data.compute()

    # Initialize storage for profiles
    profile_data = {}

    # Progress tracking
    total_combinations = len(warming_levels) * n_simulations
    with tqdm(
        total=total_combinations,
        desc="      Computing profiles",
        unit="combo",
        leave=False,
    ) as pbar:

        for wl_idx, wl in enumerate(warming_levels):
            for sim_idx, sim in enumerate(simulations):
                # Get simulation label
                sim_label = _get_simulation_label(sim, sim_idx)
                # Select data for this warming level and simulation combination
                if has_simulation:
                    subset_data = data.isel(warming_level=wl_idx, simulation=sim_idx)
                else:
                    subset_data = data.isel(warming_level=wl_idx)
                # Vectorized quantile computation using numpy
                # Reshape raw values into (n_years, hours_per_year) then compute
                # the quantile across years for each hour-of-year position
                values = subset_data.values
                n_total = len(values)
                usable = (n_total // hours_per_year) * hours_per_year
                year_hour_matrix = values[:usable].reshape(-1, hours_per_year)

                # Compute quantile targets for each of the 8760 hour positions
                quantile_targets = np.nanquantile(
                    year_hour_matrix, q, axis=0
                )  # shape: (8760,)

                # For each hour position, find the actual year whose value is
                # closest to the quantile (avoids interpolation)
                diffs = np.abs(
                    year_hour_matrix - quantile_targets[np.newaxis, :]
                )  # (n_years, 8760)
                closest_year_idx = np.nanargmin(diffs, axis=0)  # (8760,)
                profile_1d = year_hour_matrix[
                    closest_year_idx, np.arange(hours_per_year)
                ]

                # Store the profile
                key = (f"WL_{wl}", sim_label)
                profile_data[key] = profile_1d

                pbar.update(1)

    df_profile = _construct_profile_dataframe(
        profile_data=profile_data,
        warming_levels=warming_levels,
        simulations=simulations,
        sim_label_func=_get_simulation_label,
        hours_per_year=hours_per_year,
    )

    # Prepare metadata dictionary
    metadata = {
        "quantile": q,
        "method": "8760 analysis - actual data closest to quantile across 30 years",
        "description": f"Climate profile computed using actual data values closest to the {q*100:.0f}th percentile of hourly data",
    }

    # Add original data attributes if available
    if hasattr(data, "attrs"):
        if "units" in data.attrs:
            metadata["units"] = data.attrs["units"]
        if "extended_description" in data.attrs:
            metadata["extended_description"] = data.attrs["extended_description"]
        if "variable_id" in data.attrs:
            metadata["variable_name"] = data.attrs["variable_id"]
        elif hasattr(data, "name") and data.name:
            metadata["variable_name"] = data.name

    # Set all metadata using the helper function
    set_profile_metadata(df_profile, metadata)

    print(f"      ✅ Profile computation complete! Final shape: {df_profile.shape}")
    print(
        f"         With index: {df_profile.index.name}, columns: {df_profile.columns.names}"
    )
    if hasattr(data, "attrs") and "units" in data.attrs:
        print(f"         Units: {data.attrs['units']}")

    return df_profile

`get_profile_units(profile_df)`

Extract units information from a climate profile DataFrame.

Parameters:

Name	Type	Description	Default
`profile_df`	`DataFrame`	Climate profile DataFrame with units stored in attrs	required

Returns:

Type	Description
`str`	Units string, or 'Unknown' if not found

Examples:

>>> profile = get_climate_profile(variable="Air Temperature at 2m", warming_level=[2.0])
>>> units = get_profile_units(profile)
>>> print(f"Temperature units: {units}")

Source code in climakitae/explore/standard_year_profile.py

def get_profile_units(profile_df: pd.DataFrame) -> str:
    """
    Extract units information from a climate profile DataFrame.

    Parameters
    ----------
    profile_df : pd.DataFrame
        Climate profile DataFrame with units stored in attrs

    Returns
    -------
    str
        Units string, or 'Unknown' if not found

    Examples
    --------
    >>> profile = get_climate_profile(variable="Air Temperature at 2m", warming_level=[2.0])
    >>> units = get_profile_units(profile)
    >>> print(f"Temperature units: {units}")
    """
    return profile_df.attrs.get("units", "Unknown")

`get_profile_metadata(profile_df)`

Extract all metadata from a climate profile DataFrame.

Parameters:

Name	Type	Description	Default
`profile_df`	`DataFrame`	Climate profile DataFrame with metadata stored in attrs	required

Returns:

Type	Description
`dict`	Dictionary containing all available metadata

Examples:

>>> profile = get_climate_profile(variable="Air Temperature at 2m", warming_level=[2.0])
>>> metadata = get_profile_metadata(profile)
>>> print(f"Variable: {metadata.get('variable_name')}")
>>> print(f"Units: {metadata.get('units')}")
>>> print(f"Method: {metadata.get('method')}")

Source code in climakitae/explore/standard_year_profile.py

def get_profile_metadata(profile_df: pd.DataFrame) -> dict:
    """
    Extract all metadata from a climate profile DataFrame.

    Parameters
    ----------
    profile_df : pd.DataFrame
        Climate profile DataFrame with metadata stored in attrs

    Returns
    -------
    dict
        Dictionary containing all available metadata

    Examples
    --------
    >>> profile = get_climate_profile(variable="Air Temperature at 2m", warming_level=[2.0])
    >>> metadata = get_profile_metadata(profile)
    >>> print(f"Variable: {metadata.get('variable_name')}")
    >>> print(f"Units: {metadata.get('units')}")
    >>> print(f"Method: {metadata.get('method')}")
    """
    return dict(profile_df.attrs)

`set_profile_metadata(profile_df, metadata)`

Set or update metadata in a climate profile DataFrame.

Parameters:

Name	Type	Description	Default
`profile_df`	`DataFrame`	Climate profile DataFrame to update	required
`metadata`	`dict`	Dictionary containing metadata key-value pairs to set	required

Returns:

Type	Description
`None`	The function modifies the DataFrame in place

Examples:

>>> profile = get_climate_profile(variable="Air Temperature at 2m", warming_level=[2.0])
>>> new_metadata = {
...     "source": "Custom Dataset",
...     "author": "Jane Doe",
...     "notes": "This profile was generated for testing purposes."
... }
>>> set_profile_metadata(profile, new_metadata)
>>> print(profile.attrs)

Source code in climakitae/explore/standard_year_profile.py

def set_profile_metadata(profile_df: pd.DataFrame, metadata: dict) -> None:
    """
    Set or update metadata in a climate profile DataFrame.

    Parameters
    ----------
    profile_df : pd.DataFrame
        Climate profile DataFrame to update
    metadata : dict
        Dictionary containing metadata key-value pairs to set

    Returns
    -------
    None
        The function modifies the DataFrame in place

    Examples
    --------
    >>> profile = get_climate_profile(variable="Air Temperature at 2m", warming_level=[2.0])
    >>> new_metadata = {
    ...     "source": "Custom Dataset",
    ...     "author": "Jane Doe",
    ...     "notes": "This profile was generated for testing purposes."
    ... }
    >>> set_profile_metadata(profile, new_metadata)
    >>> print(profile.attrs)
    """
    if not isinstance(metadata, dict):
        raise ValueError("Metadata must be provided as a dictionary.")

    for key, value in metadata.items():
        profile_df.attrs[key] = value

Climate Vulnerability Assessment

Bases: Parameterized

Climate Analysis and Vulnerability Assessment Parameters Class.

This class defines and validates parameters for climate vulnerability analysis, supporting various climate variables, metrics, and analysis approaches.

Attributes:

Name	Type	Description
`input_locations`	`DataFrame`	Input coordinates that must have 'lat' and 'lon' columns with numeric data.
`time_start_year`	`int, default=1981`	Start year for the analysis period. Must be between 1900 and 2100.
`time_end_year`	`int, default=2010`	End year for the analysis period. Must be between 1900 and 2100.
`units`	`str, default="Celsius"`	Units for temperature measurement.
`variable`	`str, default="Air Temperature at 2m"`	Climate variable to analyze. Options include "Air Temperature at 2m", "Precipitation (total)", "NOAA Heat Index", "Effective Temperature".
`metric_calc`	`str, default="max"`	Statistical metric calculation method. Options: "min", "max", "mean", "median".
`percentile`	`float or None, default=None`	Percentile value for calculation. Must be between 0 and 100 if specified.
`heat_idx_threshold`	`float or None, default=None`	Threshold value for heat index calculations.
`one_in_x`	`int, float, list or None, default=None`	Return period(s) for extreme event analysis (e.g., 1-in-100 year event). Can be a single value or list of values.
`season`	`str, default="all"`	Season to analyze. Options: "summer", "winter", "all".
`downscaling_method`	`str, default="Dynamical"`	Climate model downscaling method. Options: "Dynamical", "Statistical".
`approach`	`str, default="Time"`	Analysis approach. Options: "Time" (time period), "Warming Level" (temperature target).
`warming_level`	`float, default=1.0`	Global warming level in degrees Celsius. Must be between 0 and 7.
`wrf_bias_adjust`	`bool, default=True`	Whether to apply bias adjustment to WRF model data.
`historical_data`	`str, default="Historical Climate"`	Type of historical data. Options: "Historical Climate", "Historical Reconstruction".
`ssp_data`	`list, default=["SSP 3-7.0"]`	Shared Socioeconomic Pathway scenarios to use. Options: "SSP 2-4.5", "SSP 3-7.0", "SSP 5-8.5".
`export_method`	`str, default="both"`	Data export method. Options: "raw", "calculate", "both", "None".
`separate_files`	`bool, default=False`	Whether to save climate variables in separate files.
`file_format`	`str, default="NetCDF"`	Output file format.
`batch_mode`	`bool, default=False`	Whether to run in batch processing mode.
`distr`	`str, default="gev"`	Statistical distribution for extreme value analysis. Options: "gev" (Generalized Extreme Value), "genpareto" (Generalized Pareto), "gamma".
`duration`	`tuple, default=UNSET`	Length of extreme event, specified as (4, 'hour'). Only implemented for hourly data.
`groupby`	`tuple, default=(1, "day")`	Group over which to look for max occurrence, specified as (1, 'day'). Only 'day' groupings supported.
`grouped_duration`	`tuple, default=UNSET`	Length of event after grouping, specified as (5, 'day'). Requires groupby to be set.
`file_name`	`str`	Base name for output files, no extension (e.g., "output", not "output.nc").

Methods:

Name Description

validate_params

Validate all parameters for consistency and compatibility.

get_names

Generate standardized names and metadata for data processing. If parameter validation fails. Common validation errors include: - Missing or non-numeric lat/lon columns in input_locations - Invalid time range (start year > end year) - Incompatible parameter combinations (e.g., multiple threshold parameters) - Unsupported variable-downscaling method combinations - Invalid approach-data type combinations

Notes

The class enforces several validation rules: - Only one threshold parameter (heat_idx_threshold, percentile, or one_in_x) can be specified - Historical Reconstruction data requires time-based approach and end year <= 2022 - NOAA Heat Index and Effective Temperature cannot use Statistical downscaling - Dynamical downscaling with time-based approach only supports SSP 3-7.0 - duration only supports "hour" units; groupby and grouped_duration only support "day" units

Examples:

>>> import pandas as pd
>>> locations = pd.DataFrame({'lat': [34.05, 36.17], 'lon': [-118.25, -115.14]})
>>> params = CavaParams(
...     input_locations=locations,
...     time_start_year=2020,
...     time_end_year=2050,
...     variable="Air Temperature at 2m",
...     percentile=95,
...     metric_calc="max"
... )

Source code in climakitae/explore/vulnerability.py

def __init__(self, **params):
    super().__init__(**params)
    self.validate_params()

`validate_params()`

Validate the parameters for vulnerability analysis.

Returns:

Type	Description
`None`

Raises:

Type	Description
`ValueError`	If any validation check fails. The error message lists all validation failures found during the validation process.

Notes

The method validates the following conditions: - Input locations DataFrame contains required 'lat' and 'lon' columns with numeric data types (float64 or int64) - Time range validity (start year must be <= end year) - Historical reconstruction data constraints: - End year must be <= 2022 - Only time-based approach is supported - Mutual exclusivity of threshold parameters (only one of heat_idx_threshold, percentile, or one_in_x can be specified) - At least one threshold parameter must be specified - Metric calculation ('min' or 'max') compatibility with percentile calculations - Variable and downscaling method compatibility (NOAA Heat Index and Effective Temperature cannot use Statistical downscaling) - SSP data constraints for WRF/Dynamical downscaling (only SSP 3-7.0 allowed for time-based approach)

Side Effects

Converts self.one_in_x to a list if it's not None and not already a list.

Source code in climakitae/explore/vulnerability.py

def validate_params(self):
    """
    Validate the parameters for vulnerability analysis.

    Returns
    -------
    None

    Raises
    ------
    ValueError
        If any validation check fails. The error message lists all validation
        failures found during the validation process.

    Notes
    -----
    The method validates the following conditions:
    - Input locations DataFrame contains required 'lat' and 'lon' columns with
      numeric data types (float64 or int64)
    - Time range validity (start year must be <= end year)
    - Historical reconstruction data constraints:
        - End year must be <= 2022
        - Only time-based approach is supported
    - Mutual exclusivity of threshold parameters (only one of heat_idx_threshold,
      percentile, or one_in_x can be specified)
    - At least one threshold parameter must be specified
    - Metric calculation ('min' or 'max') compatibility with percentile calculations
    - Variable and downscaling method compatibility (NOAA Heat Index and Effective
      Temperature cannot use Statistical downscaling)
    - SSP data constraints for WRF/Dynamical downscaling (only SSP 3-7.0 allowed
      for time-based approach)

    Side Effects
    ------------
    Converts self.one_in_x to a list if it's not None and not already a list.
    """

    errors = []

    # Change `one_in_x` type before param validation and arg generation
    one_in_x = self.one_in_x
    if one_in_x is not None:
        if not isinstance(one_in_x, list):
            self.one_in_x = [one_in_x]
        # Keep as list, don't convert to numpy array here

    # Check if input_locations is a DataFrame with lat/lon columns
    if isinstance(self.input_locations, pd.DataFrame):
        if (
            "lat" not in self.input_locations.columns
            or "lon" not in self.input_locations.columns
        ):
            errors.append("Input coordinates must have `lat` and `lon` columns.")
        elif not pd.api.types.is_numeric_dtype(
            self.input_locations["lat"]
        ) or not pd.api.types.is_numeric_dtype(self.input_locations["lon"]):
            errors.append("Input lat/lon columns must be float64 or int64 types.")

    # Validating that time start year < time end year
    if self.time_start_year > self.time_end_year:
        errors.append(
            "Start year must come before, or be the same year as, the end year."
        )

    # Validate year range for historical reconstruction
    if (
        self.historical_data == "Historical Reconstruction"
        and self.time_end_year > 2022
    ):
        errors.append(
            "End year for Historical Reconstruction data must be 2022 or earlier."
        )

    # Validate time approach with Historical Reconstruction
    if (
        self.historical_data == "Historical Reconstruction"
        and self.approach != "Time"
    ):
        errors.append(
            "Historical Reconstruction data can only be retrieved using a time-based approach."
        )

    # Validating that only one of heat_idx_threshold, percentile, or one_in_x is passed
    if (
        (self.heat_idx_threshold is not None and self.percentile is not None)
        or (self.heat_idx_threshold is not None and self.one_in_x is not None)
        or (self.percentile is not None and self.one_in_x is not None)
    ):
        errors.append(
            "Only one of heat_idx_threshold, one_in_x, and percentile can be non-None."
        )

    # Check that not all 3 customizable metric arguments are none.
    if (
        self.heat_idx_threshold is None
        and self.percentile is None
        and self.one_in_x is None
    ):
        errors.append(
            "`heat_idx_threshold`, `percentile`, and `one_in_x` cannot all be None."
        )

    # Make sure the metric is 'min' or 'max' if percentile is passed in
    if self.percentile is not None and self.metric_calc not in ["min", "max"]:
        errors.append(
            "Metric calculation must be 'min' or 'max' for percentile calculations."
        )

    # Validating variable and downscaling method combo
    if (
        self.variable in ["NOAA Heat Index", "Effective Temperature"]
        and self.downscaling_method == "Statistical"
    ):
        errors.append(
            f"`{self.variable}` cannot be used with Statistical downscaling method."
        )

    # Only allow for SSP 3-7.0 for WRF data
    if self.downscaling_method == "Dynamical" and self.approach == "Time":
        if len(self.ssp_data) > 1 or self.ssp_data[0] != "SSP 3-7.0":
            errors.append(
                "Can only use SSP 3-7.0 data for a time-based approach using WRF data."
            )

    # Validate duration, groupby, and grouped_duration params
    if self.duration is not UNSET:
        if (
            not isinstance(self.duration, tuple)
            or len(self.duration) != 2
            or self.duration[1] != "hour"
        ):
            errors.append(
                "`duration` must be a tuple with unit 'hour', e.g. (4, 'hour')."
            )

    if self.groupby is not UNSET:
        if (
            not isinstance(self.groupby, tuple)
            or len(self.groupby) != 2
            or self.groupby[1] != "day"
        ):
            errors.append(
                "`groupby` must be a tuple with unit 'day', e.g. (1, 'day')."
            )

    if self.grouped_duration is not UNSET:
        if self.groupby is UNSET:
            errors.append("To use `grouped_duration`, `groupby` must also be set.")
        if (
            not isinstance(self.grouped_duration, tuple)
            or len(self.grouped_duration) != 2
            or self.grouped_duration[1] != "day"
        ):
            errors.append(
                "`grouped_duration` must be a tuple with unit 'day', e.g. (5, 'day')."
            )

    # Raise errors if any
    if errors:
        raise ValueError("Parameter validation errors:\n\n" + "\n".join(errors))

`get_names()`

Generate names and metadata for climate data processing.

Returns:

Type Description

dict

A dictionary containing the following keys: ssp_selected : str or list The SSP (Shared Socioeconomic Pathway) data selected for analysis. variable : str The climate variable being analyzed, potentially adjusted based on downscaling method and metric calculation. variable_type : str Type of variable, either "Variable" for raw climate variables or "Derived Index" for calculated indices. var_name : str Human-readable description of the calculation being performed, including variable, metric, time period/warming level, and season. raw_name : str Standardized name for raw data storage. calc_name : str Standardized name for calculated data storage.

Notes

The function handles three main calculation types based on instance attributes: - Percentile calculations (when self.percentile is not None) - Heat index threshold calculations (when self.heat_idx_threshold is not None) - Return period calculations (when self.one_in_x is not None) For LOCA2 statistical downscaling with air temperature, the variable name is adjusted based on whether maximum or minimum metric calculation is selected. The approach can be either "Time" (using start/end years) or "Warming Level" (using a specific warming level in degrees Celsius).

Raises:

Type	Description
`ValueError`	If metric_calc is not "max" or "min" for LOCA2 air temperature data.
`ValueError`	If approach is not "Time" or "Warming Level".
`ValueError`	If an unsupported variable is used for 1-in-X calculations.

Source code in climakitae/explore/vulnerability.py

def get_names(self):
    """
    Generate names and metadata for climate data processing.

    Returns
    -------
    dict
        A dictionary containing the following keys:
        ssp_selected : str or list
            The SSP (Shared Socioeconomic Pathway) data selected for analysis.
        variable : str
            The climate variable being analyzed, potentially adjusted based on
            downscaling method and metric calculation.
        variable_type : str
            Type of variable, either "Variable" for raw climate variables or
            "Derived Index" for calculated indices.
        var_name : str
            Human-readable description of the calculation being performed,
            including variable, metric, time period/warming level, and season.
        raw_name : str
            Standardized name for raw data storage.
        calc_name : str
            Standardized name for calculated data storage.

    Notes
    -----
    The function handles three main calculation types based on instance attributes:
    - Percentile calculations (when self.percentile is not None)
    - Heat index threshold calculations (when self.heat_idx_threshold is not None)
    - Return period calculations (when self.one_in_x is not None)
    For LOCA2 statistical downscaling with air temperature, the variable name is
    adjusted based on whether maximum or minimum metric calculation is selected.
    The approach can be either "Time" (using start/end years) or "Warming Level"
    (using a specific warming level in degrees Celsius).

    Raises
    ------
    ValueError
        If metric_calc is not "max" or "min" for LOCA2 air temperature data.
    ValueError
        If approach is not "Time" or "Warming Level".
    ValueError
        If an unsupported variable is used for 1-in-X calculations.
    """

    # Re-assigning `variable` in case it needs to be re-instantiated later.
    variable = self.variable

    # Determine variable type
    if variable in ["Air Temperature at 2m", "Precipitation (total)"]:
        variable_type = "Variable"
    else:
        variable_type = "Derived Index"

    # Adjust variable name for LOCA2 data
    if (
        variable == "Air Temperature at 2m"
        and self.downscaling_method == "Statistical"
    ):
        match self.metric_calc:
            case "max":
                variable = "Maximum air temperature at 2m"
            case "min":
                variable = "Minimum air temperature at 2m"
            case _:
                raise ValueError('metric_calc needs to be either "max" or "min"')

    # Reducing potential decimal places for warming level
    self.warming_level = round(self.warming_level, 3)

    # Adding time period or warming level to filename
    match self.approach:
        case "Time":
            approach_str = f"{self.time_start_year}_to_{self.time_end_year}"
            approach_var_name = f"from {approach_str}"
        case "Warming Level":
            if float(self.warming_level).is_integer():
                wl_str = str(int(self.warming_level))
            else:
                wl_str = str(self.warming_level).replace(".", "pt")

            approach_str = f"{wl_str}degreeWL"
            approach_var_name = f"for Warming Level {self.warming_level}°C"
        case _:
            raise ValueError(
                'approach needs to be either "Time" or "Warming Level"'
            )
    # Retrieving the name of the calculation about to occur
    if self.percentile is not None:

        def ordinal(n):
            """Find ordinal name for a number (1 -> 1st, 2 -> 2nd, 3 -> 3rd)"""
            return f"{n}{('th' if 11 <= n % 100 <= 13 else {1: 'st', 2: 'nd', 3: 'rd'}.get(n % 10, 'th'))}"

        var_name = f"{ordinal(self.percentile)} percentile of Daily {self.metric_calc.capitalize()} of {variable} for {self.season} seasons {approach_var_name}"
        raw_name = "likely_seasonal_temperature_raw_data"
        calc_name = (
            f"likely_seasonal_{self.season}_{self.metric_calc}_{approach_str}"
        )

    elif self.heat_idx_threshold is not None:
        var_name = f"Days per year with {self.metric_calc.capitalize()} Daily {variable} above a heat index threshold of {self.heat_idx_threshold} {self.units} for {self.season} seasons {approach_var_name}"
        raw_name = "heat_index_raw_data"
        calc_name = f"heat_index_thresh{self.heat_idx_threshold}_{approach_str}"

    elif self.one_in_x is not None:

        event_parts = []
        if self.groupby is not UNSET:
            event_parts.append(f"groupby={self.groupby}")
        if self.duration is not UNSET:
            event_parts.append(f"duration={self.duration}")
        if self.grouped_duration is not UNSET:
            event_parts.append(f"grouped_duration={self.grouped_duration}")
        event_str = f" {', '.join(event_parts)}," if event_parts else ""

        x_vals = ", ".join([f"1-in-{x}" for x in self.one_in_x])
        var_name = f"{x_vals} year,{event_str} {self.metric_calc.capitalize()} {variable} for {self.season} seasons {approach_var_name}"

        # Making names for raw and calculated data
        x_vals_str = "-".join(map(str, self.one_in_x))

        match variable:
            case "Precipitation (total)":
                raw_name = f"one_in_{x_vals_str}_precipitation_raw_data"
                calc_name = f"one_in_{x_vals_str}_precipitation_{approach_str}"
            case (
                "Air Temperature at 2m"
                | "Maximum air temperature at 2m"
                | "Minimum air temperature at 2m"
            ):
                raw_name = f"one_in_{x_vals_str}_temperature_raw_data"
                calc_name = f"one_in_{x_vals_str}_temperature_{approach_str}"
            case "NOAA Heat Index" | "Effective Temperature":
                raw_name = f"one_in_{x_vals_str}_derived_index_raw_data"
                calc_name = f"one_in_{x_vals_str}_derived_index_{approach_str}"
            case _:
                raise ValueError(
                    f'variable "{variable}" is not supported for 1-in-X calculations'
                )

    return {
        "ssp_selected": self.ssp_data,
        "variable": variable,
        "variable_type": variable_type,
        "var_name": var_name,
        "raw_name": raw_name,
        "calc_name": calc_name,
    }

Threshold Tools

Helper functions for performing analyses related to thresholds

`calculate_ess(data, nlags=UNSET)`

Function for calculating the effective sample size (ESS) of the provided data using the autocorrelation of the data.

Parameters:

Name	Type	Description	Default
`data`	`DataArray`	Input array is assumed to be timeseries data with potential autocorrelation.	required
`nlags`	`int`	Number of lags to use in the autocorrelation function, defaults to the length of the timeseries.	`UNSET`

Returns:

Type	Description
`DataArray`	Effective sample size. Returned as a DataArray object so it can be utilized by xr.groupby and xr.resample.

References

Zwiers, F. W., and H. von Storch, 1995: Taking Serial Correlation into Account in Tests of the Mean. J. Climate, 8, 336–351, https://doi.org/10.1175/1520-0442(1995)008<0336:TSCIAI>2.0.CO;2.

Source code in climakitae/explore/threshold_tools.py

def calculate_ess(data: xr.DataArray, nlags: int = UNSET) -> xr.DataArray:
    """Function for calculating the effective sample size (ESS) of the provided data
    using the autocorrelation of the data.

    Parameters
    ----------
    data : xr.DataArray
        Input array is assumed to be timeseries data with potential autocorrelation.
    nlags : int, optional
        Number of lags to use in the autocorrelation function, defaults to the length of
        the timeseries.

    Returns
    -------
    xr.DataArray
        Effective sample size.
        Returned as a DataArray object so it can be utilized by xr.groupby and xr.resample.

    References
    ----------
    Zwiers, F. W., and H. von Storch, 1995: Taking Serial Correlation into Account in Tests of the Mean. J. Climate, 8, 336–351, https://doi.org/10.1175/1520-0442(1995)008<0336:TSCIAI>2.0.CO;2.

    """
    n = len(data)
    if nlags is UNSET:
        nlags = n
    acf = sm.tsa.stattools.acf(data, nlags=nlags, fft=True)
    sums = 0
    for k in range(1, len(acf)):
        sums = sums + (n - k) * acf[k] / n
    ess = n / (1 + 2 * sums)
    return xr.DataArray(ess, name="ess")

`get_block_maxima(da_series, extremes_type='max', duration=UNSET, groupby=UNSET, grouped_duration=UNSET, check_ess=True, block_size=1, rolling_agg='sustained')`

Convert a timeseries into block extrema (maxima or minima), defaulting to annual block size.

Resamples the input array by taking the most extreme value in each block. Optional arguments duration, groupby, and grouped_duration define the type of compound event to extract extrema from, corresponding to the event types used in get_exceedance_count.

Parameters:

Name	Type	Description	Default
`da_series`	`DataArray`	Input timeseries from which block extrema are extracted.	required
`extremes_type`	`(max, min)`	Whether to extract block maxima or minima. Default is `"max"`.	`"max"`
`duration`	`tuple[int, str]`	Length of the extreme event, specified as `(n, unit)` where unit is `"hour"`. For example, `(4, "hour")` finds events lasting at least 4 hours. Only supported for hourly data.	`UNSET`
`groupby`	`tuple[int, str]`	Temporal grouping applied before computing block extrema, specified as `(n, unit)` where unit is `"day"`. For example, `(1, "day")` aggregates to daily values. Most meaningful when combined with `grouped_duration`.	`UNSET`
`grouped_duration`	`tuple[int, str]`	Rolling window applied after `groupby`, specified as `(n, unit)` where unit is `"day"`. For example, `(5, "day")` finds events spanning at least 5 consecutive grouped periods. Requires `groupby`.	`UNSET`
`check_ess`	`bool`	If `True`, computes the effective sample size (ESS) within each block and emits a warning if the average ESS is below 25. Set to `False` to suppress this check. Default is `True`.	`True`
`block_size`	`int`	Block size in years over which to take the extremum. Default is `1`.	`1`
`rolling_agg`	`(sustained, cumulative, average)`	Aggregation method applied during rolling windows (`duration`, `grouped_duration`) and groupby resamples (`groupby`): - `"sustained"`: the extreme value that is maintained throughout the entire window (rolling min for max events, rolling max for min events). Use this when the event intensity must hold for the full duration. Example: the minimum temperature floor of a 21-day heatwave. - `"cumulative"`: the total accumulated value over the window (rolling sum). Use this for events defined by an accumulated total. Example: total precipitation of a 3-day heavy rainfall event. - `"average"`: the mean value over the window (rolling mean). Use this for events characterized by average intensity. Example: average precipitation rate of a short high-intensity burst. Default is `"sustained"`.	`"sustained"`

Returns:

Type	Description
`DataArray`	Block extrema series with the same non-time dimensions as `da_series`. The time coordinate reflects the end of each block (annual or `block_size`-year periods). Metadata describing the event configuration is stored in `.attrs`.

Raises:

Type	Description
`ValueError`	If `rolling_agg` is not one of `"sustained"`, `"cumulative"`, or `"average"`.
`ValueError`	If `extremes_type` is not `"max"` or `"min"`.
`ValueError`	If `duration` is provided but the data frequency is not hourly or the unit is not `"hour"`.
`ValueError`	If `groupby` is provided but the unit is not `"day"`.
`ValueError`	If `grouped_duration` is provided without `groupby`, or the unit is not `"day"`.
`ValueError`	If all block values are NaN (i.e. the input contains no valid observations for this variable).

Source code in climakitae/explore/threshold_tools.py

def get_block_maxima(
    da_series: xr.DataArray,
    extremes_type: str = "max",
    duration: tuple[int, str] = UNSET,
    groupby: tuple[int, str] = UNSET,
    grouped_duration: tuple[int, str] = UNSET,
    check_ess: bool = True,
    block_size: int = 1,
    rolling_agg: str = "sustained",
) -> xr.DataArray:
    """Convert a timeseries into block extrema (maxima or minima), defaulting to annual block size.

    Resamples the input array by taking the most extreme value in each block.
    Optional arguments ``duration``, ``groupby``, and ``grouped_duration`` define the
    type of compound event to extract extrema from, corresponding to the event types
    used in ``get_exceedance_count``.

    Parameters
    ----------
    da_series : xr.DataArray
        Input timeseries from which block extrema are extracted.
    extremes_type : {"max", "min"}, optional
        Whether to extract block maxima or minima. Default is ``"max"``.
    duration : tuple[int, str], optional
        Length of the extreme event, specified as ``(n, unit)`` where unit is
        ``"hour"``. For example, ``(4, "hour")`` finds events lasting at least
        4 hours. Only supported for hourly data.
    groupby : tuple[int, str], optional
        Temporal grouping applied before computing block extrema, specified as
        ``(n, unit)`` where unit is ``"day"``. For example, ``(1, "day")``
        aggregates to daily values. Most meaningful when combined with
        ``grouped_duration``.
    grouped_duration : tuple[int, str], optional
        Rolling window applied after ``groupby``, specified as ``(n, unit)``
        where unit is ``"day"``. For example, ``(5, "day")`` finds events
        spanning at least 5 consecutive grouped periods. Requires ``groupby``.
    check_ess : bool, optional
        If ``True``, computes the effective sample size (ESS) within each block
        and emits a warning if the average ESS is below 25. Set to ``False`` to
        suppress this check. Default is ``True``.
    block_size : int, optional
        Block size in years over which to take the extremum. Default is ``1``.
    rolling_agg : {"sustained", "cumulative", "average"}, optional
        Aggregation method applied during rolling windows (``duration``,
        ``grouped_duration``) and groupby resamples (``groupby``):
            - ``"sustained"``: the extreme value that is maintained throughout the
              entire window (rolling min for max events, rolling max for min events).
              Use this when the event intensity must hold for the full duration.
              *Example*: the minimum temperature floor of a 21-day heatwave.
            - ``"cumulative"``: the total accumulated value over the window (rolling
              sum). Use this for events defined by an accumulated total.
              *Example*: total precipitation of a 3-day heavy rainfall event.
            - ``"average"``: the mean value over the window (rolling mean). Use this
              for events characterized by average intensity.
              *Example*: average precipitation rate of a short high-intensity burst.
        Default is ``"sustained"``.

    Returns
    -------
    xr.DataArray
        Block extrema series with the same non-time dimensions as ``da_series``.
        The time coordinate reflects the end of each block (annual or
        ``block_size``-year periods). Metadata describing the event configuration
        is stored in ``.attrs``.

    Raises
    ------
    ValueError
        If ``rolling_agg`` is not one of ``"sustained"``, ``"cumulative"``, or
        ``"average"``.
    ValueError
        If ``extremes_type`` is not ``"max"`` or ``"min"``.
    ValueError
        If ``duration`` is provided but the data frequency is not hourly or the
        unit is not ``"hour"``.
    ValueError
        If ``groupby`` is provided but the unit is not ``"day"``.
    ValueError
        If ``grouped_duration`` is provided without ``groupby``, or the unit is
        not ``"day"``.
    ValueError
        If all block values are NaN (i.e. the input contains no valid
        observations for this variable).
    """
    valid_rolling_aggs = ["sustained", "cumulative", "average"]
    if rolling_agg not in valid_rolling_aggs:
        raise ValueError(
            f"invalid rolling_agg. expected one of the following: {valid_rolling_aggs}"
        )

    extremes_types = ["max", "min"]  # valid user options
    if extremes_type not in extremes_types:
        raise ValueError(
            "invalid extremes type. expected one of the following: %s" % extremes_types
        )

    if duration is not UNSET:
        # In this case, user is interested in extreme events lasting at least
        # as long as the length of `duration`.
        dur_len, dur_type = duration
        if dur_type != "hour" or da_series.frequency not in ["1hr", "hourly"]:
            raise ValueError(
                "Current specifications not implemented. `duration` options only implemented for `hour` frequency."
            )

        # First identify the min (max) value for each window of length `duration`
        if rolling_agg == "sustained":
            match extremes_type:
                case "max":
                    # In the case of "max" events, need to first identify the minimum value
                    # in each window of the specified duration
                    da_series = da_series.rolling(time=dur_len, center=False).min(
                        "time"
                    )
                case "min":
                    da_series = da_series.rolling(time=dur_len, center=False).max(
                        "time"
                    )
                case _:
                    raise ValueError('extremes_type needs to be either "max" or "min"')
        elif rolling_agg == "cumulative":
            da_series = da_series.rolling(time=dur_len, center=False).sum("time")
        elif rolling_agg == "average":
            da_series = da_series.rolling(time=dur_len, center=False).mean("time")

    if groupby is not UNSET:
        # In this case, select the max (min) in each group. (This option is
        # really only meaningful when coupled with the `grouped_duration` option.)
        group_len, group_type = groupby
        if group_type != "day":
            raise ValueError(
                "`groupby` specifications only implemented for 'day' groupings."
            )

        # select the max (min) in each group
        if rolling_agg == "sustained":
            match extremes_type:
                case "max":
                    da_series = da_series.resample(
                        time=f"{group_len}D", label="left"
                    ).max()
                case "min":
                    da_series = da_series.resample(
                        time=f"{group_len}D", label="left"
                    ).min()
                case _:
                    raise ValueError('extremes_type needs to be either "max" or "min"')
        elif rolling_agg == "cumulative":
            da_series = da_series.resample(time=f"{group_len}D", label="left").sum()
        elif rolling_agg == "average":
            da_series = da_series.resample(time=f"{group_len}D", label="left").mean()

    if grouped_duration is not UNSET:
        if groupby is UNSET:
            raise ValueError(
                "To use `grouped_duration` option, must first use groupby."
            )
        # In this case, identify the min (max) value of the grouped values for
        # each window of length `grouped_duration``. Must be in `days`.
        dur2_len, dur2_type = grouped_duration
        if dur2_type != "day":
            raise ValueError(
                "`grouped_duration` specification must be in days. example: `grouped_duration = (3, 'day')`."
            )

        # Rechunk time dimension if dask-backed to ensure chunks are large enough
        # for the rolling window. After groupby resampling, chunks may be smaller
        # than the rolling window size, causing a ValueError.
        if hasattr(da_series.data, "chunks"):
            da_series = da_series.chunk(time=-1)

        # Now select the min (max) from the duration period
        if rolling_agg == "sustained":
            match extremes_type:
                case "max":
                    da_series = da_series.rolling(time=dur2_len, center=False).min(
                        "time"
                    )
                case "min":
                    da_series = da_series.rolling(time=dur2_len, center=False).max(
                        "time"
                    )
                case _:
                    raise ValueError('extremes_type needs to be either "max" or "min"')
        elif rolling_agg == "cumulative":
            da_series = da_series.rolling(time=dur2_len, center=False).sum("time")
        elif rolling_agg == "average":
            da_series = da_series.rolling(time=dur2_len, center=False).mean("time")

    # Now select the most extreme value for each block in the series
    match extremes_type:
        case "max":
            bms = da_series.resample(time=f"{block_size}YE").max(keep_attrs=True)
            bms.attrs["extremes type"] = "maxima"
        case "min":
            bms = da_series.resample(time=f"{block_size}YE").min(keep_attrs=True)
            bms.attrs["extremes type"] = "minima"
        case _:
            raise ValueError('extremes_type needs to be either "max" or "min"')

    # Calculate the effective sample size of the computed event type in all blocks
    # Check the average value to ensure that it's above threshold ESS
    if check_ess:
        average_ess = -1
        if "x" in da_series.dims and "y" in da_series.dims:
            # Case for data with spatial dimensions (gridded)
            average_ess = _calc_average_ess_gridded_data(da_series, block_size)

        elif da_series.dims == ("time",):
            # Case for timeseries data (no spatial dimensions)
            average_ess = _calc_average_ess_timeseries_data(da_series, block_size)

        else:
            print(
                f"WARNING: the effective sample size can only be checked for timeseries or spatial data. You provided data with the following dimensions: {da_series.dims}."
            )

        if average_ess < 25 and average_ess > 0:
            print(
                f"WARNING: The average effective sample size in your data is {round(average_ess, 2)} per block, which is lower than a standard target of around 25. This may result in biased estimates of extreme value distributions when calculating return values, periods, and probabilities from this data. Consider using a longer block size to increase the effective sample size in each block of data."
            )

    # Common attributes
    bms = bms.assign_attrs(
        {
            "duration": duration,
            "groupby": groupby,
            "grouped_duration": grouped_duration,
            "rolling_agg": rolling_agg,
            "extreme_value_extraction_method": "block maxima",
            "block_size": f"{block_size} year",
            "timeseries_type": f"block {extremes_type} series",
        }
    )

    # Checking for null values within the blocks. This primarily occurs when values are filtered out before this function.
    # A good example of this is when trying to find 1-in-X events with precipitation. Because we filter out small noise of precip events,
    # (i.e. 10e-10), this results in a DataArray being passed into this function that has many fewer observed counts than the actual retrieve
    # climate data has. In some cases, there can even be years where no precipitation occurs. In these cases, we want to make sure to drop those
    # time blocks from the returned blocks below, hence dropping NaN values along the time dimension. This way, they do not break other functions
    # that rely on `get_block_maxima`, like `get_ks_stat` or `get_return_value`
    if bms.isnull().sum() > 0:

        # This checks if ALL of the values from `bms` are null, or if there are 0 events that occur (i.e. no precipitation counts within the DataArray).
        if bms.isnull().sum().item() == bms.size:
            raise ValueError(
                "ERROR: The given `da_series` does not include any recorded values for this variable, and we cannot create block maximums off of an empty DataArray. Please check your input data and filters to ensure that there are observed values for this variable in the given DataArray."
            )
        else:
            # Determine which non-time dims exist for detecting all-NaN time steps
            non_time_dims = [d for d in bms.dims if d != "time"]

            if non_time_dims:
                # Drop time steps where ALL values across non-time dims are NaN
                all_nan_times = bms.isnull().all(dim=non_time_dims)
            else:
                # Timeseries data (only time dim)
                all_nan_times = bms.isnull()

            n_dropped = int(all_nan_times.sum().item())
            if n_dropped > 0:
                bms = bms.where(~all_nan_times, drop=True)
                print(
                    f"Dropping {n_dropped} block maxima NaN time steps from"
                    f"{f' {bms.name}' if bms.name else ''} DataArray."
                )

    return bms

`get_ks_stat(bms, distr='gev', multiple_points=True)`

Function to perform kstest on input DataArray

Creates a dataset of ks test d-statistics and p-values from an inputed maximum series.

Parameters:

Name	Type	Description	Default
`bms`	`DataArray`	Block maximum series, can be output from the function get_block_maxima()	required
`distr`	`str`	name of distribution to use	`'gev'`
`multiple_points`	`boolean`	Whether or not the data contains multiple points (has x, y dimensions)	`True`

Returns:

Type	Description
`Dataset`

Source code in climakitae/explore/threshold_tools.py

def get_ks_stat(
    bms: xr.DataArray, distr: str = "gev", multiple_points: bool = True
) -> xr.Dataset:
    """Function to perform kstest on input DataArray

    Creates a dataset of ks test d-statistics and p-values from an inputed
    maximum series.

    Parameters
    ----------
    bms : xarray.DataArray
        Block maximum series, can be output from the function get_block_maxima()
    distr : str
        name of distribution to use
    multiple_points : boolean
        Whether or not the data contains multiple points (has x, y dimensions)

    Returns
    -------
    xarray.Dataset

    """
    distr_func = _get_distr_func(distr)
    bms_attributes = bms.attrs

    if multiple_points:
        bms = (
            bms.stack(allpoints=["y", "x"])
            .dropna(dim="allpoints")
            .squeeze()
            .groupby("allpoints")
        )

    def ks_stat(bms):
        parameters, fitted_distr = _get_fitted_distr(bms, distr, distr_func)

        match distr:
            case "gev":
                cdf = "genextreme"
                args = (parameters["c"], parameters["loc"], parameters["scale"])
            case "gumbel":
                cdf = "gumbel_r"
                args = (parameters["loc"], parameters["scale"])
            case "weibull":
                cdf = "weibull_min"
                args = (parameters["c"], parameters["loc"], parameters["scale"])
            case "pearson3":
                cdf = "pearson3"
                args = (parameters["skew"], parameters["loc"], parameters["scale"])
            case "genpareto":
                cdf = "genpareto"
                args = (parameters["c"], parameters["loc"], parameters["scale"])
            case "gamma":
                cdf = "gamma"
                args = (parameters["a"], parameters["loc"], parameters["scale"])
            case _:
                raise ValueError(
                    'invalid distribution type. expected one of the following: ["gev", "gumbel", "weibull", "pearson3", "genpareto", "gamma"]'
                )

        try:
            ks = stats.kstest(bms, cdf, args=args)
            d_statistic = ks[0]
            p_value = ks[1]
        except (ValueError, ZeroDivisionError):
            d_statistic = np.nan
            p_value = np.nan

        return d_statistic, p_value

    d_statistic, p_value = xr.apply_ufunc(
        ks_stat,
        bms,
        input_core_dims=[["time"]],
        exclude_dims=set(("time",)),
        vectorize=True,
        output_core_dims=[[], []],
    )

    d_statistic = d_statistic.rename("d_statistic")
    new_ds = d_statistic.to_dataset()
    new_ds["p_value"] = p_value

    if multiple_points:
        new_ds = new_ds.unstack("allpoints")

    new_ds["d_statistic"].attrs["stat test"] = "KS test"
    new_ds["p_value"].attrs["stat test"] = "KS test"
    new_ds.attrs = bms_attributes
    new_ds.attrs["distribution"] = "{}".format(str(distr))
    new_ds["p_value"].attrs["units"] = UNSET
    new_ds["d_statistic"].attrs["units"] = UNSET
    return new_ds

`get_return_value(bms, return_period=10, distr='gev', bootstrap_runs=100, conf_int_lower_bound=2.5, conf_int_upper_bound=97.5, multiple_points=True, extremes_type='max', dropna_time=False, dim_to_fit='time')`

Creates xarray Dataset with return values and confidence intervals from maximum series.

Parameters:

Name	Type	Description	Default
`bms`	`DataArray`	Block maximum series, can be output from the function get_block_maxima()	required
`return_period`	`float`	The recurrence interval (in years) for which to calculate the return value	`10`
`distr`	`str`	The type of extreme value distribution to fit	`'gev'`
`bootstrap_runs`	`int`	Number of bootstrap samples	`100`
`conf_int_lower_bound`	`float`	Confidence interval lower bound	`2.5`
`conf_int_upper_bound`	`float`	Confidence interval upper bound	`97.5`
`multiple_points`	`boolean`	Whether or not the data contains multiple points (has x, y dimensions)	`True`
`extremes_type`	`str`	Whether to compute max ('max') or min ('min') extremes, by default 'max'.	`'max'`
`dropna_time`	`bool`	Whether to drop NaNs along the time axis	`False`
`dim_to_fit`	`str`	Name of the dimension that the distribution is going to be fit on.	`'time'`

Returns:

Type	Description
`Dataset`	Dataset with return values and confidence intervals

Notes

This function calls _get_return_variable, which will attempt to use the block_size attribute from the BMS to get the block size. If the block_size attribute is not found, an default of 1 will be used. The block size is used to get the annual probability of the desired event when the block size is > 1 year.

Source code in climakitae/explore/threshold_tools.py

def get_return_value(
    bms: xr.DataArray,
    return_period: float = 10,
    distr: str = "gev",
    bootstrap_runs: int = 100,
    conf_int_lower_bound: float = 2.5,
    conf_int_upper_bound: float = 97.5,
    multiple_points: bool = True,
    extremes_type: str = "max",
    dropna_time: bool = False,
    dim_to_fit: str = "time",
) -> xr.Dataset:
    """Creates xarray Dataset with return values and confidence intervals from maximum series.

    Parameters
    ----------
    bms : xarray.DataArray
        Block maximum series, can be output from the function get_block_maxima()
    return_period : float
        The recurrence interval (in years) for which to calculate the return value
    distr : str
        The type of extreme value distribution to fit
    bootstrap_runs : int
        Number of bootstrap samples
    conf_int_lower_bound : float
        Confidence interval lower bound
    conf_int_upper_bound : float
        Confidence interval upper bound
    multiple_points : boolean
        Whether or not the data contains multiple points (has x, y dimensions)
    extremes_type : str, optional
        Whether to compute max ('max') or min ('min') extremes, by default 'max'.
    dropna_time: boolean
        Whether to drop NaNs along the time axis
    dim_to_fit : str
        Name of the dimension that the distribution is going to be fit on.

    Returns
    -------
    xarray.Dataset
        Dataset with return values and confidence intervals

    Notes
    -----
    This function calls _get_return_variable, which will attempt to use the block_size
    attribute from the BMS to get the block size. If the block_size attribute is not found,
    an default of 1 will be used. The block size is used to get the annual probability of
    the desired event when the block size is > 1 year.
    """
    return _get_return_variable(
        bms,
        "return_value",
        return_period,
        distr,
        bootstrap_runs,
        conf_int_lower_bound,
        conf_int_upper_bound,
        multiple_points,
        extremes_type,
        dropna_time,
        dim_to_fit,
    )

`get_return_prob(bms, threshold, distr='gev', bootstrap_runs=100, conf_int_lower_bound=2.5, conf_int_upper_bound=97.5, multiple_points=True, extremes_type='max', dropna_time=False, dim_to_fit='time')`

Creates xarray Dataset with return probabilities and confidence intervals from maximum series.

Parameters:

Name	Type	Description	Default
`bms`	`DataArray`	Block maximum series, can be output from the function get_block_maxima()	required
`threshold`	`float`	The threshold value for which to calculate the probability of exceedance	required
`distr`	`str`	The type of extreme value distribution to fit	`'gev'`
`bootstrap_runs`	`int`	Number of bootstrap samples	`100`
`conf_int_lower_bound`	`float`	Confidence interval lower bound	`2.5`
`conf_int_upper_bound`	`float`	Confidence interval upper bound	`97.5`
`multiple_points`	`boolean`	Whether or not the data contains multiple points (has x, y dimensions)	`True`
`extremes_type`	`str`	Whether to compute max ('max') or min ('min') extremes, by default 'max'.	`'max'`
`dropna_time`	`bool`	Whether to drop NaNs along the time axis	`False`
`dim_to_fit`	`str`	Name of the dimension that the distribution is going to be fit on.	`'time'`

Returns:

Type	Description
`Dataset`	Dataset with return probabilities and confidence intervals

Notes

This function calls _get_return_variable, which will attempt to use the block_size attribute from the BMS to get the block size. If the block_size attribute is not found, an default of 1 will be used. The block size is used to get the annual probability of the desired event when the block size is > 1 year.

Source code in climakitae/explore/threshold_tools.py

def get_return_prob(
    bms: xr.DataArray,
    threshold: float,
    distr: str = "gev",
    bootstrap_runs: int = 100,
    conf_int_lower_bound: float = 2.5,
    conf_int_upper_bound: float = 97.5,
    multiple_points: bool = True,
    extremes_type: str = "max",
    dropna_time: bool = False,
    dim_to_fit: str = "time",
) -> xr.Dataset:
    """Creates xarray Dataset with return probabilities and confidence intervals from maximum series.

    Parameters
    ----------
    bms : xarray.DataArray
        Block maximum series, can be output from the function get_block_maxima()
    threshold : float
        The threshold value for which to calculate the probability of exceedance
    distr : str
        The type of extreme value distribution to fit
    bootstrap_runs : int
        Number of bootstrap samples
    conf_int_lower_bound : float
        Confidence interval lower bound
    conf_int_upper_bound : float
        Confidence interval upper bound
    multiple_points : boolean
        Whether or not the data contains multiple points (has x, y dimensions)
    extremes_type : str, optional
        Whether to compute max ('max') or min ('min') extremes, by default 'max'.
    dropna_time: boolean
        Whether to drop NaNs along the time axis
    dim_to_fit : str
        Name of the dimension that the distribution is going to be fit on.

    Returns
    -------
    xarray.Dataset
        Dataset with return probabilities and confidence intervals

    Notes
    -----
    This function calls _get_return_variable, which will attempt to use the block_size
    attribute from the BMS to get the block size. If the block_size attribute is not found,
    an default of 1 will be used. The block size is used to get the annual probability of
    the desired event when the block size is > 1 year.

    """
    return _get_return_variable(
        bms,
        "return_prob",
        threshold,
        distr,
        bootstrap_runs,
        conf_int_lower_bound,
        conf_int_upper_bound,
        multiple_points,
        extremes_type,
        dropna_time=dropna_time,
        dim_to_fit=dim_to_fit,
    )

`get_return_period(bms, return_value, distr='gev', bootstrap_runs=100, conf_int_lower_bound=2.5, conf_int_upper_bound=97.5, multiple_points=True, extremes_type='max', dropna_time=False, dim_to_fit='time')`

Creates xarray Dataset with return periods and confidence intervals from maximum series.

Parameters:

Name	Type	Description	Default
`bms`	`DataArray`	Block maximum series, can be output from the function get_block_maxima()	required
`return_value`	`float`	The threshold value for which to calculate the return period of occurance	required
`distr`	`str`	The type of extreme value distribution to fit	`'gev'`
`bootstrap_runs`	`int`	Number of bootstrap samples	`100`
`conf_int_lower_bound`	`float`	Confidence interval lower bound	`2.5`
`conf_int_upper_bound`	`float`	Confidence interval upper bound	`97.5`
`multiple_points`	`boolean`	Whether or not the data contains multiple points (has x, y dimensions)	`True`
`extremes_type`	`str`	Whether to compute max ('max') or min ('min') extremes, by default 'max'.	`'max'`
`dropna_time`	`bool`	Whether to drop NaNs along the time axis	`False`
`dim_to_fit`	`str`	Name of the dimension that the distribution is going to be fit on.	`'time'`

Returns:

Type	Description
`Dataset`	Dataset with return periods and confidence intervals

Notes

This function calls _get_return_variable, which will attempt to use the block_size attribute from the BMS to get the block size. If the block_size attribute is not found, an default of 1 will be used. The block size is used to get the annual probability of the desired event when the block size is > 1 year.

Source code in climakitae/explore/threshold_tools.py

def get_return_period(
    bms: xr.DataArray,
    return_value: float,
    distr: str = "gev",
    bootstrap_runs: int = 100,
    conf_int_lower_bound: float = 2.5,
    conf_int_upper_bound: float = 97.5,
    multiple_points: bool = True,
    extremes_type: str = "max",
    dropna_time: bool = False,
    dim_to_fit: str = "time",
) -> xr.Dataset:
    """Creates xarray Dataset with return periods and confidence intervals from maximum series.

    Parameters
    ----------
    bms : xarray.DataArray
        Block maximum series, can be output from the function get_block_maxima()
    return_value : float
        The threshold value for which to calculate the return period of occurance
    distr : str
        The type of extreme value distribution to fit
    bootstrap_runs : int
        Number of bootstrap samples
    conf_int_lower_bound : float
        Confidence interval lower bound
    conf_int_upper_bound : float
        Confidence interval upper bound
    multiple_points : boolean
        Whether or not the data contains multiple points (has x, y dimensions)
    extremes_type : str, optional
        Whether to compute max ('max') or min ('min') extremes, by default 'max'.
    dropna_time: boolean
        Whether to drop NaNs along the time axis
    dim_to_fit : str
        Name of the dimension that the distribution is going to be fit on.

    Returns
    -------
    xarray.Dataset
        Dataset with return periods and confidence intervals

    Notes
    -----
    This function calls _get_return_variable, which will attempt to use the block_size
    attribute from the BMS to get the block size. If the block_size attribute is not found,
    an default of 1 will be used. The block size is used to get the annual probability of
    the desired event when the block size is > 1 year.

    """
    return _get_return_variable(
        bms,
        "return_period",
        return_value,
        distr,
        bootstrap_runs,
        conf_int_lower_bound,
        conf_int_upper_bound,
        multiple_points,
        extremes_type,
        dropna_time=dropna_time,
        dim_to_fit=dim_to_fit,
    )

`get_exceedance_count(da, threshold_value, duration1=UNSET, period=(1, 'year'), threshold_direction='above', duration2=UNSET, groupby=UNSET, smoothing=UNSET)`

Calculate the number of occurances of exceeding the specified threshold within each period.

Returns an xarray.DataArray with the same coordinates as the input data except for the time dimension, which will be collapsed to one value per period (equal to the number of event occurances in each period).

Parameters:

Name	Type	Description	Default
`da`	`DataArray`	array of some climate variable. Can have multiple scenarios, simulations, or x and y coordinates.	required
`threshold_value`	`float`	value against which to test exceedance	required
`duration1`	`tuple[int, str]`	length of exceedance in order to qualify as an event (before grouping)	`UNSET`
`period`	`tuple[int, str]`	amount of time across which to sum the number of occurances, default is (1, "year"). Specified as a tuple: (x, time) where x is an integer, and time is one of: ["day", "month", "year"]	`(1, 'year')`
`threshold_direction`	`str`	either "above" or "below", default is above.	`'above'`
`duration2`	`tuple[int, str]`	length of exceedance in order to qualify as an event (after grouping)	`UNSET`
`groupby`	`tuple[int, str]`	see examples for explanation. Typical grouping could be (1, "day")	`UNSET`
`smoothing`	`int`	option to average the result across multiple periods with a rolling average; value is either UNSET or the number of timesteps to use as the window size	`UNSET`

Returns:

Type	Description
`DataArray`

Source code in climakitae/explore/threshold_tools.py

def get_exceedance_count(
    da: xr.DataArray,
    threshold_value: float,
    duration1: tuple[int, str] = UNSET,
    period: tuple[int, str] = (1, "year"),
    threshold_direction="above",
    duration2: tuple[int, str] = UNSET,
    groupby: tuple[int, str] = UNSET,
    smoothing: int = UNSET,
) -> xr.DataArray:
    """Calculate the number of occurances of exceeding the specified threshold
    within each period.

    Returns an xarray.DataArray with the same coordinates as the input data except for
    the time dimension, which will be collapsed to one value per period (equal
    to the number of event occurances in each period).

    Parameters
    ----------
    da : xarray.DataArray
        array of some climate variable. Can have multiple
        scenarios, simulations, or x and y coordinates.
    threshold_value : float
        value against which to test exceedance
    duration1 : tuple[int, str]
        length of exceedance in order to qualify as an event (before grouping)
    period : tuple[int, str]
        amount of time across which to sum the number of occurances,
        default is (1, "year"). Specified as a tuple: (x, time) where x is an
        integer, and time is one of: ["day", "month", "year"]
    threshold_direction : str
        either "above" or "below", default is above.
    duration2 : tuple[int, str]
        length of exceedance in order to qualify as an event (after grouping)
    groupby : tuple[int, str]
        see examples for explanation. Typical grouping could be (1, "day")
    smoothing : int
        option to average the result across multiple periods with a
        rolling average; value is either UNSET or the number of timesteps to use
        as the window size

    Returns
    -------
    xarray.DataArray

    """
    # --------- Type check arguments -------------------------------------------

    # Check compatibility of periods, durations, and groupbys
    if _is_greater(duration1, groupby):
        raise ValueError(
            "Incompatible `group` and `duration1` specification. Duration1 must be shorter than group."
        )
    if _is_greater(groupby, duration2):
        raise ValueError(
            "Incompatible `group` and `duration2` specification. Duration2 must be longer than group."
        )
    if _is_greater(groupby, period):
        raise ValueError(
            "Incompatible `group` and `period` specification. Group must be longer than period."
        )
    if _is_greater(duration2, period):
        raise ValueError(
            "Incompatible `duration` and `period` specification. Period must be longer than duration."
        )

    # Check compatibility of specifications with the data frequency (hourly, daily, or monthly)
    freq = (
        (1, "hour")
        if da.frequency == "hourly"
        else ((1, "day") if da.frequency == "daily" else (1, "month"))
    )
    if _is_greater(freq, groupby):
        raise ValueError(
            "Incompatible `group` specification: cannot be less than data frequency."
        )
    if _is_greater(freq, duration2):
        raise ValueError(
            "Incompatible `duration` specification: cannot be less than data frequency."
        )
    if _is_greater(freq, period):
        raise ValueError(
            "Incompatible `period` specification: cannot be less than data frequency."
        )

    # --------- Calculate occurances -------------------------------------------

    events_da = _get_exceedance_events(
        da, threshold_value, threshold_direction, duration1, groupby
    )

    # --------- Apply specified duration requirement ---------------------------

    if duration2 is not UNSET:
        dur_len, dur_type = duration2

        if (
            groupby is not UNSET
            and groupby[1] == dur_type
            or groupby is UNSET
            and freq[1] == dur_type
        ):
            window_size = dur_len
        else:
            raise ValueError(
                "Duration options for time types (i.e. hour, day) that are different than group or frequency not yet implemented"
            )

        # The "min" operation will return 0 if any time in the window is not an
        # event, which is the behavior we want. It will only return 1 for True
        # if all values in the duration window are 1.
        events_da = events_da.rolling(time=window_size, center=False).min("time")

    # --------- Sum occurances across each period ------------------------------

    period_len, period_type = period
    period_indexer = _get_freq_string(period_type)
    exceedance_count = events_da.resample(
        time=f"{period_len}{period_indexer}", label="left"
    ).sum()

    # Optional smoothing
    if smoothing is not UNSET:
        exceedance_count = exceedance_count.rolling(time=smoothing, center=True).mean(
            "time"
        )

    # --------- Set new attributes for the counts DataArray --------------------
    exceedance_count.attrs["variable_name"] = da.name
    exceedance_count.attrs["variable_units"] = exceedance_count.units
    exceedance_count.attrs["period"] = period
    exceedance_count.attrs["duration1"] = duration1
    exceedance_count.attrs["group"] = groupby
    exceedance_count.attrs["duration2"] = duration2
    exceedance_count.attrs["threshold_value"] = threshold_value
    exceedance_count.attrs["threshold_direction"] = threshold_direction
    exceedance_count.attrs["units"] = _exceedance_count_name(exceedance_count)

    # Set name (for plotting, this will be the y-axis label)
    exceedance_count.name = "Count"

    return exceedance_count

`exceedance_plot_title(exceedance_count)`

Function to build title for exceedance plots

Helper function for making the title for exceedance plots.

Parameters:

Name	Type	Description	Default
`exceedance_count`	`DataArray`		required

Returns:

Type	Description
`string`

Examples:

'Air Temperatue at 2m: events above 35C'
'Preciptation (total): events below 10mm'

Source code in climakitae/explore/threshold_tools.py

def exceedance_plot_title(exceedance_count: xr.DataArray) -> str:
    """Function to build title for exceedance plots

    Helper function for making the title for exceedance plots.

    Parameters
    ----------
    exceedance_count : xarray.DataArray

    Returns
    -------
    string

    Examples
    --------
        'Air Temperatue at 2m: events above 35C'
        'Preciptation (total): events below 10mm'

    """
    return f"{exceedance_count.variable_name}: events {exceedance_count.threshold_direction} {exceedance_count.threshold_value}{exceedance_count.variable_units}"

`exceedance_plot_subtitle(exceedance_count)`

Function of build exceedance plot subtitle

Helper function for making the subtile for exceedance plots.

Parameters:

Name	Type	Description	Default
`exceedance_count`	`DataArray`		required

Returns:

Type	Description
`string`

Examples:

'Number of hours per year'
'Number of 4-hour events per 3-months'
'Number of days per year with conditions lasting at least 4-hours'

Source code in climakitae/explore/threshold_tools.py

def exceedance_plot_subtitle(exceedance_count: xr.DataArray) -> str:
    """Function of build exceedance plot subtitle

    Helper function for making the subtile for exceedance plots.

    Parameters
    ----------
    exceedance_count : xarray.DataArray

    Returns
    -------
    string

    Examples
    --------
        'Number of hours per year'
        'Number of 4-hour events per 3-months'
        'Number of days per year with conditions lasting at least 4-hours'

    """
    if exceedance_count.duration2 != exceedance_count.duration1:
        dur_len, dur_type = exceedance_count.duration1
        _s = "" if dur_len == 1 else "s"
        dur_str = f" with conditions lasting at least {dur_len} {dur_type}{_s}"
    else:
        dur_str = ""

    if exceedance_count.duration2 != exceedance_count.group:
        grp_len, grp_type = exceedance_count.group
        if grp_len == 1:
            grp_str = f" each {grp_type}"
        else:
            grp_str = f" every {grp_len} {grp_type}s"
    else:
        grp_str = ""

    per_len, per_type = exceedance_count.period
    if per_len == 1:
        period_str = f" each {per_type}"
    else:
        period_str = f" per {per_len}-{per_type} period"

    _subtitle = (
        _exceedance_count_name(exceedance_count) + period_str + dur_str + grp_str
    )
    return _subtitle

Thresholds

`get_threshold_data(selections)`

This function pulls data from the catalog and reads it into memory

Arguments

selections : DataParameters object holding user's selections

Returns:

Name	Type	Description
`data`	`DataArray`	data to use for creating postage stamp data

Source code in climakitae/explore/thresholds.py

def get_threshold_data(selections: DataParameters) -> xr.DataArray:
    """This function pulls data from the catalog and reads it into memory

    Arguments
    ---------
    selections : DataParameters
        object holding user's selections

    Returns
    -------
    data : xr.DataArray
        data to use for creating postage stamp data

    """
    # Read data from catalog
    data = read_catalog_from_select(selections)
    data = data.compute()  # Read into memory
    return data

Vulnerability Tables

`create_vul_table(example_loc, percentile, heat_idx_threshold, one_in_x)`

Creates a vulnerability assessment table and exports the table to CSV.

Source code in climakitae/explore/vulnerability_table.py

def create_vul_table(example_loc, percentile, heat_idx_threshold, one_in_x):
    """Creates a vulnerability assessment table and exports the table to CSV."""
    # Create empty df and instantiate variables
    df = pd.DataFrame(columns=table_vars.keys())
    lat, lon = example_loc.lat.item(), example_loc.lon.item()
    months_map = {"winter": [12, 1, 2], "summer": [6, 7, 8], "all": np.arange(1, 13)}
    time_periods = [
        (1981, 2010),
        (2021, 2050),
        (2031, 2060),
        (2061, 2090),
    ]

    # Create each column in the table, which is the historical period (1980-2010) and each WL (1.5, 2.0, 3.0, 4.0).
    for time_period in time_periods:

        metrics = []

        # Calculate each variable in the table
        for key in table_vars:
            params = table_vars[key]

            # Suppress outputs of `cava_data` function
            if suppress_output:

                with contextlib.redirect_stdout(io.StringIO()):

                    data = cava_data(
                        example_loc,
                        variable=params["variable"],
                        approach="time",
                        downscaling_method="Dynamical",
                        time_start_year=time_period[0],
                        time_end_year=time_period[1],
                        historical_data="Historical Climate",  # or "historical reconstruction"
                        ssp_data=["SSP3-7.0"],
                        metric_calc=params["metric_calc"],
                        heat_idx_threshold=params["heat_idx_threshold"]
                        and heat_idx_threshold,  # Heat index
                        one_in_x=params["one_in_x"]
                        and one_in_x,  # Thresholds tools freq. counts
                        percentile=params["percentile"] and percentile,
                        season=params["season"],
                        units="degF",
                        wrf_bias_adjust=True,
                        export_method="None",  # off-ramp, full calculate, both
                        separate_files=True,  # Toggle to determine whether or not the user wants to separate climate variable information into separate files
                        file_format="NetCDF",
                    )

            # Retrieve data and average across simulation dimension
            val = data["calc_data"][0].mean(dim="simulation").item()

            # Add val to metrics to be added into row
            metrics.append(val)

        # Create dictionary of values to be input into DataFrame
        df.loc[str(time_period)] = pd.Series(dict(zip(table_vars.keys(), metrics)))

    # Make slight modifications to DataFrame
    # df = df.T.rename(columns={"0.0": "Hist. Period (1980-2010)"})
    df = df.T

    # Write out dataframe
    df.to_csv(f"final_table_{example_loc.index.item()}.csv", index=True)

    return df