The core module

The core module.

Contains fundamental utilities and physics functionality which is shared across the code base. Submodules are the utilities and the hygro submodule.

The bounds submodule

Some functions in pyrealm are only well-behaved with given bounds but those bounds are often a little imprecise and real world data can contain extreme values. As a result, the bounds checking is deliberately not that intrusive: it warns when a variable contains out of value issues but leaves it up to the user to assess whether there is real problem and to adjust input data if needed.

The bounds module:

  • Defines a Bounds dataclass used to define bounds for a particular variable.

  • Defines a BoundsChecker class with default bounds for core variables that acts as a library for bounds checking.

  • The main use case is e.g. BoundsChecker().check("tc", np.array([10, 1000]), which will check that the alleged temperature data in °C fall within the configured bounds.

A BoundsChecker class instance is created with a predefined internal dictionary of default variables and appropriate bounds. However, users can use the update() method to override defaults or add new variables by providing a new Bounds instance.

The check() method can then be used to validate a set of values against the configured bounds for a given variable name. The check method returns the input variables, to allow values to be checked while being assigned to an attribute.

Classes:

Bounds(var_name, lower, upper, ...)

Bounds checking dataclass for variables.

BoundsChecker(*args, **kwargs)

A bounds checker for input variables.

class Bounds(var_name: str, lower: float, upper: float, interval_type: str, unit: str)

Bounds checking dataclass for variables.

Attributes:

interval_type

The interval type of the constraint ('[]', '()', '[)', '(]').

lower

A lower bound on sensible values.

unit

A string giving the expected units.

upper

An upper bound on sensible values.

var_name

A variable name, typically the form used in function arguments.

interval_type: str

The interval type of the constraint (‘[]’, ‘()’, ‘[)’, ‘(]’).

lower: float

A lower bound on sensible values.

unit: str

A string giving the expected units.

upper: float

An upper bound on sensible values.

var_name: str

A variable name, typically the form used in function arguments.

class BoundsChecker(*args: Any, **kwargs: Any)

A bounds checker for input variables.

The class provides a library of Bounds instances for core variables, keyed by the Bounds.var_name attribute. The table is populated from default values when a BoundsChecker instance is created but can be updated and extended by assigning new Bounds instances to existing or new variable name keys using the update method.

Methods:

check(var_name, values)

Check inputs fall within bounds.

update(bounds)

Update or add bounds data.

check(var_name: str, values: ndarray[tuple[Any, ...], dtype[floating]] | DataArray) ndarray[tuple[Any, ...], dtype[floating]]

Check inputs fall within bounds.

This method checks whether the provided values fall within the bounds specified for the given variable name and issues a warning when this is not the case. If the BoundsChecker class has not been configured the variable name then a warning will be given about lack of bounds checking. The method returns the input values, so that the method can be used as a pass through validator for assigning attributes.

Parameters:
  • var_name – The variable name

  • values – An np.ndarray

Returns:

The input values.

Examples

>>> vals = np.array([-15, 20, 30, 124], dtype=float)
>>> bounds_checker = BoundsChecker()
>>> bounds_checker.check("temp", vals)
array([-15.,  20.,  30., 124.])
update(bounds: Bounds) None

Update or add bounds data.

The Bounds.var_name attribute of the provided Bounds instance is used to update an existing entry for the name or add checking for a new name.

Parameters:

bounds – A Bounds instance.

The calendar submodule

Provides the Calendar and CalendarDay utility classes.

The pyrealm.core.calendar.Calendar class is currently used to support the operation of the splash submodule, which requires the Julian day, number of days in the year and year for solar calculations. The class provides iterable and indexable access to a date sequence for use in those calculations, returning individual pyrealm.core.calendar.CalendarDay instances.

It is possible that this could be replaced with xarray dt accessors if pyrealm adopts xarray data structures.

Classes:

Calendar(dates)

The Calendar class.

CalendarDay(date, year, julian_day, days_in_year)

The CalendarDay class.

class Calendar(dates: ndarray[tuple[Any, ...], dtype[datetime64]] | DataArray)

The Calendar class.

This utility class takes a numpy array of numpy.datetime64 values representing a time series of individual days and calculates the date, year, julian day and days in the year for each observation. The class object can be iterated over, yielding each date in turn, and dates within the series can also be accessed by index. In both cases, the returned object is a CalendarDay instance.

Examples

>>> days=np.arange(
...     np.datetime64("2000-01-01"),
...     np.datetime64("2002-01-01"),
...     np.timedelta64(1, "D")
... )
>>> cal = Calendar(days)
>>> cal
Calendar(2000-01-01, 2001-12-31)
>>> len(cal) == (366 + 365)
True
>>> cal[0]
CalendarDay(date=2000-01-01, year=2000, julian_day=1, days_in_year=366)
>>> for date in cal:
...     pass
>>> date
CalendarDay(date=2001-12-31, year=2001, julian_day=365, days_in_year=365)

Attributes:

dates

A numpy array containing numpy.datetime64 values.

days_in_year

A numpy array giving the number of days in the year for each datetime.

julian_day

A numpy array giving the julian day of each datetime.

year

A numpy array giving the year of each datetime.

dates: ndarray[tuple[Any, ...], dtype[datetime64]]

A numpy array containing numpy.datetime64 values.

days_in_year: ndarray[tuple[Any, ...], dtype[int64]]

A numpy array giving the number of days in the year for each datetime.

julian_day: ndarray[tuple[Any, ...], dtype[int64]]

A numpy array giving the julian day of each datetime.

year: ndarray[tuple[Any, ...], dtype[int64]]

A numpy array giving the year of each datetime.

class CalendarDay(date: datetime64, year: int, julian_day: int, days_in_year: int)

The CalendarDay class.

This dataclass holds a numpy.datetime64 datetime representing a day and the corresponding year, julian day (the integer index of the day within the year) and the total number of days in the year.

Attributes:

date

The date of the instance as numpy.datetime64.

days_in_year

The total number of days in the year.

julian_day

The julian day of the instance as an integer.

year

The year of the instance as an integer.

date: datetime64

The date of the instance as numpy.datetime64.

days_in_year: int

The total number of days in the year.

julian_day: int

The julian day of the instance as an integer.

year: int

The year of the instance as an integer.

The datasets submodule

The dataset module provides a download process for datasets from the pyrealm_build_data package. This package is only installed as part of the sdist package, so the files are accessible using importlib.resources for source installations. However, we do not want to ship the pyrealm_build_data package as part of the package binary because some of the files are quite large. This module provides:

  • The get_pyrealm_data() function to access these datasets, either by using importlib.resources for use in development environments where pyrealm_build_data is present, or by using the pooch package to download requested files to a local cache.

  • The DATASETS object, which is a pooch.Pooch instance that manages the available datasets.

  • The pyrealm_data_registry.json file that provides a dictionary mapping the available datasets by relative path to their SHA256 hashes.

  • The private _populate_pooch_registry() function to generate pyrealm_data_registry.json from a local copy of pyrealm_build_data.

The module contains an if __name__ == "__main__": section that allows the _populate_pooch_registry() function to be run using python -m pyrealm.core.datasets.

Functions:

_populate_pooch_registry([dataset_filetypes])

Generate a listing file of available pyrealm_build_data datasets.

get_pyrealm_data(filepath[, use_resources])

Get a path to pyrealm_build_data datasets, possibly downloading the dataset.

_populate_pooch_registry(dataset_filetypes: tuple[str, ...] = ('.csv', '.nc', '.json')) None

Generate a listing file of available pyrealm_build_data datasets.

This function recursively searches the pyrealm_build_data package for files with suffixes matching the dataset_filetypes argument. It then writes the pyrealm_data_registry.json file into the pyrealm.core module, which contains a dictionary keying the file path of each dataset relative to the pyrealm_build_data root to the the SHA256 hash of the file.

The command can be run manually using:

python -m  pyrealm.core.datasets
Parameters:

dataset_filetypes – A tuple of file suffixes for the types of file to be included in the registry.

get_pyrealm_data(filepath: str, use_resources: bool = False) str

Get a path to pyrealm_build_data datasets, possibly downloading the dataset.

This function returns a path to a dataset provided in the pyrealm_build_data package. These datasets are not provided in the binary build of pyrealm so are only locally available from sdist installations or from within a git clone of the package repository. A list of available files is maintained in the pyrealm_data_registry.json file within the pyrealm.core module and managed through the DATASETS file manager.

The function runs in one of two ways:

  • If use_resources=True is used or the PYREALM_USE_LOCAL_DATA environment variable is set, then the function using importlib.resources.files() to provide a path to the local copy.

  • Otherwise, the function uses the pooch.Pooch.fetch() method on the DATASETS instance to download the requested file to a local cache and then returns the path to that file.

The pyrealm package sphinx configuration for building documentation sets PYREALM_USE_LOCAL_DATA to ensure documentation builds use local data.

Parameters:
  • filepath – A string giving the path to the required dataset, relative to the root of the pyrealm_build_data package.

  • use_resources – A boolean toggle to switch to providing a path to an existing local installation of the pyrealm_build_data package, rather than downloading the requested file.

DATASETS = <pooch.core.Pooch object>

A pooch.Pooch instance used to manage the available dataset files.

The hygro submodule

The hygro submodule provides conversion functions for common hygrometric variables. The module provides conversions to vapour pressure deficit, which is the required input for the PModel from vapour pressure, specific humidity and relative humidity. The implementation is drawn from and validated against the bigleaf R package.

Functions:

calculate_enthalpy_vaporisation(tc)

Calculate the enthalpy of vaporization.

calculate_psychrometric_constant(tc, p[, ...])

Calculate the psychrometric constant.

calculate_saturation_vapour_pressure_slope(tc)

Calculate the slope of the saturation vapour pressure curve.

calculate_specific_heat(tc)

Calculate the specific heat of air.

calculate_vp_sat(ta[, core_const])

Calculate vapour pressure of saturated air.

convert_rh_to_vpd(rh, ta[, core_const, ...])

Convert relative humidity to vapour pressure deficit.

convert_sh_to_vp(sh, patm[, core_const])

Convert specific humidity to vapour pressure.

convert_sh_to_vpd(sh, ta, patm[, core_const])

Convert specific humidity to vapour pressure deficit.

convert_vp_to_vpd(vp, ta[, core_const])

Convert vapour pressure to vapour pressure deficit.

calculate_enthalpy_vaporisation(tc: ndarray[tuple[Any, ...], dtype[floating]] | DataArray) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the enthalpy of vaporization.

Calculates the latent heat of vaporization of water as a function of temperature following Henderson-Sellers (1984).

Parameters:

tc – Air temperature (°C)

Returns:

Calculated latent heat of vaporisation (J/Kg).

calculate_psychrometric_constant(tc: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, p: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the psychrometric constant.

Calculates the psychrometric constant (\(\lambda\), Pa/K) given the temperature and atmospheric pressure following Allen et al. (1998) and Tsilingiris (2008).

Parameters:
  • tc – Air temperature (°C)

  • p – Atmospheric pressure (Pa)

  • core_const – An instance of CoreConst giving the settings to be used in conversions.

Returns:

The calculated psychrometric constant

calculate_saturation_vapour_pressure_slope(tc: ndarray[tuple[Any, ...], dtype[floating]] | DataArray) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the slope of the saturation vapour pressure curve.

Calculates the slope of the saturation pressure temperature curve, following equation 13 of Allen et al. (1998).

Parameters:

tc – The air temperature (°C)

Returns:

The calculated slope in kPa °C-1.

calculate_specific_heat(tc: ndarray[tuple[Any, ...], dtype[floating]] | DataArray) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the specific heat of air.

Calculates the specific heat of air at a constant pressure (\(c_{pm}\), J/kg/K) following Tsilingiris (2008). This equation is only valid for temperatures between 0 and 100 °C.

Parameters:

tc – Air temperature (°C)

Returns:

The specific heat of air values.

calculate_vp_sat(ta: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate vapour pressure of saturated air.

This function calculates the vapour pressure of saturated air in kPa at a given temperature in °C, using the Magnus equation:

\[P = a \exp\(\frac{b - T}{T + c}\),\]

where \(a,b,c\) are defined in magnus_coef.

Parameters:
  • ta – The air temperature in °C.

  • core_const – An instance of CoreConst giving the parameters for conversions.

Returns:

Saturated air vapour pressure in kPa.

Examples

>>> # Saturated vapour pressure at 21°C
>>> import numpy as np
>>> temp = np.array([21])
>>> calculate_vp_sat(temp).round(6)
array([2.480904])
>>> from pyrealm.constants import CoreConst
>>> allen = CoreConst(magnus_option='Allen1998')
>>> calculate_vp_sat(temp, core_const=allen).round(6)
array([2.487005])
>>> alduchov = CoreConst(magnus_option='Alduchov1996')
>>> calculate_vp_sat(temp, core_const=alduchov).round(6)
array([2.481888])
convert_rh_to_vpd(rh: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, ta: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, core_const: CoreConst = CoreConst(), bounds_checker: BoundsChecker = BoundsChecker()) ndarray[tuple[Any, ...], dtype[floating]]

Convert relative humidity to vapour pressure deficit.

Parameters:
  • rh – The relative humidity (proportion in (0,1))

  • ta – The air temperature in °C

  • core_const – An instance of CoreConst giving the settings to be used in conversions.

  • bounds_checker – A BoundsChecker instance used to validate inputs.

Returns:

The vapour pressure deficit in kPa

Examples

>>> import numpy as np
>>> from pyrealm.constants import CoreConst
>>> import sys; sys.stderr = sys.stdout
>>> rh = np.array([0.7])
>>> temp = np.array([21])
>>> convert_rh_to_vpd(rh, temp).round(7)
array([0.7442712])
>>> allen = CoreConst(magnus_option='Allen1998')
>>> convert_rh_to_vpd(rh, temp, core_const=allen).round(7)
array([0.7461016])
convert_sh_to_vp(sh: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, patm: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Convert specific humidity to vapour pressure.

Parameters:
  • sh – The specific humidity in kg kg-1

  • patm – The atmospheric pressure in kPa

  • core_const – An instance of CoreConst giving the settings to be used in conversions.

Returns:

The vapour pressure in kPa

Examples

>>> import numpy as np
>>> sh = np.array([0.006])
>>> patm = np.array([99.024])
>>> convert_sh_to_vp(sh, patm).round(7)
array([0.9517451])
convert_sh_to_vpd(sh: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, ta: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, patm: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Convert specific humidity to vapour pressure deficit.

Parameters:
  • sh – The specific humidity in kg kg-1

  • ta – The air temperature in °C

  • patm – The atmospheric pressure in kPa

  • core_const – An instance of CoreConst giving the settings to be used in conversions.

Returns:

The vapour pressure deficit in kPa

Examples

>>> import numpy as np
>>> from pyrealm.constants import CoreConst
>>> sh = np.array([0.006])
>>> temp = np.array([21])
>>> patm = np.array([99.024])
>>> convert_sh_to_vpd(sh, temp, patm).round(6)
array([1.529159])
>>> allen = CoreConst(magnus_option='Allen1998')
>>> convert_sh_to_vpd(sh, temp, patm, core_const=allen).round(5)
array([1.53526])
convert_vp_to_vpd(vp: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, ta: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Convert vapour pressure to vapour pressure deficit.

Parameters:
  • vp – The vapour pressure in kPa

  • ta – The air temperature in °C

  • core_const – An instance of CoreConst giving the settings to be used in conversions.

Returns:

The vapour pressure deficit in kPa

Examples

>>> import numpy as np
>>> from pyrealm.constants import CoreConst
>>> vp = np.array([1.9])
>>> temp = np.array([21])
>>> convert_vp_to_vpd(vp, temp).round(7)
array([0.5809042])
>>> allen = CoreConst(magnus_option='Allen1998')
>>> convert_vp_to_vpd(vp, temp, core_const=allen).round(7)
array([0.5870054])

The pressure submodule

The pressure submodule contains core functions for calculating atmospheric pressure.

Functions:

calculate_patm(elv[, core_const])

Calculate atmospheric pressure from elevation.

calculate_patm(elv: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate atmospheric pressure from elevation.

Calculates atmospheric pressure as a function of elevation with reference to the standard atmosphere. The elevation-dependence of atmospheric pressure is computed by assuming a linear decrease in temperature with elevation and a mean adiabatic lapse rate (Eqn 3, Berberan-Santos et al., 2009):

\[p(z) = p_0 ( 1 - L z / K_0) ^{ G M / (R L) },\]
Parameters:
  • elv – Elevation above sea-level (\(z\), metres above sea level.)

  • core_const – Instance of CoreConst.

Returns:

A numeric value for \(p\) in Pascals.

Examples

>>> # Standard atmospheric pressure, in Pa, corrected for 1000 m.a.s.l.
>>> round(calculate_patm(1000), 2)
90241.54

The solar submodule

The solar submodule provides functions to calculate photosynthetic photon flux density (PPFD), daily solar radiation fluxes and other radiative values.

Several of the calculations share a set of intermediate values (\(r_u, r_v, r_w\)), which are calculated using calculate_ru_rv_intermediates() and calculate_rw_intermediate(). The functions for these calculations have a public function that uses standard input variables, and a second private version of the function that accepts precalculated values of \(r_u, r_v, r_w\). These private functions can be used to avoid repeated calculation of intermediate values where multiple functions need to be run.

Classes:

SolarPositions(latitude, ...], ...)

Solar values for observation locations and times.

Functions:

calculate_daily_solar_radiation(...[, ...])

Calculate daily extraterrestrial solar radiation.

calculate_day_angle(ordinal_date)

Calculate the solar day angle.

calculate_daytime_net_radiation(...[, ...])

Calculate daily net radiation.

calculate_distance_factor(nu[, ...])

Calculates distance factor.

calculate_equation_of_time(day_angle[, coef])

Calculate the equation of time.

calculate_heliocentric_longitudes(...[, ...])

Calculate heliocentric longitude and anomaly.

calculate_local_hour_angle(current_time, ...)

Calculate the local hour angle.

calculate_net_longwave_radiation(...[, coef])

Calculate net longwave radiation.

calculate_net_radiation_crossover_hour_angle(...)

Calculates the net radiation crossover hour angle.

calculate_nighttime_net_radiation(...[, ...])

Calculates nighttime net radiation.

calculate_ppfd(sunshine_fraction, elevation, ...)

Calculates photosynthetic photon flux density (PPFD).

calculate_ppfd_from_tau_rd(transmissivity, ...)

Calculate photosynthetic photon flux density from intermediates.

calculate_ru_rv_intermediates(declination, ...)

Calculates intermediate values for use in solar radiation calcs.

calculate_rw_intermediate(transmissivity, ...)

Calculate the rw intermediate variable.

calculate_solar_declination(ordinal_date)

Calculate solar declination angle.

calculate_solar_declination_angle(lambda_[, ...])

Calculate the solar declination angle.

calculate_solar_elevation(latitude, ...)

Calculate the solar elevation angle.

calculate_solar_noon(longitude, equation_of_time)

Calculate the solar noon for a given location.

calculate_sunset_hour_angle(declination, ...)

Calculate the sunset hour angle.

calculate_transmissivity(sunshine_fraction, ...)

Calculate atmospheric transmissivity.

class SolarPositions(latitude: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.floating]], longitude: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.floating]], datetime: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.datetime64]], core_const: ~pyrealm.constants.core_const.CoreConst = <factory>)

