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 {class}`~pyrealm.core.bounds.Bounds` dataclass used to define bounds for a particular variable.

  • Defines a {class}`~pyrealm.core.bounds.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 {meth}`~pyrealm.core.bounds.BoundsChecker.update` method to overide defaults or add new variables by providing a new Bounds instance.

The {meth}`~pyrealm.core.bounds.BoundsChecker.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 {class}`~pyrealm.core.bounds.Bounds` instances for core variables, keyed by the {attr}`Bounds.var_name<pyrealm.core.bounds.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[_ScalarT]]) ndarray[tuple[Any, ...], dtype[_ScalarT]]

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 {attr}`Bounds.var_name<pyrealm.core.bounds.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]])

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 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:

calc_enthalpy_vaporisation(tc)

Calculate the enthalpy of vaporization.

calc_psychrometric_constant(tc, p[, core_const])

Calculate the psychrometric constant.

calc_saturation_vapour_pressure_slope(tc)

Calculate the slope of the saturation vapour pressure curve.

calc_specific_heat(tc)

Calculate the specific heat of air.

calc_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.

calc_enthalpy_vaporisation(tc: ndarray[tuple[Any, ...], dtype[float64]]) ndarray[tuple[Any, ...], dtype[float64]]

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).

calc_psychrometric_constant(tc: ndarray[tuple[Any, ...], dtype[float64]], p: ndarray[tuple[Any, ...], dtype[float64]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[float64]]

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

calc_saturation_vapour_pressure_slope(tc: ndarray[tuple[Any, ...], dtype[float64]]) ndarray[tuple[Any, ...], dtype[float64]]

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.

calc_specific_heat(tc: ndarray[tuple[Any, ...], dtype[float64]]) ndarray[tuple[Any, ...], dtype[float64]]

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.

calc_vp_sat(ta: ndarray[tuple[Any, ...], dtype[float64]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[float64]]

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])
>>> calc_vp_sat(temp).round(6)
array([2.480904])
>>> from pyrealm.constants import CoreConst
>>> allen = CoreConst(magnus_option='Allen1998')
>>> calc_vp_sat(temp, core_const=allen).round(6)
array([2.487005])
>>> alduchov = CoreConst(magnus_option='Alduchov1996')
>>> calc_vp_sat(temp, core_const=alduchov).round(6)
array([2.481888])
convert_rh_to_vpd(rh: ndarray[tuple[Any, ...], dtype[float64]], ta: ndarray[tuple[Any, ...], dtype[float64]], core_const: CoreConst = CoreConst(), bounds_checker: BoundsChecker = BoundsChecker()) ndarray[tuple[Any, ...], dtype[float64]]

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])
>>> rh_percent = np.array([70])
>>> convert_rh_to_vpd(rh_percent, temp).round(7) 
pyrealm... contains values outside the expected range (0,1). Check units?
array([-171.1823864])
convert_sh_to_vp(sh: ndarray[tuple[Any, ...], dtype[float64]], patm: ndarray[tuple[Any, ...], dtype[float64]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], ta: ndarray[tuple[Any, ...], dtype[float64]], patm: ndarray[tuple[Any, ...], dtype[float64]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], ta: ndarray[tuple[Any, ...], dtype[float64]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[float64]]

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:

calc_patm(elv[, core_const])

Calculate atmospheric pressure from elevation.

calc_patm(elv: ndarray[tuple[Any, ...], dtype[float64]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[float64]]

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(calc_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 nightime 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.float64]], longitude: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]], 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[float64]]

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[float64]]

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

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

An array of decimal hour values calculated from the datetimes.

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

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

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

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

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

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

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

The latitude of the location in degrees.

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

The latitude of the location in radians, calculated automatically.

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

An array of local meridians given the longitude.

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

