Skip to content

Tools Module

Computational tools for climate indices and batch processing.

Overview

The climakitae.tools module provides functions for: - Batch processing — Process multiple locations efficiently - Climate indices — Calculate heat stress, effective temperature, and other climate metrics - Derived variables — Compute secondary climate variables from primary ones

Batch Processing

batch_select(approach, selections, points, load_data=False, progress_bar=True)

Conducts batch mode analysis on a series of points for a given metric.

Parameters:

Name Type Description Default
approach str
required
selections DataParameters

Selections object that describes the area of interest. The area_subset and cached_area attributes are automatically overwritten.

 required
points ndarray

An array at lat/lon points to gather the specified data at.

 required
load_data boolean

A boolean that tells the function whether or not to load the data into memory.

 False
progress_bar boolean

A boolean that determines whether progress bar is displayed.

 True

Returns:

Name Type Description
cells_of_interest DataArray

Gridcells that the points lie within, aggregated together into one DataArray. It can or cannot be loaded into memory, depending on load_data.
Source code in climakitae/tools/batch.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
def batch_select(
    approach: str,
    selections: DataParameters,
    points: np.ndarray,
    load_data: bool = False,
    progress_bar: bool = True,
) -> xr.DataArray:
    """Conducts batch mode analysis on a series of points for a given metric.

    Parameters
    ----------
    approach : str
    selections : DataParameters
        Selections object that describes the area of interest. The `area_subset` and `cached_area` attributes are automatically overwritten.
    points : np.ndarray
        An array at lat/lon points to gather the specified data at.
    load_data : boolean
        A boolean that tells the function whether or not to load the data into memory.
    progress_bar : boolean
        A boolean that determines whether progress bar is displayed.

    Returns
    -------
    cells_of_interest: xr.DataArray
        Gridcells that the points lie within, aggregated together into one DataArray. It can or cannot be loaded into memory, depending on `load_data`.

    """
    print(f"Batch retrieving all {len(points)} points passed in...\n")

    # Add selections attributes to cover the entire domain since we don't know exactly where the selected points lie.
    selections.area_subset = "none"
    selections.cached_area = ["entire domain"]

    data = selections.retrieve()

    if approach == "Time":
        # Remove leap days, if applicable
        data = data.sel(time=~((data.time.dt.month == 2) & (data.time.dt.day == 29)))

    # Find the closest gridcells for each of the passed in points and concatenate them on a new 'points' dimension to go from 2D grid to 1D series of points
    cells_of_interest = get_closest_gridcells(data, points[:, 0], points[:, 1])

    # Load in the cells of interest into memory, if desired.
    if load_data:
        cells_of_interest = load(cells_of_interest, progress_bar=progress_bar)

    return cells_of_interest

Climate Indices

Functions for deriving indices

effective_temp(T)

Compute effective temperature Effective Temp = (1/2)(yesterday's effective temp) + (1/2)(today's actual temp) To make sense of the expansion, today's ET is consist of a portion of the actual temperature of each day up to today--half of today's temp, 1/4 of yesterday's temp, 1/8 of the day before yesterday's temp etc, thus it's "an exponentially smoothed temperature" as stated in the glossary of the reference. This derivation only considers 4 days of temperature data in the computation of EFT: today, yesterday, the day before yesterday, and two days before yesterday. Thus, the first 3 timesteps of the EFT will be NaN.

Parameters:

Name Type Description Default
T DataArray

Daily air temperature in any units

 required

Returns:

Name Type Description
eft DataArray

Effective temperature
References

https://www.nationalgas.com/document/132516/download

