Module mod17.science

Functions

def climatology365(series: Sequence[+T_co], dates: Sequence[+T_co]) ‑> Sequence[+T_co]

Computes a 365-day climatology for different locations from a time series of length T. Ignores leap days. The climatology could then be indexed using ordinals generated by ordinals365().

Parameters

series : numpy.ndarray

T x … array of data

dates : list or tuple

Sequence of datetime.datetime or datetime.date instances

Returns

numpy.ndarray

 

Expand source code

def climatology365(series: Sequence, dates: Sequence) -> Sequence:
    '''
    Computes a 365-day climatology for different locations from a time series
    of length T. Ignores leap days. The climatology could then be indexed
    using ordinals generated by `ordinals365()`.

    Parameters
    ----------
    series : numpy.ndarray
        T x ... array of data
    dates : list or tuple
        Sequence of datetime.datetime or datetime.date instances

    Returns
    -------
    numpy.ndarray
    '''
    # Get first and last day of the year (DOY)
    ordinal = np.array([
        # Finally, subtract 1 from each day in a leap year after Leap Day
        (doy - 1) if ((dates[i].year % 4 == 0) and doy >= 60) else doy
        for i, doy in enumerate([
            # Next, fill in 0 wherever Leap Day occurs
            0 if (dates[i].year % 4 == 0 and doy == 60) else doy
            for i, doy in enumerate([
                # First, convert datetime.datetime to ordinal day-of-year (DOY)
                int(dt.strftime('%j')) for dt in dates
            ])
        ])
    ])
    with warnings.catch_warnings():
        warnings.simplefilter('ignore')
        return np.array([
            np.nanmean(series[ordinal == day,...], axis = 0)
            for day in range(1, 366)
        ])

def daynight_partition(arr_24hr, updown, reducer='mean', missing=(0, 0))

Partitions a 24-hour time series array into daytime and nighttime values, then calculates the mean in each group. Daytime is defined as when the sun is above the horizon; nighttime is the complement.

NOTE: If the sun never rises, the "daytime" average is computed over 24 hours, anyway. In these cases, "nighttime" average is identical to the "daytime" average (i.e., 24 hours of night). If the sun is always above the horizon (i.e., sun never sets) on a given day, the "daytime" average is computed as expected (over 24 hours) and the nighttime average is missing and is instead filled with the second element of the missing argument's sequence.

Parameters

arr_24hr : numpy.ndarray

A size (24 x …) array; the first axis must have 24 elements corresponding to the measurement in each hour

updown : numpy.ndarray

A size (2 x …) array, compatible with arr_24hr, where the first axis has the hour of sunrise and sunset, in that order, for each element

reducer : str

One of "mean" or "sum" indicating whether an average or a total of the daytime/ nighttime values should be calculated; e.g., for "mean", the hourly values from daytime hours are added up and divided by the length of the day (in hours).

missing : tuple

Values to use when the sun is always below or above the horizon for 24 hours (i.e., never rises or never sets); first value is ignored (may be used to fill missing daylight hours in some future version), and the second value is used to fill in missing nighttime hours (Default: (0, 0)).

Returns

numpy.ndarray

A size (2 x …) array where the first axis enumerates the daytime and nighttime mean values, respectively

Expand source code

def daynight_partition(arr_24hr, updown, reducer = 'mean', missing = (0, 0)):
    '''
    Partitions a 24-hour time series array into daytime and nighttime values,
    then calculates the mean in each group. Daytime is defined as when the sun
    is above the horizon; nighttime is the complement.

    NOTE: If the sun never rises, the "daytime" average is computed over 24
    hours, anyway. In these cases, "nighttime" average is identical to the
    "daytime" average (i.e., 24 hours of night). If the sun is always above
    the horizon (i.e., sun never sets) on a given day, the "daytime" average
    is computed as expected (over 24 hours) and the nighttime average is
    missing and is instead filled with the second element of the `missing`
    argument's sequence.

    Parameters
    ----------
    arr_24hr : numpy.ndarray
        A size (24 x ...) array; the first axis must have 24 elements
        corresponding to the measurement in each hour
    updown: numpy.ndarray
        A size (2 x ...) array, compatible with arr_24hr, where the first axis
        has the hour of sunrise and sunset, in that order, for each element
    reducer : str
        One of "mean" or "sum" indicating whether an average or a total of the
        daytime/ nighttime values should be calculated; e.g., for "mean", the
        hourly values from daytime hours are added up and divided by the
        length of the day (in hours).
    missing : tuple
        Values to use when the sun is always below or above the horizon for 24
        hours (i.e., never rises or never sets); first value is ignored (may
        be used to fill missing daylight hours in some future version), and
        the second value is used to fill in missing nighttime hours
        (Default: `(0, 0)`).

    Returns
    -------
    numpy.ndarray
        A size (2 x ...) array where the first axis enumerates the daytime and
        nighttime mean values, respectively
    '''
    assert reducer in ('mean', 'sum'),\
        'Argument "reducer" must be one of: "mean", "sum"'
    # Prepare single-valued output array
    arr_daytime = np.zeros(arr_24hr.shape[1:])
    arr_nighttime = arr_daytime.copy()
    daylight_hrs = arr_daytime.copy().astype(np.int16)
    # Do sunrise and sunset define an interval? (Sunset > Sunrise)?
    inside_interval = np.apply_along_axis(lambda x: x[1] > x[0], 0, updown)
    # Or is the sun never up?
    never_up = np.logical_and(updown[0,...] == -1, updown[1,...] == -1)
    # Iteratively sum daytime VPD and temperature values
    for hr in range(0, 24):
        # Given only hour of sunrise/set on a 24-hour clock...
        #   if sun rises and sets on same day: SUNRISE <= HOUR <= SUNSET;
        #   if sun sets on next day: either SUNRISE <= HOUR or HOUR <= SUNSET;
        sun_is_up = np.logical_or( # Either...
            np.logical_and(inside_interval, # ...Rises and sets same day
                np.logical_and(updown[0,...] <= hr, hr <= updown[1,...])),
            np.logical_and(~inside_interval, # ...Sets on next day
                np.logical_or(updown[0,...] <= hr, hr <= updown[1,...])))
        # For simplicity, compute a 24-hour mean even if the sun never rises;
        #   there's no way to know what the "correct" daytime value is
        mask = np.logical_or(never_up, sun_is_up)
        np.add(np.where(
            mask, arr_24hr[hr,...], 0), arr_daytime, out = arr_daytime)
        np.add(np.where(
            ~mask, arr_24hr[hr,...], 0), arr_nighttime, out = arr_nighttime)
        # Keep track of the denominator (hours) for calculating the mean;
        #   note that this over-estimates actual daylight hours by 1 hour
        #   but results in the correct denominator for the sums above
        np.add(np.where(mask, 1, 0), daylight_hrs, out = daylight_hrs)
    arr_24hr = None
    # Calculate mean quantities
    if reducer == 'mean':
        arr_daytime = np.divide(arr_daytime, daylight_hrs)
        arr_nighttime = np.divide(arr_nighttime, 24 - daylight_hrs)
        # For sites where the sun is always above/ below the horizon, set missing
        #   nighttime values to zero
        arr_nighttime[~np.isfinite(arr_nighttime)] = missing[1]
    return np.stack((arr_daytime, arr_nighttime))