The longitude of the location in degrees.

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

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[float64]]

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

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

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

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

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]]) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], crossover_hour_angle: ndarray[tuple[Any, ...], dtype[float64]], declination: ndarray[tuple[Any, ...], dtype[float64]], latitude: ndarray[tuple[Any, ...], dtype[float64]], transmissivity: ndarray[tuple[Any, ...], dtype[float64]], distance_factor: ndarray[tuple[Any, ...], dtype[float64]], shortwave_albedo: float = CoreConst().shortwave_albedo, solar_constant: float = CoreConst().solar_constant, day_seconds: float = CoreConst().day_seconds) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], solar_eccentricity: float = CoreConst().solar_eccentricity) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], coef: tuple[float, ...] = CoreConst.equation_of_time_coef) ndarray[tuple[Any, ...], dtype[float64]]

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]], n_days: ndarray[tuple[Any, ...], dtype[int64]], solar_eccentricity: float = CoreConst().solar_eccentricity, solar_perihelion: float = CoreConst().solar_perihelion) tuple[ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]]]

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[float64]], solar_noon: ndarray[tuple[Any, ...], dtype[float64]]) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], temperature: ndarray[tuple[Any, ...], dtype[float64]], coef: tuple[float, float] = CoreConst().net_longwave_radiation_coef) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], transmissivity: ndarray[tuple[Any, ...], dtype[float64]], distance_factor: ndarray[tuple[Any, ...], dtype[float64]], declination: ndarray[tuple[Any, ...], dtype[float64]], latitude: ndarray[tuple[Any, ...], dtype[float64]], shortwave_albedo: float = CoreConst().shortwave_albedo, solar_constant: float = CoreConst().solar_constant) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], crossover_hour_angle: ndarray[tuple[Any, ...], dtype[float64]], sunset_hour_angle: ndarray[tuple[Any, ...], dtype[float64]], declination: ndarray[tuple[Any, ...], dtype[float64]], latitude: ndarray[tuple[Any, ...], dtype[float64]], transmissivity: ndarray[tuple[Any, ...], dtype[float64]], distance_factor: ndarray[tuple[Any, ...], dtype[float64]], shortwave_albedo: float = CoreConst().shortwave_albedo, solar_constant: float = CoreConst().solar_constant, day_seconds: float = CoreConst().day_seconds) ndarray[tuple[Any, ...], dtype[float64]]

Calculates nightime 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[float64]], elevation: ndarray[tuple[Any, ...], dtype[float64]], latitude: ndarray[tuple[Any, ...], dtype[float64]], ordinal_date: ndarray[tuple[Any, ...], dtype[int64]], n_days: ndarray[tuple[Any, ...], dtype[int64]], const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], daily_solar_radiation: ndarray[tuple[Any, ...], dtype[float64]], visible_light_albedo: float = CoreConst().visible_light_albedo, swdown_to_ppfd_factor: float = CoreConst().swdown_to_ppfd_factor) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], latitude: ndarray[tuple[Any, ...], dtype[float64]]) tuple[ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]]]

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[float64]], distance_factor: ndarray[tuple[Any, ...], dtype[float64]], shortwave_albedo: float = CoreConst().shortwave_albedo, solar_constant: float = CoreConst().solar_constant) ndarray[tuple[Any, ...], dtype[float64]]

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]]) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], solar_obliquity: float = CoreConst().solar_obliquity) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], declination: ndarray[tuple[Any, ...], dtype[float64]], hour_angle: ndarray[tuple[Any, ...], dtype[float64]]) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], equation_of_time: ndarray[tuple[Any, ...], dtype[float64]], standard_longitude: ndarray[tuple[Any, ...], dtype[float64]] = np.array([0])) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], latitude: ndarray[tuple[Any, ...], dtype[float64]]) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], elevation: ndarray[tuple[Any, ...], dtype[float64]], coef: tuple[float, ...] = CoreConst.transmissivity_coef) ndarray[tuple[Any, ...], dtype[float64]]

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(timing[, ...])

A calculator class for annual means and totals from time series data.

Functions:

broadcast_time(values, shape)

Broadcast an array along the time (zeroth) axis.

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

A calculator class for annual means and totals from time series data.