Source code in climakitae/tools/indices.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
def effective_temp(T: xr.DataArray) -> xr.DataArray:
    """Compute effective temperature
    Effective Temp = (1/2)*(yesterday's effective temp) + (1/2)*(today's actual temp)
    To make sense of the expansion, today's ET is consist of a portion of the actual temperature of each day up to today--half of today's temp, 1/4 of yesterday's temp, 1/8 of the day before yesterday's temp etc, thus it's "an exponentially smoothed temperature" as stated in the glossary of the reference.
    This derivation only considers 4 days of temperature data in the computation of EFT: today, yesterday, the day before yesterday, and two days before yesterday. Thus, the first 3 timesteps of the EFT will be NaN.

    Parameters
    ----------
    T : xr.DataArray
        Daily air temperature in any units

    Returns
    -------
    eft : xr.DataArray
        Effective temperature

    References
    ----------
    https://www.nationalgas.com/document/132516/download

    """
    # Get "yesterday" temp by shifting the time index back one time step (1 day)
    # Get "day before" temp by shifting the time index back two time steps (2 days)
    # Get "2 days before yesterday" temp by shifting the time index back three time steps (3 days)
    T_minus1 = T.shift(time=1)
    T_minus2 = T.shift(time=2)
    T_minus3 = T.shift(time=3)

    # Compute EFT, using 3 days back
    # Effective temp for 2 days before yesterday is set to the temperature of that day
    eft_minus3 = T_minus3
    eft_minus2 = eft_minus3 * 0.5 + T_minus2 * 0.5
    eft_minus1 = eft_minus2 * 0.5 + T_minus1 * 0.5
    eft = eft_minus1 * 0.5 + T * 0.5

    # Assign same attributes as input data
    # Or else, the output data will have no attributes :(
    eft.attrs = T.attrs

    return eft

noaa_heat_index(T, RH)

Compute the NOAA Heat Index. Heat Index quantifies the perceived "real feel" of air temperature on the human body, including the impact of humidity. See references for more information on this derived variable.

Parameters:

Name Type Description Default
T DataArray

Temperature in deg F

 required
RH DataArray

Relative Humidity in percentage (0-100)

 required

Returns:

Name Type Description
HI DataArray

Heat index per timestep
References

NOAA: https://www.wpc.ncep.noaa.gov/html/heatindex_equation.shtml NCAR NCL documentation: https://www.ncl.ucar.edu/Document/Functions/Heat_stress/heat_index_nws.shtml

Source code in climakitae/tools/indices.py 
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
def noaa_heat_index(T: xr.DataArray, RH: xr.DataArray) -> xr.DataArray:
    """Compute the NOAA Heat Index.
    Heat Index quantifies the perceived "real feel" of air temperature on the human body,
    including the impact of humidity. See references for more information on this derived variable.

    Parameters
    ----------
    T : xr.DataArray
        Temperature in deg F
    RH : xr.DataArray
        Relative Humidity in percentage (0-100)

    Returns
    -------
    HI : xr.DataArray
        Heat index per timestep

    References
    ----------
    NOAA: https://www.wpc.ncep.noaa.gov/html/heatindex_equation.shtml
    NCAR NCL documentation: https://www.ncl.ucar.edu/Document/Functions/Heat_stress/heat_index_nws.shtml

    """
    T = T.reindex_like(RH)  # Need to have the same dimension/coordinate orders
    HI = (
        -42.379
        + 2.04901523 * T
        + 10.14333127 * RH
        - 0.22475541 * T * RH
        - 0.00683783 * T * T
        - 0.05481717 * RH * RH
        + 0.00122874 * T * T * RH
        + 0.00085282 * T * RH * RH
        - 0.00000199 * T * T * RH * RH
    )

    # Adjust for high temperature, low relative humidity
    # 80 < T < 112 (deg F)
    # RH < 13%
    adj_highT_lowRH = ((13 - RH) / 4) * ((17 - abs(T - 95)) / 17) ** (
        1 / 2
    )  # Adjustment
    HI_highT_lowRH = HI - adj_highT_lowRH  # Subtract adjustment from HI

    # Adjust for low temperature, high relative humidity
    # 80 < T < 87 (deg F)
    # RH > 85%
    adj_lowT_highRH = ((RH - 85) / 10) * ((87 - T) / 5)  # Adjustment
    HI_lowT_highRH = HI + adj_lowT_highRH  # Add adjustment from HI

    # Use different equation if heat index if the heat index value < 80
    low_HI = 0.5 * (T + 61.0 + ((T - 68.0) * 1.2) + (RH * 0.094))

    # Adjust heat index depending on different conditions for RH, T, and valid range of HI
    HI = xr.where((RH < 13) & (T > 80) & (T < 112), HI_highT_lowRH, HI)
    HI = xr.where(((RH > 85) & (T < 87) & (T > 80)), HI_lowT_highRH, HI)
    # Per NOAA: use the simple formula first; only apply Rothfusz when simple >= 80
    HI = xr.where((low_HI < 80), low_HI, HI)

    # Following NCAR documentation (see function references), for temperature values less than 40F, the HI is set to the ambient temperature.
    HI = xr.where((T < 40), T, HI)

    # Reassign coordinate attributes
    # For some reason, these get improperly assigned in the xr.where step
    for coord in list(HI.coords):
        HI[coord].attrs = RH[coord].attrs

    # Assign units attribute
    HI.attrs["units"] = "degF"
    return HI