Solar values for observation locations and times.

This class encapsulates the calculation of key solar position data for observations given arrays of the latitude and longitudes of locations and the datetimes of observation.

Example

>>> import numpy as np
>>> sp = SolarPositions(
...     latitude=-35.058333,
...     longitude=147.34167,
...     datetime=np.array([np.datetime64("2024-08-12T10:30:32")]),
... )
>>> sp.decimal_time.round(5)
array([10.50889])
>>> sp.solar_elevation.round(5)
array([0.60252])

Methods:

get_local_standard_meridian()

Calculate local meridian from longitude.

Attributes:

core_const

A core constants instance.

datetime

An array of np.datetime64 values corresponding to observations at the location (local time).

day_angle

The solar day angle of the observations (\(\Gamma\), radians).

decimal_time

An array of decimal hour values calculated from the datetimes.

declination

The declination for the observations ( (\(\delta\), radians).

equation_of_time

The equation of time value for the observations (\(E_t\), minutes).

hour_angle

The local hour angle for the observations (\(h\), radians).

latitude

The latitude of the location in degrees.

latitude_rad

The latitude of the location in radians, calculated automatically.

local_standard_meridian

An array of local meridians given the longitude.

longitude

The longitude of the location in degrees.

longitude_rad

The longitude of the location in radians, calculated automatically.

ordinal_date

An array of ordinal dates calculated from the datetimes.

solar_elevation

The solar elevation for the observations (\(\alpha\), radians) .

solar_noon

The solar noon for the observations (\(t_0\), decimal hour).

get_local_standard_meridian() ndarray[tuple[Any, ...], dtype[floating]]

Calculate local meridian from longitude.

This calculates the approximate local standard meridian given the longitude. The longitudes are simply mapped to the central longitude of the 15° band in which they fall. This is a _very_ approximate estimate of the local time zone.

core_const: CoreConst

A core constants instance.

datetime: ndarray[tuple[Any, ...], dtype[datetime64]]

An array of np.datetime64 values corresponding to observations at the location (local time).

day_angle: ndarray[tuple[Any, ...], dtype[floating]]

The solar day angle of the observations (\(\Gamma\), radians).

decimal_time: ndarray[tuple[Any, ...], dtype[floating]]

An array of decimal hour values calculated from the datetimes.

declination: ndarray[tuple[Any, ...], dtype[floating]]

The declination for the observations ( (\(\delta\), radians).

equation_of_time: ndarray[tuple[Any, ...], dtype[floating]]

The equation of time value for the observations (\(E_t\), minutes).

hour_angle: ndarray[tuple[Any, ...], dtype[floating]]

The local hour angle for the observations (\(h\), radians).

latitude: ndarray[tuple[Any, ...], dtype[floating]]

The latitude of the location in degrees.

latitude_rad: ndarray[tuple[Any, ...], dtype[floating]]

The latitude of the location in radians, calculated automatically.

local_standard_meridian: ndarray[tuple[Any, ...], dtype[floating]]

An array of local meridians given the longitude.

longitude: ndarray[tuple[Any, ...], dtype[floating]]

The longitude of the location in degrees.

longitude_rad: ndarray[tuple[Any, ...], dtype[floating]]

The longitude of the location in radians, calculated automatically.

ordinal_date: ndarray[tuple[Any, ...], dtype[int64]]

An array of ordinal dates calculated from the datetimes.

solar_elevation: ndarray[tuple[Any, ...], dtype[floating]]

The solar elevation for the observations (\(\alpha\), radians) .

solar_noon: ndarray[tuple[Any, ...], dtype[floating]]

The solar noon for the observations (\(t_0\), decimal hour).

calculate_daily_solar_radiation(distance_factor: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, sunset_hour_angle: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, declination: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, latitude: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, day_seconds: float = CoreConst().day_seconds, solar_constant: float = CoreConst().solar_constant) ndarray[tuple[Any, ...], dtype[floating]]

Calculate daily extraterrestrial solar radiation.

This function calculates the daily extraterrestrial solar radition (\(R_d\), J m-2) using Eq. 1.10.3 of Duffie and Beckman (2013):

\[R_d = \left( \frac{n_s}{\pi} \right) G_{sc} \, d_r \left(r_u \, h_s + r_v \sin(h_s) \right),\]

where \(r_u\) and \(r_v\) are the outputs of calculate_ru_rv_intermediates().

Parameters:
  • distance_factor – The distance factor (\(d_r\), unitless)

  • sunset_hour_angle – Local hour angle, (\(h_s\), degrees)

  • declination – solar declination (\(\delta\), degrees)

  • latitude – Site latitude (\(\phi\), degrees)

  • day_seconds – Number of seconds in one solar day (\(n_s\), seconds), defaulting to CoreConst.day_seconds.

  • solar_constant – The solar constant (\(G_{sc}\), W m-2), defaulting to CoreConst.solar_constant.

Returns:

An array of daily solar radiation, J/m^2

calculate_day_angle(ordinal_date: ndarray[tuple[Any, ...], dtype[int64]] | DataArray) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the solar day angle.

Calculates the solar day angle (\(\Gamma\), radians) for ordinal dates (‘Julian dates’) using Eqn A18 of De Pury and Farquhar (1997).

\[\Gamma = \frac{2\pi (N - 1)}{365}\]
Parameters:

ordinal_date – The ordinal date for which to calculate the day angle.

Returns:

An array of solar day angles in radians.

calculate_daytime_net_radiation(net_longwave_radiation: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, crossover_hour_angle: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, declination: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, latitude: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, transmissivity: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, distance_factor: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, shortwave_albedo: float = CoreConst().shortwave_albedo, solar_constant: float = CoreConst().solar_constant, day_seconds: float = CoreConst().day_seconds) ndarray[tuple[Any, ...], dtype[floating]]

Calculate daily net radiation.

Calculates the daily net radiation (\(R_{nd}\), J m-2) as:

\[R_{nd} = \left( \frac{n_s}{\pi} \right) \left( h_n (r_w \, r_u - R_{nl}) + r_w \, r_v \, sin(h_n)) \right)\]

where \(r_u\), \(r_v\) and \(r_w\) are the outputs of calculate_ru_rv_intermediates() and calculate_rw_intermediate().

Parameters:
  • net_longwave_radiation – Net longwave radiation, (\(R_{nl}\), W m-2)

  • crossover_hour_angle – Crossover hour angle, (\(h_n\), degrees)

  • declination – solar declination (\(\delta\), degrees)

  • latitude – Site latitude (\(\phi\), degrees)

  • transmissivity – Bulk transmissivity (\(\tau\), unitless)

  • distance_factor – The distance factor (\(d_r\), unitless)

  • shortwave_albedo – The shortwave albedo (\(A_{sw}\), unitless), defaulting to CoreConst.shortwave_albedo.

  • solar_constant – The solar constant, (\(G_{sc}\), W/m^2), defaulting to CoreConst.solar_constant.

  • day_seconds – Number of seconds in one solar day (\(n_s\), seconds), defaulting to CoreConst.day_seconds.

Result:

An array of daily net radiation.

calculate_distance_factor(nu: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, solar_eccentricity: float = CoreConst().solar_eccentricity) ndarray[tuple[Any, ...], dtype[floating]]

Calculates distance factor.

This function calculates the distance factor \(d_r\) using the method of Berger et al. (1993).

\[d_r = \left( 1.0 \mathbin{/} \left(\frac{1.0 - e^2} {1.0 + e \cos\left(\nu) \right)} \right) \right)^2\]
Parameters:
  • nu – Heliocentric true anomaly (\(\nu\), degrees)

  • solar_eccentricity – Solar eccentricity (\(e\)), defaulting to CoreConst.solar_eccentricity.