This class is used to calculate annual means and totals from time series data. An instance is created by providing a set of timings for the times series data, either as a one-dimensional array of datetimes or as an AcclimationModel instance from a SubdailyPModel, which provides validated datetimes at subdaily temporal resolutions.

The calculation process accounts for observations that span year boundaries, such as fortnightly data, by calculating the duration of each observation within each year. The process also handles unequal sampling intervals - such as monthly data - by calculating the actual duration of observations. However, with uneven sampling, the duration of the last interval is unknown and so an explicit endpoint must be provided.

The indexing of annual subsets of observations, along with the appropriate weightings for observations values in calculating annual values, are calculated when the class is created and then used by the get_annual_means and get_annual_totals methods. Both methods return values for all years sampled: the year_completeness attribute records what fraction of a year has been sampled to give a particular value.

Note

The class handles a wide range of different possible sampling frequencies and calculates weights for observations using the duration of observations with second precision. With uneven durations - such as monthly data - this will give slightly different values to simple means assuming equal duration.

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(datetimes, endpoint=np.datetime64('2003-01'))
>>> avc.year_completeness
array([1., 1., 1.])

Methods:

get_annual_means(values[, within_growing_season])

Get annual means from an array of values.

get_annual_totals(values[, ...])

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.

growing_season

The initial input array of growing season data.

growing_season_by_year

A list of arrays giving the growing season subarrays for each 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.

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[float64]], within_growing_season: bool = False) ndarray[tuple[Any, ...], dtype[floating]]

Get annual means from an array of values.

Average values are calculated weighted by the __duration__ of each observation, including weighting partial observations than span year boundaries. If within_growing_season is True, the weights for observations outside of the observations marked as the growing season are set to zero. The calculation handles missing values.

The annual means can be computed for a range of different sites by using n-D array for values. In this case time is assumed along axis 0 and so an n-D array will be returned with annual values along axis 0.

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(
...     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_growing_season – Should the mean only include values within the growing season.

get_annual_totals(values: ndarray[tuple[Any, ...], dtype[float64]], within_growing_season: bool = False) ndarray[tuple[Any, ...], dtype[floating]]

Get annual totals from an array of values.

The contribution of each observation to the total is weighted by the __fractional__ duration of each observation within each year in order to partition the sum across years correctly. If within_growing_season is True, the weights for observations not identified as growing season values are set to zero. The method handles missing data (np.nan) but obviously the resulting annual total will be reduced.

The annual totals can be computed for a range of different sites by using an n-D array for values. In this case time is assumed along axis 0 and so an n-D array will be returned with annual values along axis 0.

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(
...     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_growing_season – Should the mean only include values within the growing season.

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[float64]]]

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

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

The initial input array of growing season data.

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

A list of arrays giving the growing season subarrays for each 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.

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

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[float64]], cf: Sequence | ndarray[tuple[Any, ...], dtype[_ScalarT]]) ndarray[tuple[Any, ...], dtype[float64]]

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[float64]], initial_values: ndarray[tuple[Any, ...], dtype[float64]] | None = None, alpha: float = 0.067, allow_holdover: bool = False) ndarray[tuple[Any, ...], dtype[float64]]

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

The water submodule contains core functions for calculating the density and viscosity of water given the air temperature and atmospheric pressure.

Functions:

calc_density_h2o(tc, patm[, core_const, safe])

Calculate water density.

calc_density_h2o_chen(tc, p[, core_const])

Calculate the density of water using Chen et al 2008.

calc_density_h2o_fisher(tc, patm[, core_const])

Calculate water density.

calc_viscosity_h2o(tc, patm[, core_const, ...])

Calculate the viscosity of water.

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

Calculate the viscosity of water.

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.

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

Calculate water density.

Calculates the density of water as a function of temperature and atmospheric pressure (in kg/m^3). This function uses either the method provided by Fisher and Dial Jr (1975) (calc_density_h2o_fisher()) or Chen et al. (2008) (calc_density_h2o_chen()).