fosberg_fire_index(t2_F, rh_percent, windspeed_mph)

Compute the Fosberg Fire Weather Index. Use hourly weather as inputs. Ensure that the input variables are in the correct units (see below).

Parameters:

Name Type Description Default
t2_F DataArray

Air temperature in units of Fahrenheit

 required
rh_percent DataArray

Relative humidity in units of 0-100 (percent)

 required
windspeed_mph DataArray

Windspeed in units of miles per hour

 required

Returns:

Name Type Description
FFWI DataArray

Fosberg Fire Weather Index computed for each grid cell
References

https://a.atmos.washington.edu/wrfrt/descript/definitions/fosbergindex.html https://github.com/sharppy/SHARPpy/blob/main/sharppy/sharptab/fire.py https://www.spc.noaa.gov/exper/firecomp/INFO/fosbinfo.html

Source code in climakitae/tools/indices.py 
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
def fosberg_fire_index(
    t2_F: xr.DataArray, rh_percent: xr.DataArray, windspeed_mph: xr.DataArray
) -> xr.DataArray:
    """Compute the Fosberg Fire Weather Index.
    Use hourly weather as inputs.
    Ensure that the input variables are in the correct units (see below).

    Parameters
    ----------
    t2_F : xr.DataArray
        Air temperature in units of Fahrenheit
    rh_percent : xr.DataArray
        Relative humidity in units of 0-100 (percent)
    windspeed_mph : xr.DataArray
        Windspeed in units of miles per hour

    Returns
    -------
    FFWI: xr.DataArray
        Fosberg Fire Weather Index computed for each grid cell

    References
    ----------
    https://a.atmos.washington.edu/wrfrt/descript/definitions/fosbergindex.html
    https://github.com/sharppy/SHARPpy/blob/main/sharppy/sharptab/fire.py
    https://www.spc.noaa.gov/exper/firecomp/INFO/fosbinfo.html

    """
    # Compute the equilibrium moisture constant
    m_low, m_mid, m_high = _equilibrium_moisture_constant(h=rh_percent, T=t2_F)

    # For RH < 10%, use the low m value.
    # For RH >= 10%, use the mid value
    m = xr.where(rh_percent < 10, m_low, m_mid)
    # For RH > 50%, use the high m value.
    m = xr.where(rh_percent > 50, m_high, m)

    # Compute the moisture dampening coefficient
    n = _moisture_dampening_coeff(m)

    # Compute the index
    U = windspeed_mph
    # If the value falls out of [0-100] range clip the value
    tmp = (n * ((1 + U**2) ** 0.5)) / 0.3002
    FFWI = tmp.clip(min=0.0, max=100.0)

    # Reassign coordinate attributes
    # For some reason, these get improperly assigned in the xr.where step
    for coord in list(FFWI.coords):
        FFWI[coord].attrs = t2_F[coord].attrs

    # Add descriptive attributes
    FFWI.name = "Fosberg Fire Weather Index"
    FFWI.attrs["units"] = "[0 to 100]"

    return FFWI

Derived Variables

Functions for deriving frequently used variables

compute_hdd_cdd(t2, hdd_threshold, cdd_threshold)

Compute heating degree days (HDD) and cooling degree days (CDD)

Parameters:

Name Type Description Default
t2 DataArray

Air temperature at 2m gridded data

 required