def degree_lengths(phi, a=6378137, b=6356752.3142)

Returns the approximate length of degrees of latitude and longitude.

Source

https://en.wikipedia.org/wiki/Latitude

phi : Number Latitude, in degrees a : Number Radius of the Earth (major axis) in meters b : Number Length of minor axis of the Earth in meters

Returns

tuple

Length of a degree of (longitude, latitude), respectively

Expand source code

def degree_lengths(phi, a = 6378137, b = 6356752.3142):
    '''
    Returns the approximate length of degrees of latitude and longitude.
    Source:

        https://en.wikipedia.org/wiki/Latitude

    phi : Number
        Latitude, in degrees
    a : Number
        Radius of the Earth (major axis) in meters
    b : Number
        Length of minor axis of the Earth in meters

    Returns
    -------
    tuple
        Length of a degree of (longitude, latitude), respectively
    '''
    e2 = ((a**2) - (b**2)) / (a**2)
    # Approximate length of a degree of latitude
    lat_m = 111132.954 - (559.822 * np.cos(2 * np.deg2rad(phi))) +\
        (1.175 * np.cos(4 * np.deg2rad(phi)))
    lng_m = (np.pi * a * np.cos(np.deg2rad(phi))) / (
        180 * np.sqrt(1 - (e2 * np.sin(np.deg2rad(phi))**2)))
    return (lng_m, lat_m)

def nash_sutcliffe(predicted: Sequence[+T_co], observed: Sequence[+T_co], norm: bool = False) ‑> float

Computes the Nash-Sutcliffe efficiency of a model's predictions.

$$\mathrm{NSE} = 1 - rac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - ar{y}_i)^2}$$

Parameters

predicted : Sequence

Predicted values

observed : Sequence

Observed values

norm : bool

True to return the normalized Nash-Sutcliffe efficiency (Default: False)

Returns

float

 

Expand source code

def nash_sutcliffe(
        predicted: Sequence, observed: Sequence, norm: bool = False) -> float:
    '''
    Computes the Nash-Sutcliffe efficiency of a model's predictions.

    $$
    \mathrm{NSE} = 1 -
        \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \bar{y}_i)^2}
    $$

    Parameters
    ----------
    predicted : Sequence
        Predicted values
    observed : Sequence
        Observed values
    norm : bool
        True to return the normalized Nash-Sutcliffe efficiency
        (Default: False)

    Returns
    -------
    float
    '''
    nse = 1 - (np.nansum(np.power(observed - predicted, 2)) /\
        np.nansum(np.power(observed - np.nanmean(observed), 2)))
    if norm:
        return 1 / (2 - nse) # Return normalized NSE
    return nse

def ordinals365(dates: Sequence[+T_co]) ‑> Sequence[+T_co]

Returns a length-T sequence of ordinals on [1,365]. Can be used for indexing a 365-day climatology; see climatology365().

Parameters

dates : list or tuple

Sequence of datetime.datetime or datetime.date instances

Returns

list

 

Expand source code

def ordinals365(dates: Sequence) -> Sequence:
    '''
    Returns a length-T sequence of ordinals on [1,365]. Can be used for
    indexing a 365-day climatology; see `climatology365()`.

    Parameters
    ----------
    dates : list or tuple
        Sequence of datetime.datetime or datetime.date instances

    Returns
    -------
    list
    '''
    return [
        t - 1 if (year % 4 == 0 and t >= 60) else t
        for t, year in [(int(t.strftime('%j')), t.year) for t in dates]
    ]