The constants attribute water_density_method can be used to set which of the fisher or chen methods is used.

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

>>> round(calc_density_h2o(20, 101325), 3)
np.float64(998.206)
calc_density_h2o_chen(tc: ndarray[tuple[Any, ...], dtype[float64]], p: ndarray[tuple[Any, ...], dtype[float64]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[float64]]

Calculate the density of water using Chen et al 2008.

This function calculates the density of water at a given temperature and pressure (kg/m^3) following Chen et al. (2008).

Warning

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

Parameters:
  • tc – Air temperature (°C)

  • p – Atmospheric pressure (Pa)

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

Returns:

Water density in kg/m^3

Raises:

ValueError – if the inputs have incompatible shapes.

Examples

>>> round(calc_density_h2o_chen(20, 101325), 3)
np.float64(998.25)
calc_density_h2o_fisher(tc: ndarray[tuple[Any, ...], dtype[float64]], patm: ndarray[tuple[Any, ...], dtype[float64]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[float64]]

Calculate water density.

Calculates the density of water as a function of temperature and atmospheric pressure (kg/m^3), 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.

Returns:

Water density in kg/m^3.

Raises:

ValueError – if the inputs have incompatible shapes.

Examples

>>> round(calc_density_h2o_fisher(20, 101325), 3)
np.float64(998.206)
calc_viscosity_h2o(tc: ndarray[tuple[Any, ...], dtype[float64]], patm: ndarray[tuple[Any, ...], dtype[float64]], core_const: CoreConst = CoreConst(), simple: bool = False) ndarray[tuple[Any, ...], dtype[float64]]

Calculate the viscosity of water.

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

Parameters:
  • tc – air temperature (°C)

  • patm – atmospheric pressure (Pa)

  • core_const – Instance of CoreConst

  • simple – Use the simple formulation.

Returns:

A float giving the viscosity of water (mu, Pa s)

Examples

>>> # Density of water at 20 degrees C and standard atmospheric pressure:
>>> round(calc_viscosity_h2o(20, 101325), 7)
np.float64(0.0010016)
calc_viscosity_h2o_matrix(tc: ndarray[tuple[Any, ...], dtype[float64]], patm: ndarray[tuple[Any, ...], dtype[float64]], core_const: CoreConst = CoreConst(), simple: bool = False) ndarray[tuple[Any, ...], dtype[float64]]

Calculate the viscosity of water.

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

Parameters:
  • tc – air temperature (°C)

  • patm – atmospheric pressure (Pa)

  • core_const – Instance of CoreConst

  • simple – Use the simple formulation.

Returns:

A float giving the viscosity of water (mu, Pa s)

Examples

>>> # Viscosity of water at 20 degrees C and standard atmospheric pressure:
>>> round(calc_viscosity_h2o(20, 101325), 7)
np.float64(0.0010016)
calculate_water_molar_volume(tc: ndarray[tuple[Any, ...], dtype[float64]], patm: ndarray[tuple[Any, ...], dtype[float64]], 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.
>>> round(calculate_water_molar_volume(0, 101235), 3)
np.float64(18.015)
convert_water_mm_to_moles(water_mm: ndarray[tuple[Any, ...], dtype[float64]], tc: ndarray[tuple[Any, ...], dtype[float64]], patm: ndarray[tuple[Any, ...], dtype[float64]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[float64]]

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.
>>> round(convert_water_mm_to_moles(water_mm=1, tc=0, patm=101325), 3)
np.float64(55.508)
convert_water_moles_to_mm(water_moles: ndarray[tuple[Any, ...], dtype[float64]], tc: ndarray[tuple[Any, ...], dtype[float64]], patm: ndarray[tuple[Any, ...], dtype[float64]], core_const: CoreConst = CoreConst()) ndarray[tuple[Any, ...], dtype[float64]]

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
>>> round(convert_water_moles_to_mm(water_moles=1, tc=0, patm=101325), 3)
np.float64(0.018)