hdd_threshold int

Standard temperature in Fahrenheit.

 required
cdd_threshold int

Standard temperature in Fahrenheit.

 required

Returns:

Type Description
tuple of xr.DataArray

(hdd, cdd)
Source code in climakitae/tools/derived_variables.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def compute_hdd_cdd(
    t2: xr.DataArray, hdd_threshold: int, cdd_threshold: int
) -> tuple[xr.DataArray, xr.DataArray]:
    """Compute heating degree days (HDD) and cooling degree days (CDD)

    Parameters
    ----------
    t2 : xr.DataArray
        Air temperature at 2m gridded data
    hdd_threshold : int, optional
        Standard temperature in Fahrenheit.
    cdd_threshold : int, optional
        Standard temperature in Fahrenheit.

    Returns
    -------
    tuple of xr.DataArray
        (hdd, cdd)

    """

    # Check that temperature data was passed to function, throw error if not
    if t2.name not in ["t2", "Air Temperature at 2m"]:
        raise Exception(
            "Invalid input data, please provide Air Temperature at 2m data to CDD/HDD calculation"
        )

    # Subtract t2 from the threshold inputs
    hdd_deg_less_than_standard = hdd_threshold - t2
    cdd_deg_less_than_standard = cdd_threshold - t2

    # Compute HDD: Find positive difference (i.e. days < 65 degF)
    hdd = hdd_deg_less_than_standard.clip(0, None)
    # Replace negative values with 0
    hdd.name = "Heating Degree Days"
    hdd.attrs["hdd_threshold"] = (
        str(hdd_threshold) + " degF"
    )  # add attribute of threshold value

    # Compute CDD: Find negative difference (i.e. days > 65 degF)
    cdd = (-1) * cdd_deg_less_than_standard.clip(None, 0)
    # Replace positive values with 0
    cdd.name = "Cooling Degree Days"
    cdd.attrs["cdd_threshold"] = (
        str(cdd_threshold) + " degF"
    )  # add attribute of threshold value

    return (hdd, cdd)

compute_hdh_cdh(t2, hdh_threshold, cdh_threshold)

Compute heating degree hours (HDH) and cooling degree hours (CDH)

Parameters:

Name Type Description Default
t2 DataArray

Air temperature at 2m gridded data

 required
hdh_threshold int

Standard temperature in Fahrenheit.

 required
cdh_threshold int

Standard temperature in Fahrenheit.

 required

Returns:

Type Description
tuple of xr.DataArray

(hdh, cdh)
Source code in climakitae/tools/derived_variables.py 
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
def compute_hdh_cdh(
    t2: xr.DataArray, hdh_threshold: int, cdh_threshold: int
) -> tuple[xr.DataArray, xr.DataArray]:
    """Compute heating degree hours (HDH) and cooling degree hours (CDH)

    Parameters
    ----------
    t2 : xr.DataArray
        Air temperature at 2m gridded data
    hdh_threshold : int, optional
        Standard temperature in Fahrenheit.
    cdh_threshold : int, optional
        Standard temperature in Fahrenheit.

    Returns
    -------
    tuple of xr.DataArray
        (hdh, cdh)

    """

    # Check that temperature data was passed to function, throw error if not
    if t2.name not in ["t2", "Air Temperature at 2m"]:
        raise Exception(
            "Invalid input data, please provide Air Temperature at 2m data to CDH/HDH calculation"
        )

    # Calculate heating and cooling hours
    cooling_hours = t2.where(
        t2 > cdh_threshold
    )  # temperatures above threshold, require cooling
    heating_hours = t2.where(
        t2 < hdh_threshold
    )  # temperatures below threshold, require heating

    # Compute CDH: count number of hours and resample to daily (max 24 value)
    cdh = cooling_hours.resample(time="1D").count(dim="time").squeeze()
    cdh.name = "Cooling Degree Hours"
    cdh.attrs["cdh_threshold"] = str(cdh_threshold) + " degF"

    # Compute HDH: count number of hours and resample to daily (max 24 value)
    hdh = heating_hours.resample(time="1D").count(dim="time").squeeze()
    hdh.name = "Heating Degree Hours"
    hdh.attrs["hdh_threshold"] = str(hdh_threshold) + " degF"

    return (hdh, cdh)