Returns:

An array of distance factors

calculate_equation_of_time(day_angle: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, coef: tuple[float, ...] = CoreConst.equation_of_time_coef) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the equation of time.

Calculates values of the equation of time (\(E_t\), minutes) from the day angle as:

\[E_t = f \left( a + b \cos(\Gamma) + c \sin(\Gamma) + d \cos(2\Gamma) + e \sin(2\Gamma) \right),\]

where \(a,b,c,d,e,f\) are the coefficients of the equation provided in the coef argument. The implementation is taken from Eqn 1.4.1 of Iqbal (1983). Note that Eqn A17 of De Pury and Farquhar (1997) provides an implementation that contains errors.

Parameters:
  • day_angle – The day angle in radians (\(\Gamma\), radians)

  • coef – A tuple of coefficients for the equation of time, defaulting to the value of CoreConst.equation_of_time_coef.

Returns:

An array of equation of time values (minutes)

calculate_heliocentric_longitudes(ordinal_date: ndarray[tuple[Any, ...], dtype[int64]] | DataArray, n_days: ndarray[tuple[Any, ...], dtype[int64]] | DataArray, solar_eccentricity: float = CoreConst().solar_eccentricity, solar_perihelion: float = CoreConst().solar_perihelion) tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]

Calculate heliocentric longitude and anomaly.

This function calculates the heliocentric true anomaly (nu, \(\nu\), degrees) and true longitude (lambda, \(\lambda\), degrees), given the ordinal date in the year and the number of days in the year, following Berger (1978).

Parameters:
  • ordinal_date – The ordinal date

  • n_days – The number of days in the year

  • solar_eccentricity – The solar eccentricity (\(e\)), defaulting to CoreConst.solar_eccentricity.

  • solar_perihelion – The solar perihelion (\(\omega\)), defaulting to CoreConst.solar_perihelion.

Returns:

A tuple of NDArrays containing nu and lambda_.

calculate_local_hour_angle(current_time: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, solar_noon: ndarray[tuple[Any, ...], dtype[floating]] | DataArray) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the local hour angle.

