Source code for ocgis.calc.library.index.freeze_thaw

import datetime as dt

import numpy as np
from ocgis import env
from ocgis.calc import base
from ocgis.util.units import get_are_units_equal_by_string_or_cfunits


[docs]class FreezeThaw(base.AbstractUnivariateSetFunction, base.AbstractParameterizedFunction): key = 'freezethaw' description = "Number of freeze-thaw events, where freezing and thawing occurs once a threshold of degree days below or above 0C is reached. A complete cycle (freeze-thaw-freeze) will return a value of 2. " long_name = "Number of freeze-thaw events" standard_name = "freeze-thaw" required_units = ['K', 'C'] parms_definition = {'threshold': float}
[docs] def calculate(self, values, threshold=15): """ Return the number of freeze-thaw transitions. A value of 2 corresponds to a complete cycle (frozen-thawed-frozen). :param threshold: The number of degree-days above or below the freezing point after which the ground is considered frozen or thawed. """ assert (len(values.shape) == 3) # Check temporal resolution t = self.field['time'].get_value()[:2] step = t[1] - t[0] if type(step) == dt.timedelta: assert step == dt.timedelta(days=1) else: if self.field['time'].units.startswith('days'): assert step == 1 else: raise NotImplementedError # Unit conversion units = self.field.data_variables[0].units if get_are_units_equal_by_string_or_cfunits(units, 'C', try_cfunits=env.USE_CFUNITS): tas = values elif get_are_units_equal_by_string_or_cfunits(units, 'K', try_cfunits=env.USE_CFUNITS): tas = values - 273.15 out = np.apply_along_axis(freezethaw1d, 0, tas, threshold=threshold) return np.ma.masked_invalid(out)
def freezethaw1d(x, threshold): """ Return the number of freeze-thaw transitions. Parameters ---------- x : ndarray The daily temperature series (C). threshold : float The threshold in degree-days above or below the freezing point at which we consider the soil thawed or frozen. Returns ------- out : int The number of times the medium thawed or froze, so that a complete freeze-thaw cycle corresponds to a count of 2. """ # Masked values are just compressed. if isinstance(x, np.ma.MaskedArray): if x.mask.all(): return np.nan else: x = x.compressed() # This avoids issues when the threshold is reached right at the first value. x = np.concatenate([[0, ], x]) # Compute the cumulative degree days relative to the freezing point. cx = np.cumsum(x) # Find the places where the temperature crosses the freezing point (FP). over = x >= 0 cross = np.concatenate([[0, ], np.nonzero(np.diff(over) != 0)[0]]) cycles = [0, ] for ci in cross: # Skip FP crossing if it occurs before the threshold is reached. if ci < np.abs(cycles[-1]): continue # Otherwise reset the cumulative sum starting from the crossing. d = cx[ci:] - cx[ci] # Find the first place where t is exceeded (from above or below) w = np.nonzero(np.abs(d) >= threshold)[0] if len(w) > 0: w = w[0] # <-- This is where it occurs. s = np.sign(d[w]) # <-- This indicates if its thawing or freezing. # Test for the alternance of freeze and thaw events. # Store only an event if its different from the last. if s != np.sign(cycles[-1]): cycles.append(s * (w + ci)) # Return the number of transitions from frozen to thawed or vice-versa return float(len(cycles) - 2) # There are two "artificial" transitions