compute_dewpointtemp(temperature, rel_hum)

Calculate dew point temperature

Parameters:

Name Type Description Default
temperature DataArray

Temperature in Kelvin (K)

 required
rel_hum DataArray

Relative humidity (0-100 scale)

 required

Returns:

Name Type Description
dew_point DataArray

Dew point (K)
Source code in climakitae/tools/derived_variables.py 
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
def compute_dewpointtemp(
    temperature: xr.DataArray, rel_hum: xr.DataArray
) -> xr.DataArray:
    """Calculate dew point temperature

    Parameters
    ----------
    temperature : xr.DataArray
        Temperature in Kelvin (K)
    rel_hum : xr.DataArray
        Relative humidity (0-100 scale)

    Returns
    -------
    dew_point : xr.DataArray
        Dew point (K)

    """
    es = 0.611 * np.exp(
        5423 * ((1 / 273) - (1 / temperature))
    )  # calculates saturation vapor pressure
    e_vap = (es * rel_hum) / 100.0  # calculates vapor pressure
    tdps = (
        (1 / 273) - 0.0001844 * np.log(e_vap / 0.611)
    ) ** -1  # calculates dew point temperature, units = K

    # Assign descriptive name
    tdps.name = "dew_point_derived"
    tdps.attrs["units"] = "K"
    return tdps

compute_specific_humidity(tdps, pressure, name='q2_derived')

Compute specific humidity.

Parameters:

Name Type Description Default
tdps DataArray

Dew-point temperature, in K

 required
pressure DataArray

Air pressure, in Pascals

 required
name str

Name to assign to output DataArray

 'q2_derived'

Returns:

Name Type Description
spec_hum DataArray

Specific humidity
Source code in climakitae/tools/derived_variables.py 
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
def compute_specific_humidity(
    tdps: xr.DataArray, pressure: xr.DataArray, name: str = "q2_derived"
) -> xr.DataArray:
    """Compute specific humidity.

    Parameters
    ----------
    tdps : xr.DataArray
        Dew-point temperature, in K
    pressure : xr.DataArray
        Air pressure, in Pascals
    name : str, optional
        Name to assign to output DataArray

    Returns
    -------
    spec_hum : xr.DataArray
        Specific humidity

    """

    # Calculate vapor pressure, unit is in kPa
    e = 0.611 * np.exp((2500000 / 461) * ((1 / 273) - (1 / tdps)))

    # Calculate specific humidity, unit is g/g, pressure has to be divided by 1000 to get to kPa at this step
    q = (0.622 * e) / (pressure / 1000)

    # Convert from g/g to g/kg for more understandable value
    q = q * 1000

    # Assign descriptive name
    q.name = name
    q.attrs["units"] = "g/kg"
    return q

compute_relative_humidity(pressure, temperature, mixing_ratio, name='rh_derived')

Compute relative humidity. Variable attributes need to be assigned outside of this function because the metpy function removes them

Parameters:

Name Type Description Default
pressure DataArray

Pressure in hPa

 required
temperature DataArray

Temperature in Celsius

 required
mixing_ratio DataArray

Dimensionless mass mixing ratio in g/kg

 required
name str

Name to assign to output DataArray

 'rh_derived'

Returns:

Name Type Description
rel_hum DataArray