The local hour angle (\(h\), radians) is a measure of time, expressed as an angle indicating the position of the sun relative to solar noon. This function calculates the local hour angle following equation A15 of De Pury and Farquhar (1997).

\[h = \pi \frac{t - t_{0}}{12}\]
Parameters:
  • current_time – Current time values in decimal hours (\(t\)).

  • solar_noon – Solar noon time values in decimal hours (\(t_0\)).

Returns:

The local hour angle in radians

calculate_net_longwave_radiation(sunshine_fraction: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, temperature: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, coef: tuple[float, float] = CoreConst().net_longwave_radiation_coef) ndarray[tuple[Any, ...], dtype[floating]]

Calculate net longwave radiation.

This function calculates the net longwave radiation (\(R_{nl}\), W m-2) following Eqn. 11 and Table 1 of Prentice et al. (1993):

\[R_{nl} = (b + (1.0 - b) n_i) (A - t_c)\]
Parameters:
  • sunshine_fraction – Sunshine fraction of observations (\(n_i\), unitless)

  • temperature – Temperature of observations (\(t_c\), °C).

  • coef – Coefficients \(b\) and \(A\) of the equation,, defaulting to CoreConst.net_longwave_radiation_coef.

Returns:

An array of net longwave radiation.

calculate_net_radiation_crossover_hour_angle(net_longwave_radiation: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, transmissivity: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, distance_factor: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, declination: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, latitude: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, shortwave_albedo: float = CoreConst().shortwave_albedo, solar_constant: float = CoreConst().solar_constant) ndarray[tuple[Any, ...], dtype[floating]]

Calculates the net radiation crossover hour angle.

This function calculates the net radiation crossover hour angle (\(h_n\), degrees) as:

\[h_n = \arccos(\frac{R_{nl} - r_w r_u}{r_w r_v}),\]

where \(r_u\), \(r_v\) and \(r_w\) are the outputs of calculate_ru_rv_intermediates() and calculate_rw_intermediate(). The function is clamped within \([-180°, 180°]\).