Relative humidity
Source https://www.weather.gov/media/epz/wxcalc/mixingRatio.pdf
Source code in climakitae/tools/derived_variables.py 
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
def compute_relative_humidity(
    pressure: xr.DataArray,
    temperature: xr.DataArray,
    mixing_ratio: xr.DataArray,
    name: str = "rh_derived",
) -> xr.DataArray:
    """Compute relative humidity.
    Variable attributes need to be assigned outside of this function because the metpy function removes them

    Parameters
    ----------
    pressure : xr.DataArray
        Pressure in hPa
    temperature : xr.DataArray
        Temperature in Celsius
    mixing_ratio : xr.DataArray
        Dimensionless mass mixing ratio in g/kg
    name : str, optional
        Name to assign to output DataArray

    Returns
    -------
    rel_hum : xr.DataArray
        Relative humidity

    Source: https://www.weather.gov/media/epz/wxcalc/mixingRatio.pdf

    """

    # Calculates saturated vapor pressure
    e_s = 6.11 * 10 ** (7.5 * (temperature / (237.7 + temperature)))

    # calculate saturation mixing ratio, unit is g/kg
    w_s = 621.97 * (e_s / (pressure - e_s))

    # Calculates relative humidity, unit is 0 to 100
    rel_hum = 100 * (mixing_ratio / w_s)

    # Reset unrealistically low relative humidity values
    # Lowest recorded relative humidity value in CA is 0.8%
    rel_hum = xr.where(rel_hum > 0.5, rel_hum, 0.5)

    # Reset values above 100 to 100
    rel_hum = xr.where(rel_hum < 100, rel_hum, 100)

    # Reassign coordinate attributes
    # For some reason, these get improperly assigned in the xr.where step
    for coord in list(rel_hum.coords):
        rel_hum[coord].attrs = temperature[coord].attrs

    # Assign descriptive name
    rel_hum.name = name
    rel_hum.attrs["units"] = "[0 to 100]"
    return rel_hum

compute_wind_mag(u10, v10, name='wind_speed_derived')

Compute wind magnitude at 10 meters

Parameters:

Name Type Description Default
u10 DataArray

Zonal velocity at 10 meters height in m/s

 required
v10 DataArray

Meridonal velocity at 10 meters height in m/s

 required
name str

Name to assign to output DataArray

 'wind_speed_derived'

Returns:

Name Type Description
wind_mag DataArray

Wind magnitude
Source code in climakitae/tools/derived_variables.py 
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
def compute_wind_mag(
    u10: xr.DataArray, v10: xr.DataArray, name: str = "wind_speed_derived"
) -> xr.DataArray:
    """Compute wind magnitude at 10 meters

    Parameters
    ----------
    u10 : xr.DataArray
        Zonal velocity at 10 meters height in m/s
    v10 : xr.DataArray
        Meridonal velocity at 10 meters height in m/s
    name : str, optional
        Name to assign to output DataArray

    Returns
    -------
    wind_mag : xr.DataArray
        Wind magnitude

    """
    wind_mag = np.sqrt(np.square(u10) + np.square(v10))
    wind_mag.name = name
    wind_mag.attrs["units"] = "m s-1"
    return wind_mag

compute_wind_dir(u10, v10, name='wind_direction_derived')

Compute wind direction at 10 meters

Parameters:

Name Type Description Default
u10 DataArray

Zonal velocity at 10 meters height in m/s

 required
v10 DataArray

Meridional velocity at 10 meters height in m/s

 required
name str

Name to assign to output DataArray

 'wind_direction_derived'

Returns:

Name Type Description
wind_dir DataArray

Wind direction, in [0, 360] degrees, with 0/360 defined as north, by meteorological convention
Notes 
source: https://sites.google.com/view/raybellwaves/cheat-sheets/xarray
Source code in climakitae/tools/derived_variables.py 
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
def compute_wind_dir(
    u10: xr.DataArray, v10: xr.DataArray, name: str = "wind_direction_derived"
) -> xr.DataArray:
    """Compute wind direction at 10 meters

    Parameters
    ----------
    u10 : xr.DataArray
        Zonal velocity at 10 meters height in m/s
    v10 : xr.DataArray
        Meridional velocity at 10 meters height in m/s
    name : str, optional
        Name to assign to output DataArray

    Returns
    -------
    wind_dir : xr.DataArray
        Wind direction, in [0, 360] degrees, with 0/360 defined as north, by meteorological convention

    Notes
    -----
        source: https://sites.google.com/view/raybellwaves/cheat-sheets/xarray

    """

    wind_dir = np.mod(90 - np.arctan2(-v10, -u10) * (180 / np.pi), 360)
    wind_dir.name = name
    wind_dir.attrs["units"] = "degrees"
    return wind_dir

compute_sea_level_pressure(psfc, t2, q2, elevation, lapse_rate=0.0065, average_t2=True, name='slp_derived')

Calculate sea level pressure from hourly surface pressure, temperature, and mixing ratio.

This function uses the basic method derived from the hydrostatic balance equation and the equation of state (Hess 1979). The SLP calculation method used here may not produce satisfactory results in locations with high terrain.

By default this method uses a standard lapse rate of 6.5°K/km when calculating the sea level virtual temperature (see Pauley 1998). Users should consider what lapse rate is appropriate for their location.

An option is provided to use a 12-hour average temperature when computing the lapse rate; this option is expected to produce more moderate SLP values that are less influenced by extreme temperatures.

Parameters:

Name Type Description Default
psfc DataArray

Hourly surface pressure in Pascals

 required
t2 DataArray

Hourly surface air temperature in Kelvin

 required
q2 DataArray

Hourly surface mixing ratio

 required
elevation DataArray

Elevation in meters

 required
lapse_rate Union[float, DataArray]

Lapse rate in K/m. Default is 0.0065 K/m

 0.0065
average_t2 bool

True to use 12-hour mean temperature. Default is True.

 True
name str

Name to assign to output DataArray

 'slp_derived'

Returns:

Type Description
DataArray

Sea level pressure in Pascals
Notes

Virtual temperature is computed in the following way: T_virtual = ((1 + 1.609 q2) / (1 + q2)) * t2 T_virtual_mean = (2 * T_virtual + lapse_rate * elevation) / 2

Sea level pressure is calculated as: slp = psfc * np.exp(elevation / ((Rd * T_virtual_mean)/g)) where Rd is the specific gas constant for dry air and g is the acceleration due to gravity.

References

Hess, S. L., 1979: Introduction to Theoretical Meteorology. Robert E. Krieger Publishing Company, 362 pp. Pauley, P. M., 1998: An Example of Uncertainty in Sea Level Pressure Reduction. Wea. Forecasting, 13, 833–850, https://doi.org/10.1175/1520-0434(1998)013<0833:AEOUIS>2.0.CO;2.

Source code in climakitae/tools/derived_variables.py 
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
def compute_sea_level_pressure(
    psfc: xr.DataArray,
    t2: xr.DataArray,
    q2: xr.DataArray,
    elevation: xr.DataArray,
    lapse_rate: Union[float, xr.DataArray] = 0.0065,
    average_t2: bool = True,
    name: str = "slp_derived",
) -> xr.DataArray:
    """Calculate sea level pressure from hourly surface pressure, temperature, and mixing ratio.

    This function uses the basic method derived from the hydrostatic balance equation
    and the equation of state (Hess 1979). The SLP calculation method used here may not produce
    satisfactory results in locations with high terrain.

    By default this method uses a standard lapse rate of 6.5°K/km when calculating the
    sea level virtual temperature (see Pauley 1998). Users should consider what lapse rate is
    appropriate for their location.

    An option is provided to use a 12-hour average temperature when computing the lapse rate;
    this option is expected to produce more moderate SLP values that are less influenced by
    extreme temperatures.

    Parameters
    ----------
    psfc : xr.DataArray
        Hourly surface pressure in Pascals
    t2 : xr.DataArray
        Hourly surface air temperature in Kelvin
    q2 : xr.DataArray
        Hourly surface mixing ratio
    elevation : xr.DataArray
        Elevation in meters
    lapse_rate : Union[float, xr.DataArray]
        Lapse rate in K/m. Default is 0.0065 K/m
    average_t2 : bool, optional
        True to use 12-hour mean temperature. Default is True.
    name : str, optional
        Name to assign to output DataArray

    Returns
    -------
    xr.DataArray
        Sea level pressure in Pascals

    Notes
    -----
    Virtual temperature is computed in the following way:
    T_virtual = ((1 + 1.609 q2) / (1 + q2)) * t2
    T_virtual_mean = (2 * T_virtual + lapse_rate * elevation) / 2

    Sea level pressure is calculated as:
    slp = psfc * np.exp(elevation / ((Rd * T_virtual_mean)/g))
       where Rd is the specific gas constant for dry air
       and g is the acceleration due to gravity.

    References
    ----------
    Hess, S. L., 1979: Introduction to Theoretical Meteorology. Robert E. Krieger Publishing Company, 362 pp.
    Pauley, P. M., 1998: An Example of Uncertainty in Sea Level Pressure Reduction. Wea. Forecasting, 13, 833–850, https://doi.org/10.1175/1520-0434(1998)013<0833:AEOUIS>2.0.CO;2.
    """
    # Get mean virtual temperature
    if average_t2:
        logger.info("compute_sea_level_pressure: Using 12-timestep mean temperature.")
        if "time" in t2.dims:
            t2 = t2.rolling(time=12).mean()
        elif "time_delta" in t2.dims:
            t2 = t2.rolling(time_delta=12).mean()
        else:
            raise KeyError(
                "No time or time_delta axis found in t2. Use `average_t2=False` for data without time axis."
            )

    t_virtual_sfc = ((1 + 1.609 * q2) / (1 + q2)) * t2
    t_virtual_mean = (2 * t_virtual_sfc + lapse_rate * elevation) / 2

    # Adjust pressure with hypsometric equation
    Rd = 287.04  # gas constant for dry air, J kg−1 K−1
    g = 9.81  # acceleration due to gravity, m s-2

    h = (Rd * t_virtual_mean) / g
    slp = psfc * np.exp(elevation / h)
    slp.name = name
    slp.attrs["units"] = "Pa"
    return slp