Parameters:
  • net_longwave_radiation – Net longwave radiation, (\(R_{nl}\), W m-2)

  • transmissivity – Bulk transmissivity (\(\tau\), unitless)

  • declination – solar declination (\(\delta\), degrees)

  • latitude – Site latitude (\(\phi\), degrees)

  • distance_factor – The distance factor (\(d_r\), unitless)

  • shortwave_albedo – The shortwave albedo (\(A_{sw}, unitless), defaulting to :attr:`CoreConst.shortwave_albedo<pyrealm.constants.CoreConst.shortwave_albedo>\).

  • solar_constant – The solar constant, (\(G_{sc}\), W m-2), defaulting to CoreConst.solar_constant.

Returns:

The net radiation crossover hour angle.

calculate_nighttime_net_radiation(net_longwave_radiation: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, crossover_hour_angle: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, sunset_hour_angle: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, declination: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, latitude: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, transmissivity: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, distance_factor: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, shortwave_albedo: float = CoreConst().shortwave_albedo, solar_constant: float = CoreConst().solar_constant, day_seconds: float = CoreConst().day_seconds) ndarray[tuple[Any, ...], dtype[floating]]

Calculates nighttime net radiation.

This function calculates nighttime net radiation (\(R_{nn}\), J m-2) as:

\[R_{nn} = \left( r_w \, r_v \, (\sin(h_s) - \sin(h_n)) + r_w \, r_u \, (h_s - h_n) - R_{nl} \, (\pi - h_n) \right) \left( \frac{n_s}{\pi} \right),\]

where \(r_u\), \(r_v\) and \(r_w\) are the outputs of calculate_ru_rv_intermediates() and calculate_rw_intermediate().

Parameters:
  • net_longwave_radiation – Net longwave radiation, (\(R_{nl}\), W m-2)

  • crossover_hour_angle – Crossover hour angle, (\(h_n\), degrees)

  • sunset_hour_angle – sunset hour angle, (\(h_s\), degrees)

  • declination – solar declination (\(\delta\), degrees)

  • latitude – Site latitude (\(\phi\), degrees)

  • transmissivity – Bulk transmissivity (\(\tau\), unitless)

  • distance_factor – The distance factor (\(d_r\), unitless)

  • shortwave_albedo – The shortwave albedo (\(A_{sw}\), unitless), defaulting to CoreConst.shortwave_albedo.

  • solar_constant – The solar constant (\(G_{sc}\), W m-2), defaulting to CoreConst.solar_constant.

  • day_seconds – Number of seconds in one solar day (\(n_s\), seconds), defaulting to CoreConst.day_seconds.

Returns:

An array of nighttime net radiation in J m-2

calculate_ppfd(sunshine_fraction: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, elevation: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, latitude: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, ordinal_date: ndarray[tuple[Any, ...], dtype[int64]] | DataArray, n_days: ndarray[tuple[Any, ...], dtype[int64]] | DataArray, const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculates photosynthetic photon flux density (PPFD).

This function calculates photosynthetic photon flux density (PPFD, mol m-2) for a given ordinal date in a year for a location at a given latitude and elevation. The PPFD value is moderated by the sunshine fraction for the day to account for effects of cloud cover.

Parameters:
  • sunshine_fraction – Daily sunshine fraction of observations (\(s_f\), unitless).

  • elevation – Elevation of observations, metres.

  • latitude – The latitude of observations (degrees).

  • ordinal_date – The ordinal_date of the observation.

  • n_days – Days in the calendar year.

  • const – A CoreConst object.

Returns:

An array of photosynthetic photon flux density.

Example

>>> # Define variable values
>>> # Evaluate function
>>> calculate_ppfd(
...    sunshine_fraction=np.array([1.0]),
...    elevation=np.array([142]),
...    latitude=np.array([37.7]),
...    ordinal_date=np.array([172]),
...    n_days=np.array([366])
... )
array([62.04230021])
calculate_ppfd_from_tau_rd(transmissivity: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, daily_solar_radiation: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, visible_light_albedo: float = CoreConst().visible_light_albedo, swdown_to_ppfd_factor: float = CoreConst().swdown_to_ppfd_factor) ndarray[tuple[Any, ...], dtype[floating]]

Calculate photosynthetic photon flux density from intermediates.

This function calculates photosynthetic photon flux density (PPFD, µmol m-2 s-1) from precalculated values for daily solar radiation and transimissivity.

\[\mathrm{PPFD} = (1.0 \times 10^{-6}) f_{PPFD} (1.0 - A_{vis}) \tau R_{d}\]
Parameters:
  • transmissivity – The bulk transmissivity (\(\tau\), unitless)

  • daily_solar_radiation – Daily solar radiation (\(R_d\), J m-2)

  • visible_light_albedo – The visible light albedo (\(A_{vis}\), unitless), defaulting to CoreConst.visible_light_albedo.

  • swdown_to_ppfd_factor – Conversion factor from W m-2 of sunlight to µmol m-2 s-1 of PPFD (\(f_{PPFD}\)), defaulting to CoreConst.swdown_to_ppfd_factor.

Returns:

An array of photosynthetic photon flux density, µmol m-2 s-1.

calculate_ru_rv_intermediates(declination: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, latitude: ndarray[tuple[Any, ...], dtype[floating]] | DataArray) tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]

Calculates intermediate values for use in solar radiation calcs.

This function calculates \(r_u\) and \(r_v\) which are dimensionless intermediate values calculated from the solar declination angle (\(\delta\)) and the observation latitude (\(\phi\)).

\[ \begin{align*} r_u &= \sin(\delta) \sin(\phi) \\ r_v &= \cos(\delta) \cos(\phi) \end{align*} \]
Parameters:
  • declination – Solar declination (\(\delta\), degrees)

  • latitude – Observation latitude (\(\phi\), degrees)

Returns:

A tuple of \(r_u\) and \(r_v\).

calculate_rw_intermediate(transmissivity: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, distance_factor: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, shortwave_albedo: float = CoreConst().shortwave_albedo, solar_constant: float = CoreConst().solar_constant) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the rw intermediate variable.

This function calculates the widely used variable substitute rw (\(r_w\), W m-2) as:

\[r_w = (1.0 - A_{sw}) \, \tau \, G_{sc} \, d_r\]
Parameters:
  • transmissivity – Bulk transmissivity, (\(\tau\), unitless).

  • distance_factor – Distance factor (\(d_r\), unitless)

  • shortwave_albedo – The shortwave albedo (\(A_{sw}\), unitless), defaulting to CoreConst.shortwave_albedo.

  • solar_constant – The solar constant (\(G_{sc}\), W m-2), defaulting to CoreConst.solar_constant.

Returns:

An array of the intermediate variable \(r_w\) in W m-2.

calculate_solar_declination(ordinal_date: ndarray[tuple[Any, ...], dtype[int64]] | DataArray) ndarray[tuple[Any, ...], dtype[floating]]

Calculate solar declination angle.

Calculates the solar declination angle (\(\delta\), radians) from the ordinal date, following Eqn A14 of De Pury and Farquhar (1997).

\[\delta = -23.4 \frac{\pi}{180} \cos \left(\frac{2 \pi (td + 10)}{365}\right)\]
Parameters:

ordinal_date – The ordinal dates of observations (\(t_d\), days).

Returns:

Solar declination angles in radians.

calculate_solar_declination_angle(lambda_: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, solar_obliquity: float = CoreConst().solar_obliquity) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the solar declination angle.

This function calculates the solar declination angle (\(\delta\), degrees) following Woolf (1968).

\[\delta = \arcsin(\sin(\lambda)\sin(\epsilon))\]
Parameters:
  • lambda – The heliocentric longitude (\(\lambda\), degrees)

  • solar_obliquity – Solar obliquity (\(\epsilon\)), defaulting to CoreConst.solar_obliquity.

Returns:

An array of solar declination angle delta

calculate_solar_elevation(latitude: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, declination: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, hour_angle: ndarray[tuple[Any, ...], dtype[floating]] | DataArray) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the solar elevation angle.

The solar elevation angle (\(\alpha\), radians) is the angle between the horizon and the sun, giving the height of the sun in the sky at a given time, following Eqn A13 of De Pury and Farquhar (1997):

\[\alpha = \arcsin(\sin(\phi) \sin(\delta) + \cos(\phi) \cos(\delta) \cos(h))\]
Parameters:
  • latitude – Observation latitudes (\(\phi\), radians).

  • declination – Solar declination angles (\(\delta\), radians).

  • hour_angle – Observation hour angles (\(\h\), radians).

Returns:

Solar elevation angles in radians.

calculate_solar_noon(longitude: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, equation_of_time: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, standard_longitude: ndarray[tuple[Any, ...], dtype[floating]] | DataArray = np.array([0])) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the solar noon for a given location.

The solar noon (\(t_0\), decimal hour) is the time of day when the sun is at its highest point in the sky for a given location. This function calculates the solar noon by adjusting the standard noon time (12:00 PM) given the local longitude and the equation of time for the day of observation, following Eqn. A16 of De Pury and Farquhar (1997).

\[t_{0} = 12 + \frac{4 \cdot -(L_{e} - L_{s}) - E_{t}}{60}\]
Parameters:
  • longitude – The local longitude in degrees (\(L_e\), degrees).

  • equation_of_time – The equation of time given the ordinal date of the observation (\(E_t\), minutes).

  • standard_longitude – The standard meridian for the observation (\(L_s\), degrees), defaulting to the Greenwich meridan.

Returns:

The solar noon time in decimal hours.

calculate_sunset_hour_angle(declination: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, latitude: ndarray[tuple[Any, ...], dtype[floating]] | DataArray) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the sunset hour angle.

This function calculates the sunset hour angle \(h_s\) in degrees following Eq. 3.22 of Stine and Geyer (2001) as:

\[hs = \arccos(-1.0 \frac{r_u}{r_v}),\]

where \(r_u\) and \(r_v\) are the outputs of calculate_ru_rv_intermediates(). The function is clamped within \([-180°, 180°]\).

Parameters:
  • declination – The solar declination angle (\(\delta\), degrees)

  • latitude – The site latitude (\(\phi\), degrees)

Returns:

An array of local hour angle, degrees

calculate_transmissivity(sunshine_fraction: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, elevation: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, coef: tuple[float, ...] = CoreConst.transmissivity_coef) ndarray[tuple[Any, ...], dtype[floating]]

Calculate atmospheric transmissivity.

This function calculates atmospheric transmissivity (\(\tau\), unitless) following Eqn. 11 of Linacre (1968) and Eq 2 of Allen (1996).

\[\tau = (c + d \cdot s_f) (1.0 + f \cdot H)\]
Parameters:
  • sunshine_fraction – Daily sunshine fraction of observations, (\(s_f\), unitless)

  • elevation – Elevation of observations, (\(H\), metres)

  • coef – Coefficients of the equation (\(c, d, f\)), , defaulting to CoreConst.transmissivity_coef.

Returns:

An array of bulk transmissivity

The time_series submodule

This module provides general tools for working with time series data. It currently provides the AnnualValueCalculator, which is used to calculate annual means and totals of time series data, and broadcast_time(), which is used to broadcast arrays over the time axis.

Classes:

AnnualValueCalculator(data_shape, timing[, ...])

Annual means and totals from time series data.

Functions:

broadcast_time(values, shape)

Broadcast an array along the time (zeroth) axis.

class AnnualValueCalculator(data_shape: tuple[int, ...], timing: AcclimationModel | ndarray[tuple[Any, ...], dtype[datetime64]], subset_mask: ndarray[tuple[Any, ...], dtype[bool]] | None = None, endpoint: datetime64 | None = None)

Annual means and totals from time series data.

This class calculates annual means and totals from time series data and is designed to handle time series data sampled at a wide range of intervals. If time series are sampled at intervals that do not map neatly onto years - such as weekly or fortnightly data - or do not have precisely uneven sampling - such as monthly data, then calculating annual means and sums becomes awkward.

This class handles the issue by taking a set of datetimes for observations and mapping them onto years. Observations that span year boundaries are included in both years and the class calculates two forms of weightings:

  • Duration weights provide the actual duration of each observation within a year, and are used to calculate weighted means

  • Proportion weights provide the proportion of an observation in a given year and are used to calculate annual sums.

The class must be created for a specific target dataset with a known array shape (the data_shape argument), with the first axis representing time. The timings argument then sets the timing of observations along that axis. Lastly, the optional subset_mask argument can be used to define a subset of observations to be used when calculating totals and means. This subset mask must be an array that can be broadcast to the data shape. The datetimes do not have to completely sample all years: the year_completeness attribute records what fraction of a year is covered by the datetimes.

As an example, 10 years of monthly data for a 5 x 5 spatial grid of sites could have:

  • a data shape of (120, 5, 5),

  • a timings array with shape (120,), and then

  • a subset mask with shape (120, 1, 1), providing a single subset mask for use at all 25 grid cells. Alternatively, a subset mask with shape (120, 5, 5) would provide site specific subsets.

With uneven sampling - such as monthly - an explicit endpoint must be provided to set the duration of the final observation.

Once an instance is created, then the AnnualValueCalculator.get_annual_totals() and AnnualValueCalculator.get_annual_means() methods can be used to calculate the actual summary statistics for arrays of values that match the data shape.

The timings argument can either directly provide an array of numpy.datetime64 values, or use an existing AcclimationModel object, which contains such an array. The subset_mask argument must be a boolean array.

Example

>>> # Three years of monthly data - 36 observations
>>> datetimes = np.arange(
...     np.datetime64('2000-01'),
...     np.datetime64('2003-01'),
...     np.timedelta64(1, "M")
... )
>>> # Monthly data is uneven - requires an explicit endpoint.
>>> avc = AnnualValueCalculator(
...    data_shape=(36,),
...    timing=datetimes,
...    endpoint=np.datetime64('2003-01')
... )
>>> avc.year_completeness
array([1., 1., 1.])
Parameters:
  • timing – A AcclimationModel instance or a one-dimensional array of numpy.datetime64 values.

  • data_shape – A tuple giving the shape of the data to be used with an instance.

  • subset – An array of numpy.bool_ values to be used in excluding observations from summary stat calculations within a year.

  • endpoint – A single numpy.datetime64 value, required to provide an explicit endpoint for observations with uneven frequency.

Methods:

get_annual_means(values[, within_subset])

Get annual means from an array of values.

get_annual_totals(values[, within_subset])

Get annual totals from an array of values.

Attributes:

datetimes

The start datetime of observations taking from the initial timings

duration_seconds

The duration of each observation in seconds.

duration_weights

A list of arrays giving the number of seconds that each observation within a year contributes to that year.

endpoint

A datetime giving of the end of the last observation.

fractional_weights

A list of arrays giving the fraction of each observation within a year that falls in the year.

indexing

Pairs of integers giving start and end indices to extract consecutive years of data from the time series.

n_obs

The number of observations in the time series.

shape

The shape of data value arrays accepted by an instance.

subset_mask

The initial input array for the subset mask.

subset_mask_by_year

A list of arrays giving the subset mask subarrays for each year.

year_completeness

Provides the fractional coverage of observations for each year.

year_n_days

The total number of days in each year in the time series.

year_n_growing_days

The total number of growing days for each year in the time series.

year_total_seconds

The total number of seconds for each year in the time series.

years

The covered years as np.datetime64 at year precision.

get_annual_means(values: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, within_subset: bool = False) ndarray[tuple[Any, ...], dtype[floating]]

Get annual means from an array of values.

Annual mean values are calculated using a weighted mean of values falling within each year, using the total __duration__ of each observation within the year as a weight. If an observation spans two years, the observation is included in the mean for both years, weighted by the total duration within the year. The method handles missing data (np.nan).

The input values must be an array with a shape that matches the data_shape used to create the AnnualValueCalculator instance. Annual means are calculated along the first axis of the input array and array dimensions are preserved along all but that first axis. For example, monthly data values over 10 years for a 3 x 3 grid of cells might have shape (120, 3, 3) - the resulting array of mean values would have shape (10, 3, 3).

The within_subset argument can be used to calculate mean values only from observations included in the subset_mask defined when creating the AnnualValueCalculator instance.

Note

The method returns values for incomplete years, but obviously these may be be biased as a result.

Example

>>> # Three years of monthly data
>>> datetimes = np.arange(
...     np.datetime64('2000-01'),
...     np.datetime64('2003-01'),
...     np.timedelta64(1, "M")
... )
>>> # Monthly data is uneven - requires an explicit endpoint.
>>> avc = AnnualValueCalculator(
...    data_shape=(36,),
...    timing=datetimes,
...    endpoint=np.datetime64('2003-01')
... )
>>> # Note that the means are weighted by the actual durations of months.
>>> avc.get_annual_means(np.arange(0, 36)).round(4)
array([ 5.5137, 17.526 , 29.526 ])
Parameters:
  • values – The data to summarize by year

  • within_subset – Should the mean only include values within the subset mask.

get_annual_totals(values: ndarray[tuple[Any, ...], dtype[floating]] | DataArray, within_subset: bool = False) ndarray[tuple[Any, ...], dtype[floating]]

Get annual totals from an array of values.

Annual totals are calculated by splitting the values into subarrays of values by year and then calculating the sum along the first axis. When observations span two years, they are included in both sets of annual data, but the sum is weighted by the fraction of the observation in each year to partition the value between years. The method handles missing data (np.nan) but obviously the resulting annual total will be reduced.

The input values must be an array with a shape that matches the data_shape used to create the AnnualValueCalculator instance. Annual totals are calculated along the first axis of the input array and array dimensions are preserved along all but that first axis. For example, monthly data values over 10 years for a 3 x 3 grid of cells might have shape (120, 3, 3) - the resulting array of total values would have shape (10, 3, 3).

The within_subset argument can be used to calculate totals values only from observations included in the subset_mask defined when creating the AnnualValueCalculator instance.

Note

The method returns values for incomplete years, but obviously these will be an underestimate. A simple correction would be to divide the totals by the year completeness to scale back up to full annual values.

Example

>>> # Three years of monthly data with incomplete years at start and end
>>> datetimes = np.arange(
...     np.datetime64('2000-07'),
...     np.datetime64('2003-07'),
...     np.timedelta64(1, "M")
... )
>>> # Monthly data is uneven - requires an explicit endpoint.
>>> avc = AnnualValueCalculator(
...    data_shape=(36,),
...    timing=datetimes,
...    endpoint=np.datetime64('2003-07')
... )
>>> # Note that the means are weighted by the actual durations of months.
>>> avc.get_annual_totals(np.arange(0, 36)).round(4)
array([ 15., 138., 282., 195.])
>>> # Year completeness: 184/366 days in 2000, 181/365 days in 2003.
>>> avc.year_completeness.round(4)
array([0.5027, 1.    , 1.    , 0.4959])
Parameters:
  • values – The data to summarize by year

  • within_subset – Should the mean only include values within the subset mask.

datetimes: ndarray[tuple[Any, ...], dtype[datetime64]]

The start datetime of observations taking from the initial timings

duration_seconds: ndarray[tuple[Any, ...], dtype[int64]]

The duration of each observation in seconds.

duration_weights: list[ndarray[tuple[Any, ...], dtype[int64]]]

A list of arrays giving the number of seconds that each observation within a year contributes to that year.

endpoint: datetime64

A datetime giving of the end of the last observation.

fractional_weights: list[ndarray[tuple[Any, ...], dtype[floating]]]

A list of arrays giving the fraction of each observation within a year that falls in the year.

indexing: list[tuple[int, int]]

Pairs of integers giving start and end indices to extract consecutive years of data from the time series.

n_obs: int

The number of observations in the time series.

shape: tuple[int, ...]

The shape of data value arrays accepted by an instance.

subset_mask: ndarray[tuple[Any, ...], dtype[bool]]

The initial input array for the subset mask.

subset_mask_by_year: list[ndarray[tuple[Any, ...], dtype[bool]]]

A list of arrays giving the subset mask subarrays for each year.

year_completeness: ndarray[tuple[Any, ...], dtype[floating]]

Provides the fractional coverage of observations for each year.

year_n_days: ndarray[tuple[Any, ...], dtype[floating]]

The total number of days in each year in the time series.

year_n_growing_days: ndarray[tuple[Any, ...], dtype[floating]]

The total number of growing days for each year in the time series. If the growing_season input varies within days, these values can contain non-integer values.

year_total_seconds: ndarray[tuple[Any, ...], dtype[int64]]

The total number of seconds for each year in the time series.

years: ndarray[tuple[Any, ...], dtype[datetime64]]

The covered years as np.datetime64 at year precision.

broadcast_time(values: ndarray[tuple[Any, ...], dtype[_ScalarT]], shape: tuple[int, ...]) ndarray[tuple[Any, ...], dtype[_ScalarT]]

Broadcast an array along the time (zeroth) axis.

The values array must be broadcastable to the full shape, however it does not need the full set of dimensions as defined by shape. The returned array will have the full set of dimensions, and be broadcast along just the zeroth axis.

Example

>>> broadcast_time(np.ones((1,3)), (2,3))
array([[1., 1., 1.],
       [1., 1., 1.]])
>>> broadcast_time(np.ones(3), (2,2,3)).shape
(2, 1, 3)
Parameters:
  • values – The array to broadcast.

  • shape – The full n-dimensional shape, where the first value is the length of the time axis to broadcast over.

The utilities submodule

The utilities submodule provides shared utility functions used to:

  • check input arrays to functions have congruent shapes

  • summarize object attributes in a tabular format

  • share efficient implementations of core components, such as Horner-form polynomial calculation.

Functions:

check_input_shapes(*args[, shape])

Check sets of input variables have congruent shapes with equal dimensions.

evaluate_horner_polynomial(x, cf)

Evaluates a polynomial with coefficients cf at x using Horner's method.

exponential_moving_average(values[, ...])

Apply an exponential moving average to a variable.

summarize_attrs(obj, attrs[, dp, repr_head])

Print a summary table of object attributes.

check_input_shapes(*args: float | int | generic | ndarray | None, shape: tuple[int, ...] | None = None) tuple

Check sets of input variables have congruent shapes with equal dimensions.

This helper function validates inputs to check that they are either scalars or arrays and that any arrays have the same number of dimensions and are broadcastable to the same shape. It returns a tuple of the common shape of the arguments, which is (1,) if all the arguments are scalar.

Parameters:
  • *args – A set of numpy arrays or scalar values

  • shape – An optional expected result for the common shape.

Returns:

The common shape of any array inputs or 1 if all inputs are scalar.

Raises:

ValueError – if the inputs contain arrays of differing shapes or dimensions.

Examples

>>> check_input_shapes(np.array([1,2,3]), 5)
(3,)
>>> check_input_shapes(4, 5)
(1,)
>>> check_input_shapes(np.ones((1,2,3)), np.ones((1,2,4)))
Traceback (most recent call last):
...
ValueError: Inputs contain arrays of different shapes.
>>> check_input_shapes(np.ones((1,2,3)), np.ones((2,3)))
Traceback (most recent call last):
...
ValueError: Inputs contain arrays of different dimensions.
evaluate_horner_polynomial(x: ndarray[tuple[Any, ...], dtype[floating]], cf: Sequence | ndarray[tuple[Any, ...], dtype[floating]]) ndarray[tuple[Any, ...], dtype[floating]]

Evaluates a polynomial with coefficients cf at x using Horner’s method.

Horner’s method is a fast way to evaluate polynomials, especially for large degrees, that can be evaluated efficiently using the following rearrangement to avoid taking large powers.

\[ \begin{align*} p(x) &= 5 + 4x + 3x^2 + 2x^3\\ &= 5 + x(4 + x(3 + 2x)) \end{align*} \]
Parameters:
  • x – The values at which to evaluate the polynomial

  • cf – The coefficients of the polynomial, ordered from the lowest (constant) to the highest degree.

exponential_moving_average(values: ndarray[tuple[Any, ...], dtype[floating]], initial_values: ndarray[tuple[Any, ...], dtype[floating]] | None = None, alpha: float = 0.067, allow_holdover: bool = False) ndarray[tuple[Any, ...], dtype[floating]]

Apply an exponential moving average to a variable.

This function implements an exponential moving average for an input array, using the parameter alpha (\(\alpha\)) to control the speed of convergence of the smoothed values (\(S\)) to the input values (\(X\)):

\[S_{t} = S_{t-1}(1 - \alpha) + X_{t} \alpha\]

For \(t_{0}\), the first value in the input values is used (\(S_{0} = X_{0}\)), unless alternative initial values are provided. The values array can have multiple dimensions but the smoothing is only applied along the first dimension, which is assumed to represent time.

By default, the values array must not contain any missing values (numpy.nan). However, when allow_holdover is set then the function attempts to handle missing data by holding over the last non-missing smoothed value. This obviously cannot replace missing data at the start of a sequence, but after this if the current input value is missing, then the previous smoothed value will be recycled until it can next be updated from the input data.

Current input data (\(X_{t}\))

NA

not NA

Previous smoothed data (\(S_{t-1}\))

NA

NA

X_{t}

not NA

\(S_{t-1}\)

\(S_{t-1}(1-a) + X_{t}a\)

Parameters:
  • values – The values to apply the memory effect to.

  • initial_values – Values to use at \(S_{0}\) in place of \(X_{0}\).

  • alpha – The relative weight applied to the most recent observation.

  • allow_holdover – Allow missing values to be filled by holding over earlier values.

Returns:

An array of the same shape as values with the memory effect applied.

summarize_attrs(obj: object, attrs: tuple[tuple[str, str], ...], dp: int = 2, repr_head: bool = True) None

Print a summary table of object attributes.

This helper function prints a simple table of the mean, min, max and nan count for named attributes from an object for use in class summary functions.

Parameters:
  • obj – An object with attributes to summarize

  • attrs – A tuple of 2-tuples giving attribute names and units.

  • dp – The number of decimal places used in rounding summary stats.

  • repr_head – A boolean indicating whether to show the object representation before the summary table.

The water submodule

This module provides implementations of functions for the calculation of water density and viscosity. The methods vary greatly in their complexity also in whether they correct for the effects of atmospheric pressure. The most precise implementations are several orders of magnitude slower than the fastest, and differences in prediction are usually slight.

The various methods are organised into two registry dictionaries that provide a simple lookup method to select the implementation for use through settings to a CoreConst instance. The two decorators register_density_method() and register_viscosity_method() are used to add methods to the respective DENSITY_METHODS and VISCOSITY_METHODS registries. In order to maintain a consistent API for calling methods, all methods within a registry must use the same function signature even if the method does not use all of the parameters: functions that do not correct for atmospheric pressure must still require that argument.

These registries can be extended by users to provide their own calculations of density and viscosity.

Functions:

calculate_density_h2o(tc, patm[, ...])

Calculate water density.

calculate_density_h2o_chen(tc, patm[, ...])

Calculate the density of water using Chen et al 2008.

calculate_density_h2o_fisher(tc, patm[, ...])

Calculate water density.

calculate_density_h2o_jones_harris_eq6(tc, patm)

Calculate density of water following Jones and Harris Eqn 6.

calculate_density_h2o_jones_harris_eq8(tc, patm)

Calculate density of water following Jones and Harris Eqn 8.

calculate_density_h2o_kell(tc, patm[, ...])

Calculate density of water by Kell's method.

calculate_viscosity_h2o(tk, patm[, core_const])

Calculate the viscosity of water.

calculate_viscosity_h2o_daubert_danner(tk, patm)

Calculate viscosity of water using the Daubert Danner form.

calculate_viscosity_h2o_girifalco(tk, patm)

Calculate viscosity of water using the Girifalco form.

calculate_viscosity_h2o_huber(tk, patm[, ...])

Calculate the viscosity of water.

calculate_viscosity_h2o_reid(tk, patm[, ...])

Calculate viscosity of water using the Reid form.

calculate_viscosity_h2o_viswanath_natarajan(tk, ...)

Calculate viscosity of water using the Viswanath and Natarajan form.

calculate_viscosity_h2o_vogel(tk, patm[, ...])

Calculate viscosity of water using the Vogel-Fulcher-Tammann equation.

calculate_water_molar_volume(tc, patm[, ...])

Calculate the volume of a mole of water at a given temperature and pressure.

convert_water_mm_to_moles(water_mm, tc, patm)

Convert water in mm per square meter to moles.

convert_water_moles_to_mm(water_moles, tc, patm)

Convert water in moles to mm per square meter.

register_density_method(method_name)

Registration decorator for water density functions.

register_viscosity_method(method_name)

Registration decorator for water viscosity functions.

Data:

DENSITY_FUNCTION_SIGNATURE

The expected signature for registered water density functions, as tuples of parameter name and type annotation.

DENSITY_METHODS

A registry for functions calculating water density.

VISCOSITY_FUNCTION_SIGNATURE

The expected signature for registered water viscosity functions, as tuples of parameter name and type annotation.

VISCOSITY_METHODS

A registry for functions calculating water viscosity.

calculate_density_h2o(tc: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst(), safe: bool = True) ndarray[tuple[Any, ...], dtype[floating]]

Calculate water density.

Calculates the density of water as a function of temperature and atmospheric pressure (in kg/m^3), using the method specified in CoreConst.water_density_method, defaulting to jones_harris_eq6 (calculate_density_h2o_jones_harris_eq6()).

Parameters:
  • tc – air temperature, °C

  • patm – atmospheric pressure, Pa

  • core_const – Instance of CoreConst

  • safe – Prevents the function from estimating density below -30°C, where the functions are numerically unstable.

Returns:

Water density in kg/m^3.

Raises:

ValueError – if tc contains values below -30°C and safe is True, or if the inputs have incompatible shapes.

Examples

>>> calculate_density_h2o(20, 101325).round(3)
np.float64(998.201)
calculate_density_h2o_chen(tc: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the density of water using Chen et al 2008.

This function calculates the density of water (\(\rho\), kg/m^3) as a function of temperature (\(T\), °C) and atmospheric pressure (\(P\), Pa) following Chen et al. (2008).

Warning

The predictions from this function are numerically unstable around -58°C.

Parameters:
  • tc – Air temperature (°C)

  • patm – Atmospheric pressure (Pa)

  • core_const – Instance of CoreConst, providing the polynomial coefficients for the Chen et al. (2008) equations.

Examples

>>> calculate_density_h2o_chen(
...     np.array([20]), np.array([101325])
... ).round(3)
array([998.25])
calculate_density_h2o_fisher(tc: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate water density.

Calculates the density of water (\(\rho\), kg/m^3) as a function of temperature (\(T\), °C) and atmospheric pressure (\(P\), Pa), using the Tumlirz Equation and coefficients calculated by Fisher and Dial Jr (1975).

Warning

The predictions from this function are unstable around -45°C.

Parameters:
  • tc – air temperature, °C

  • patm – atmospheric pressure, Pa

  • core_const – Instance of CoreConst, providing the polynomial coefficients for the Fisher and Dial Jr (1975) equations.

Examples

>>> calculate_density_h2o_fisher(20, 101325).round(3)
np.float64(998.206)
calculate_density_h2o_jones_harris_eq6(tc: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate density of water following Jones and Harris Eqn 6.

This function calculates water density as a function of temperature (\(T\), °C) following Eqn 6. in (Jones and Harris, 1992).

\[\rho_{as} = a + bT + cT^2 + dT^3 + eT^4\]

with coefficients \(a,b,c,d,e\) defined in CoreConst.density_jones_harris_rho This is the density of air-saturated water using the ITS-90 definition of temperature and is intended for use in the range 5-40°C. This method does not correct for atmospheric pressure.

Parameters:
  • tc – Temperature in °C

  • patm – Atmospheric pressure in Pa

  • core_const – An instance of CoreConst providing coefficients

Examples

>>> calculate_density_h2o_jones_harris_eq6(
...     np.array([20]), np.array([101325])
... ).round(3)
array([998.201])
calculate_density_h2o_jones_harris_eq8(tc: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate density of water following Jones and Harris Eqn 8.

This function calculates water density (\(\rho\), kg/m3) as a function of temperature (\(T\), °C) and pressure (\(P\), Pa) following Eqn 8. in (Jones and Harris, 1992).

\[\rho_{asc} = \rho_{as}( 1+ \kappa_T(P / 1000 - 101.325))\]

where \(\rho_{as}\) is calculated as in calculate_density_h2o_jones_harris_eq6() and:

\[\kappa_T = a + bT + cT^2 + dT^3 + eT^4\]

with coefficients \(a,b,c,d,e\) defined in CoreConst.density_jones_harris_kappa. This is the density of air-saturated water using the ITS-90 definition of temperature and is intended for use in the range 5-40°C.

Parameters:
  • tc – Temperature in °C

  • patm – Atmospheric pressure in Pa

  • core_const – An instance of CoreConst providing coefficients

Examples

>>> calculate_density_h2o_jones_harris_eq8(
...     np.array([20]), np.array([101325])
... ).round(3)
array([998.201])
calculate_density_h2o_kell(tc: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate density of water by Kell’s method.

This function calculates water density as a function of temperature (\(T\), °C) following Eqn 16. in (Kell, 1975). This method does not correct for atmospheric pressure.

\[\rho = \frac{a + bT + cT^2 + dT^3 + eT^4 + fT^5}{1 + gT}\]

with coefficients \(a,b,c,d,e,f,g\) defined in CoreConst.density_kell

Parameters:
  • tc – Temperature in °C

  • patm – Atmospheric pressure in Pa

  • core_const – An instance of CoreConst providing coefficients

Examples

>>> calculate_density_h2o_kell(np.array([20]), np.array([101325])).round(3)
array([997.936])
calculate_viscosity_h2o(tk: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the viscosity of water.

Calculates the viscosity of water (\(\mu\), Pa s) as a function of temperature (\(T\), K) and atmospheric pressure (\(P\), Pa), using the method specified in CoreConst.water_viscosity_method, defaulting to vogel (calculate_viscosity_h2o_vogel()).

Parameters:
  • tk – air temperature (K)

  • patm – atmospheric pressure (Pa)

  • core_const – Instance of CoreConst

Examples

>>> # Density of water at 20 °C and standard atmospheric pressure:
>>> calculate_viscosity_h2o(293.15, 101325).round(8)
np.float64(0.00100353)
calculate_viscosity_h2o_daubert_danner(tk: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate viscosity of water using the Daubert Danner form.

This function calculates water viscosity (\(\mu\), Pa s) as a function of temperature (\(T\), K) using the Daubert and Danner form from Equation 4.41 of (Viswanath, 2007).

\[\mu = \exp \left( a + \frac{b}{T} + c \ln T + dT^10 \right)\]

with coefficients \(a,b,c,d\) defined in CoreConst.viscosity_daubert_danner

This method does not correct for atmospheric pressure.

Parameters:
  • tk – air temperature (K)

  • patm – Atmospheric pressure in Pa

  • core_const – An instance of CoreConst providing coefficients

Examples

>>> calculate_viscosity_h2o_daubert_danner(
...    tk=np.array([293.15]), patm=np.array([101325])
... ).round(8)
array([0.00103321])
calculate_viscosity_h2o_girifalco(tk: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate viscosity of water using the Girifalco form.

This function calculates water viscosity (\(\mu\), Pa s) as a function of temperature (\(T\), K) using the Girifalco form and coefficients from Equation 4.27 of (Viswanath, 2007).

\[\log \mu = a + \frac{b}{T} + \frac{c}{T^2}\]

with coefficients \(a,b,c\) defined in CoreConst.viscosity_girifalco

This method does not correct for atmospheric pressure.

Parameters:
  • tk – air temperature (K)

  • patm – Atmospheric pressure in Pa

  • core_const – An instance of CoreConst providing coefficients

Examples

>>> calculate_viscosity_h2o_girifalco(
...    tk=np.array([293.15]), patm=np.array([101325])
... ).round(8)
array([0.00100742])
calculate_viscosity_h2o_huber(tk: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the viscosity of water.

Calculates the viscosity of water (\(\mu\), Pa s) as a function of temperature and atmospheric pressure (Huber et al., 2009).

Parameters:
  • tk – air temperature (K)

  • patm – atmospheric pressure (Pa)

  • core_const – Instance of CoreConst

  • simple – Use the simple formulation.

Examples

>>> # Density of water at 20 degrees C and standard atmospheric pressure:
>>> calculate_viscosity_h2o_huber(
...    tk=np.array([293.15]), patm=np.array([101325])
... ).round(8)
array([0.0010016])
calculate_viscosity_h2o_reid(tk: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate viscosity of water using the Reid form.

This function calculates water viscosity (\(\mu\), Pa s) as a function of temperature (\(T\), K) using a polynomial form from Equation 4.39 of (Viswanath, 2007) and coefficients taken from Reid et al.

\[\log \mu = a + \frac{b}{T} + cT + dT^2\]

with coefficients \(a,b,c,d\) defined in CoreConst.viscosity_reid

This method does not correct for atmospheric pressure.

Parameters:
  • tk – air temperature (K)

  • patm – Atmospheric pressure in Pa

  • core_const – An instance of CoreConst providing coefficients

Examples

>>> calculate_viscosity_h2o_reid(
...    tk=np.array([293.15]), patm=np.array([101325])
... ).round(8)
array([0.00101766])
calculate_viscosity_h2o_viswanath_natarajan(tk: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate viscosity of water using the Viswanath and Natarajan form.

This function calculates water viscosity (\(\mu\), Pa s) as a function of temperature (\(T\), K) using the Viswanath and Natarajan form and coefficients from Equation 4.18 of (Viswanath, 2007).

\[\log \mu = a + \frac{b}{c - T}\]

with coefficients \(a,b,c\) defined in CoreConst.viscosity_viswanath_natarajan

This method does not correct for atmospheric pressure.

Parameters:
  • tk – air temperature (K)

  • patm – Atmospheric pressure in Pa

  • core_const – An instance of CoreConst providing coefficients

Examples

>>> calculate_viscosity_h2o_viswanath_natarajan(
...     tk=np.array([293.15]), patm=np.array([101325])
... ).round(8)
array([0.00100576])
calculate_viscosity_h2o_vogel(tk: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate viscosity of water using the Vogel-Fulcher-Tammann equation.

This function calculates water viscosity (\(\mu\), Pa s) as a function of temperature (\(T\), K) using the Vogel-Fulcher-Tammann equation with coefficients from (Viswanath and Natavajan, 1988).

\[\mu = a e^{\frac{b}{T - c}\]

with coefficients \(a,b,c\) defined in CoreConst.viscosity_vogel

This method does not correct for atmospheric pressure.

Parameters:
  • tk – air temperature (K)

  • patm – Atmospheric pressure in Pa

  • core_const – An instance of CoreConst providing coefficients

Examples

>>> calculate_viscosity_h2o_vogel(
...     tk=np.array([293.15]), patm=np.array([101325])
... ).round(8)
array([0.00100353])
calculate_water_molar_volume(tc: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Calculate the volume of a mole of water at a given temperature and pressure.

Parameters:
  • tc – air temperature (°C)

  • patm – atmospheric pressure (Pa)

  • core_const – Instance of CoreConst

Returns:

Water molar volume in mol cm-3, or equivalently mol/mL

Examples

>>> # A mole of water at standard temperature and pressure occupies ~18 cm3.
>>> calculate_water_molar_volume(0, 101235).round(3)
np.float64(18.015)
convert_water_mm_to_moles(water_mm: ndarray[tuple[Any, ...], dtype[floating]], tc: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Convert water in mm per square meter to moles.

This function converts water volumes expressed as mm per m2 into a number of moles of water. It accounts for the changing density of water with temperature and pressure.

Parameters:
  • water_mm – Water volume in mm per square meter

  • tc – air temperature (°C)

  • patm – atmospheric pressure (Pa)

  • core_const – Instance of CoreConst

Returns:

Moles of water (-)

Examples

>>> # At 0°C and 101325 Pa, one mole of water is ~18 g (18 cm3, 0.018 mm m-2).
>>> # So, 1 mm m2 = 1 / 0.018 = ~55 moles.
>>> convert_water_mm_to_moles(water_mm=1, tc=0, patm=101325).round(3)
np.float64(55.508)
convert_water_moles_to_mm(water_moles: ndarray[tuple[Any, ...], dtype[floating]], tc: ndarray[tuple[Any, ...], dtype[floating]], patm: ndarray[tuple[Any, ...], dtype[floating]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[floating]]

Convert water in moles to mm per square meter.

This function converts water volumes expressed as moles into mm per m2. It accounts for the changing density of water with temperature and pressure.

Parameters:
  • water_moles – Water volume in moles

  • tc – air temperature (°C)

  • patm – atmospheric pressure (Pa)

  • core_const – Instance of CoreConst

Returns:

Water volume in mm per m2

Examples

>>> # At 0°C and 101325 Pa, one mole of water is ~18 g (18 cm3, 0.018 mm m-2).
>>> # So, 1 mol = 0.018 mm
>>> convert_water_moles_to_mm(water_moles=1, tc=0, patm=101325).round(3)
np.float64(0.018)
register_density_method(method_name: str) Callable

Registration decorator for water density functions.

Functions decorated with register_density_method are automatically added to the DENSITY_METHODS registry when imported, using the provided name as a key.

Parameters:

method_name – A short name used as a key for the function in the registry.

register_viscosity_method(method_name: str) Callable

Registration decorator for water viscosity functions.

Functions decorated with register_viscosity_method are automatically added to the VISCOSITY_METHODS registry when imported, using the provided name as a key.

Parameters:

method_name – A short name used as a key for the function in the registry.

DENSITY_FUNCTION_SIGNATURE: tuple[tuple[str, str], ...] = (('tc', 'NDArray[np.floating]'), ('patm', 'NDArray[np.floating]'), ('core_const', 'CoreConst'))

The expected signature for registered water density functions, as tuples of parameter name and type annotation.

DENSITY_METHODS: dict[str, Callable] = {'chen': <function calculate_density_h2o_chen>, 'fisher': <function calculate_density_h2o_fisher>, 'jones_harris_eq6': <function calculate_density_h2o_jones_harris_eq6>, 'jones_harris_eq8': <function calculate_density_h2o_jones_harris_eq8>, 'kell': <function calculate_density_h2o_kell>}

A registry for functions calculating water density. All registered functions are expected to have the same signature.

VISCOSITY_FUNCTION_SIGNATURE: tuple[tuple[str, str], ...] = (('tk', 'NDArray[np.floating]'), ('patm', 'NDArray[np.floating]'), ('core_const', 'CoreConst'))

The expected signature for registered water viscosity functions, as tuples of parameter name and type annotation.

VISCOSITY_METHODS: dict[str, Callable] = {'daubert_danner': <function calculate_viscosity_h2o_daubert_danner>, 'girifalco': <function calculate_viscosity_h2o_girifalco>, 'huber': <function calculate_viscosity_h2o_huber>, 'reid': <function calculate_viscosity_h2o_reid>, 'viswanath_natarajan': <function calculate_viscosity_h2o_viswanath_natarajan>, 'vogel': <function calculate_viscosity_h2o_vogel>}

A registry for functions calculating water viscosity. All registered functions are expected to have the same signature.

The xarray submodule

Utilities for handling xarray inputs to functions that expect arrays.

Functions:

get_common_dims(*arrays[, init_dims])

Get the full list of dimensions across all DataArray arguments.

is_arraytype(var)

Checks if a variable is an instance of ArrayType.

xarray_inputs()

Converts any xarray.DataArray inputs to numpy arrays.

get_common_dims(*arrays: ndarray[tuple[Any, ...], dtype[_ScalarT]] | DataArray, init_dims: list[Hashable] | None = None) list[Hashable]

Get the full list of dimensions across all DataArray arguments.

This needs to be called when there are arrays with multiple dtypes that cannot be combined in a single call to xarray_inputs.

Parameters:
  • *arrays – The variables to convert into numpy arrays.

  • init_dims – An optional list of dims to start with.

Returns:

A list of dimension names.

Examples

>>> input_1 = xr.DataArray([], dims="a")
>>> input_2 = xr.DataArray([[]], dims=["b", "a"])
>>> get_common_dims(input_1, input_2)
['a', 'b']
>>> get_common_dims(input_1, input_2, init_dims=["c"])
['c', 'a', 'b']
is_arraytype(var: Any) TypeGuard[ndarray[tuple[Any, ...], dtype[T]] | DataArray]

Checks if a variable is an instance of ArrayType.

xarray_inputs(array1: ArrayType[T], /, *, kwargs: None = None, dims: list[Hashable] | None = None, dims_full: bool = True) ndarray[tuple[Any, ...], dtype[T]]
xarray_inputs(array1: ArrayType[T], array2: ArrayType[T], /, *other_arrays: ArrayType[T], kwargs: None = None, dims: list[Hashable] | None = None, dims_full: bool = True) tuple[ndarray[tuple[Any, ...], dtype[T]], ...]
xarray_inputs(*arrays: ArrayType[T], kwargs: dict[str, ArrayType[T]], dims: list[Hashable] | None = None, dims_full: bool = True) tuple[tuple[ndarray[tuple[Any, ...], dtype[T]], ...], dict[str, ndarray[tuple[Any, ...], dtype[T]]]]

Converts any xarray.DataArray inputs to numpy arrays.

This allows functions that expect numpy arrays to be used directly with xarray DataArrays, simplifying compatibility between data types.

All DataArray inputs will be expanded to have the same set of dimensions. Where the expanded dimensions will have a length of one. Note that the order of dimensions will depend on the order of inputs - the first input will initially define the order and additional dimensions in later inputs will be appended.

No checking of shape consistency is performed - use check_input_shapes for this.

Parameters:
  • *arrays – The variables to convert into numpy arrays.

  • kwargs – An optional dictionary of variables to convert to numpy arrays.

  • dims – An optional list of dimension names to expand DataArrays to.

  • dims_full – If dims is provided this indicates whether this is the full list of dimensions (True) or can be added to by arrays (default: True).

Returns:

The stripped array(s). As a tuple if more than one is provided. As (tuple, dict) if kwargs is provided, with the dict containing the converted kwargs.

Examples

>>> input = xr.DataArray([1, 2, 3])
>>> array = xarray_inputs(input)
>>> type(array)
<class 'numpy.ndarray'>
>>> array
array([1, 2, 3])
ArrayType: TypeAlias = numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[~T]] | xarray.core.dataarray.DataArray

Type for array inputs. A union of numpy arrays and xarray DataArrays.