compute_geostrophic_wind(geopotential_height, gridlabel='d01')

Calculate the geostrophic wind at a single point on a constant pressure surface.

Currently only implemented for data on the WRF grid. This code follows the MetPy code for calculating the geostrophic wind on an unevenly spaced grid.

Parameters:

Name Type Description Default
geopotential_height DataArray

Geopotential height in meters on WRF grid. May include multiple pressure levels

 required

Returns:

Type Description
tuple[DataArray]

Earth-relative U and V components of the geostrophic wind.
References

Hess, S. L., 1979: Introduction to Theoretical Meteorology. Robert E. Krieger Publishing Company, 362 pp. MetPy, 2026: geostrophic_wind. Accessed 10 Feb 2026, https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.geostrophic_wind.html#geostrophic-wind.

Source code in climakitae/tools/derived_variables.py 
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
def compute_geostrophic_wind(
    geopotential_height: xr.DataArray, gridlabel="d01"
) -> tuple[xr.DataArray]:
    """Calculate the geostrophic wind at a single point on a constant pressure surface.

    Currently only implemented for data on the WRF grid. This code follows the
    MetPy code for calculating the geostrophic wind on an unevenly spaced grid.

    Parameters
    ----------
    geopotential_height : xr.DataArray
        Geopotential height in meters on WRF grid. May include multiple pressure levels

    Returns
    -------
    tuple[xr.DataArray]
        Earth-relative U and V components of the geostrophic wind.

    References
    ----------
    Hess, S. L., 1979: Introduction to Theoretical Meteorology. Robert E. Krieger Publishing Company, 362 pp.
    MetPy, 2026: geostrophic_wind. Accessed 10 Feb 2026, https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.geostrophic_wind.html#geostrophic-wind.
    """
    lat_to_radian = geopotential_height.lat.data * np.pi / 180
    omega = 7292115e-11  # rad/s
    g = 9.81  # m/s2
    f = 2 * omega * np.sin(lat_to_radian)
    norm_factor = g / f

    dhdx, dhdy = _get_spatial_derivatives(geopotential_height)

    # These components are u and v on the WRF grid
    geo_u, geo_v = -norm_factor * dhdy, norm_factor * dhdx

    # Rotate these components to an earth-relative E/W orientation
    geo_u_earth, geo_v_earth = _get_rotated_geostrophic_wind(geo_u, geo_v, gridlabel)

    # Update attributes for results
    geo_u_earth.name = "u"
    geo_u_earth.attrs["long_name"] = "Geostrophic Wind U Component"
    geo_v_earth.name = "v"
    geo_v_earth.attrs["long_name"] = "Geostrophic Wind V Component"

    return geo_u_earth, geo_v_earth