Merge pull request #191 from USDA-ARS-NWRC/susong1999

Snow - Fix early return in Susong1999 model.

.travis.yml

@@ -20,6 +20,7 @@ env:
     - CIBW_MANYLINUX_X86_64_IMAGE=manylinux2014
     - CIBW_TEST_REQUIRES=nose
     - CIBW_TEST_COMMAND="NOT_ON_GOLD_HOST=yup nosetests -vv --exe smrf"
+    - CIBW_TEST_EXTRAS=tests
     - CIBW_BUILD="cp3*-manylinux_x86_64 cp3*-macosx_x86_64"
     - CIBW_SKIP="?p27* pp* cp35*"
     - CIBW_BUILD_VERBOSITY=1
@@ -45,7 +46,7 @@ unittest: &unittest
       packages:
         - eccodes
         - gcc@7
-  install: 
+  install:
     - python3 --version
     - python3 -m pip install tox-travis cython
   script: tox
@@ -99,7 +100,7 @@ jobs:
       python: 3.6
       install: skip
       script: travis/build-sdist.sh
-          
+
     # Deploy on linux
     - <<: *ci-build-wheels
       name: Build and deploy Linux wheels

requirements_dev.txt

@@ -10,4 +10,4 @@ pydata-sphinx-theme
 sphinxcontrib-bibtex>=1.0
 sphinxcontrib-websupport>=1.0.1
 pygit2
-git+https://github.com/USDA-ARS-NWRC/goldmeister.git@main
+git+https://github.com/USDA-ARS-NWRC/goldmeister.git@main

setup.cfg

@@ -1,2 +1,5 @@
 [flake8]
-max-line-length = 79
+max-line-length = 79
+per-file-ignores =
+    # imported but unused
+    __init__.py: F401

setup.py

@@ -142,6 +142,9 @@ def finalize_options(self):
             'sphinxcontrib-bibtex',
             'sphinxcontrib-websupport',
             'pydata-sphinx-theme'
+        ],
+        'tests': [
+            'mock',
         ]
     },
     use_scm_version={

smrf/distribute/precipitation.py

@@ -3,7 +3,7 @@
 import numpy as np
 
 from smrf.distribute import image_data
-from smrf.envphys import precip, snow, storms
+from smrf.envphys import precip, Snow, storms
 from smrf.utils import utils
 
 
@@ -368,7 +368,7 @@ def distribute_for_marks2017(self, data, precip_temp, ta, time):
                 self._logger.debug('''Calculating new snow density for
                                     storm #{0}'''.format(self.storm_id+1))
                 # determine the precip phase and den
-                snow_den, perc_snow = snow.calc_phase_and_density(
+                snow_den, perc_snow = Snow.phase_and_density(
                     precip_temp,
                     self.precip,
                     nasde_model=self.nasde_model)
@@ -417,7 +417,7 @@ def distribute_for_susong1999(self, data, ppt_temp, time):
             self.precip = utils.set_min_max(self.precip, self.min, self.max)
 
             # determine the precip phase and den
-            snow_den, perc_snow = snow.calc_phase_and_density(
+            snow_den, perc_snow = Snow.phase_and_density(
                 ppt_temp,
                 self.precip,
                 nasde_model=self.nasde_model)

smrf/envphys/__init__.py

@@ -0,0 +1 @@
+from .snow import Snow

smrf/envphys/nasde_model/__init__.py

@@ -0,0 +1,3 @@
+from .marks_2017 import Marks2017
+from .piecewise_suosong_1999 import PiecewiseSusong1999
+from .susong_1999 import Susong1999

smrf/envphys/nasde_model/marks_2017.py

@@ -0,0 +1,207 @@
+import numpy as np
+
+from .utils import check_temperature
+from .piecewise_suosong_1999 import PiecewiseSusong1999
+
+
+class Marks2017:
+    """
+    A new accumulated snow density model that accounts for compaction. The
+    model builds upon :func:`~smrf.envphys.snow.piecewise_susong1999` by
+    adding effects from compaction. Of four mechanisms for compaction, this
+    model accounts for compaction by destructive metamorphism and overburden.
+    These two processes are accounted for by calculating a proportionality
+    using data from Kojima, Yosida and Mellor. The overburden is in part
+    estimated using total storm deposition, where storms are defined in
+    :func:`~smrf.envphys.storms.tracking_by_station`. Once this is determined
+    the final snow density is applied through the entire storm only varying
+    with hourly temperature.
+
+    Snow Density:
+        .. math::
+            \\rho_{s} = \\rho_{ns} + (\\Delta \\rho_{c} + \\Delta
+            \\rho_{m}) \\rho_{ns}
+
+    Overburden Proportionality:
+        .. math::
+            \\Delta \\rho_{c} = 0.026 e^{-0.08 (T_{z} - T_{snow})}
+            SWE*  e^{-21.0 \\rho_{ns}}
+
+    Metamorphism Proportionality:
+        .. math::
+            \\Delta \\rho_{m} = 0.01 c_{11} e^{-0.04 (T_{z} - T_{snow})}
+
+            c_{11} = c_min + (T_{z} - T_{snow}) C_{factor} + 1.0
+
+    Constants:
+        .. math::
+            C_{factor} = 0.0013
+
+            Tz = 0.0
+
+            ex_{max} = 1.75
+
+            exr = 0.75
+
+            ex_{min} = 1.0
+
+            c1r = 0.043
+
+            c_{min} = 0.0067
+
+            c_{fac} = 0.0013
+
+            T_{min} = -10.0
+
+            T_{max} = 0.0
+
+            T_{z} = 0.0
+
+            T_{r0} = 0.5
+
+            P_{cr0} = 0.25
+
+            P_{c0} = 0.75
+    """
+
+    # These should be used but unsure why not
+    # ex_max = 1.75
+    # exr = 0.75
+    # ex_min = 1.0
+    # #c1_min = 0.026
+    # #c1_max = 0.069
+    # c1r = 0.043
+    # c_min = 0.0067
+    # cfac = 0.0013
+    T_MIN = -10.0
+    T_MAX = 0.0
+    TZ = 0.0
+    # Tr0 = 0.5
+    # Pcr0 = 0.25
+    # Pc0 = 0.75
+
+    RHO_WATER = 1000.0
+
+    @staticmethod
+    def run(tpp, pp):
+        """
+        Args:
+            tpp: Numpy array of a single hour of temperature, use dew point if
+                available [degrees C].
+            pp: Numpy array representing the total amount of precipitation
+                deposited during a storm in millimeters
+
+        Returns:
+            dictionary:
+              - **rho_s** (*numpy.array*) - Density of the fresh snow in
+                  kg/m^3.
+              - **swe** (*numpy.array*) - Snow water equivalent.
+              - **pcs** (*numpy.array*) - Percent of the precipitation that is
+                  snow in values 0.0-1.0.
+              - **rho_ns** (*numpy.array*) - Density of the uncompacted snow,
+                  which is equivalent to the output from
+                  :func:`~smrf.envphys.snow.piecewise_susong1999`.
+              - **d_rho_c** (*numpy.array*) - Proportional coefficient for
+                  compaction from overburden.
+              - **d_rho_m** (*numpy.array*) - Proportional coefficient for
+                  compaction from melt.
+              - **rho_s** (*numpy.array*) - Final density of the snow [kg/m^3].
+              - **rho** (*numpy.array*) - Density of the precipitation, which
+                  continuously ranges from low density snow to pure liquid
+                  water
+                  (50-1000 kg/m^3).
+              - **zs** (*numpy.array*) - Snow height added from the
+                  precipitation
+
+        """
+
+        rho_ns = d_rho_m = d_rho_c = zs = rho_s = swe = pcs = np.zeros(
+            tpp.shape
+        )
+        rho = np.ones(tpp.shape)
+
+        if np.any(pp > 0):
+
+            # check the precipitation temperature
+            tpp, tsnow = check_temperature(
+                tpp, t_max=Marks2017.T_MAX, t_min=Marks2017.T_MIN
+            )
+
+            # Calculate the percent snow and initial uncompacted density
+            result = PiecewiseSusong1999.run(
+                tpp, pp, t_max=Marks2017.T_MAX, t_min=Marks2017.T_MIN
+            )
+            pcs = result['pcs']
+            rho_orig = result['rho_s']
+
+            swe = pp * pcs
+
+            # Calculate the density only where there is swe
+            swe_ind = swe > 0.0
+            if np.any(swe_ind):
+
+                s_pcs = pcs[swe_ind]
+                s_pp = pp[swe_ind]
+                s_swe = swe[swe_ind]
+                # transforms to a 1D array, will put back later
+                s_tsnow = tsnow[swe_ind]
+
+                # transforms to a 1D array, will put back later
+                s_rho_ns = rho_orig[swe_ind]
+
+                # Convert to a percentage of water
+                s_rho_ns = s_rho_ns/Marks2017.RHO_WATER
+
+                # proportional total storm mass compaction
+                s_d_rho_c = (0.026 * np.exp(-0.08 * (Marks2017.TZ - s_tsnow))
+                             * s_swe * np.exp(-21.0 * s_rho_ns))
+
+                ind = (s_rho_ns * Marks2017.RHO_WATER) >= 100.0
+                c11 = np.ones(s_rho_ns.shape)
+
+                # c11[ind] = (c_min + ((Tz - s_tsnow[ind]) * cfac)) + 1.0
+
+                c11[ind] = np.exp(-0.046*(s_rho_ns[ind]*1000.0-100.0))
+
+                s_d_rho_m = 0.01 * c11 * np.exp(
+                    -0.04 * (Marks2017.TZ - s_tsnow)
+                )
+
+                # Snow density, depth & combined liquid and snow density
+                s_rho_s = s_rho_ns + ((s_d_rho_c + s_d_rho_m) * s_rho_ns)
+
+                s_zs = s_swe / s_rho_s
+
+                # a mixture of snow and liquid
+                s_rho = s_rho_s.copy()
+                mix = (s_swe < s_pp) & (s_pcs > 0.0)
+                if np.any(mix):
+                    s_rho[mix] = \
+                        (s_pcs[mix] * s_rho_s[mix]) + (1.0 - s_pcs[mix])
+
+                s_rho[s_rho > 1.0] = 1.0
+
+                # put the values back into the larger array
+                rho_ns[swe_ind] = s_rho_ns
+                d_rho_m[swe_ind] = s_d_rho_m
+                d_rho_c[swe_ind] = s_d_rho_c
+                zs[swe_ind] = s_zs
+                rho_s[swe_ind] = s_rho_s
+                rho[swe_ind] = s_rho
+                pcs[swe_ind] = s_pcs
+
+        # convert densities from proportions, to kg/m^3 or mm/m^2
+        rho_ns = rho_ns * Marks2017.RHO_WATER
+        rho_s = rho_s * Marks2017.RHO_WATER
+        rho = rho * Marks2017.RHO_WATER
+
+        return {
+            'swe': swe,
+            'pcs': pcs,
+            'rho_ns': rho_ns,
+            'd_rho_c': d_rho_c,
+            'd_rho_m': d_rho_m,
+            'rho_s': rho_s,
+            'rho': rho,
+            'zs': zs
+        }

smrf/envphys/nasde_model/piecewise_suosong_1999.py

@@ -0,0 +1,87 @@
+from .utils import calc_percent_snow, check_temperature
+
+
+class PiecewiseSusong1999:
+    """
+    Follows :func:`~smrf.envphys.snow.susong1999` but is the piecewise form
+    of table shown there. This model adds to the former by accounting for
+    liquid water effect near 0.0 Degrees C.
+
+    The table was estimated by Danny Marks in 2017 which resulted in the
+    piecewise equations below:
+
+        Percent Snow:
+            .. math::
+
+                \\%_{snow} = \\Bigg \\lbrace{
+                    \\frac{-T}{T_{r0}} P_{cr0} + P_{c0}, \\hfill
+                    -0.5^{\\circ} C \\leq T \\leq 0.0^{\\circ} C
+                     \\atop
+                     \\frac{-T_{pp}}{T_{max} + 1.0} P_{c0} + P_{c0},
+                     \\hfill 0.0^\\circ C \\leq T \\leq T_{max}
+                    }
+
+        Snow Density:
+            .. math::
+                \\rho_{s} = 50.0 + 1.7 * (T_{pp} + 15.0)^{ex}
+
+
+                ex = \\Bigg \\lbrace{
+                        ex_{min} + \\frac{T_{range} + T_{snow} - T_{max}}
+                        {T_{range}} * ex_{r}, \\hfill ex < 1.75
+                        \\atop
+                        1.75, \\hfill, ex > 1.75
+                        }
+
+    """
+    EX_MAX = 1.75
+    EXR = 0.75
+    EX_MIN = 1.0
+
+    TZ = 0.0
+
+    @staticmethod
+    def run(tpp, _precipitation, t_max=0.0, t_min=-10.0, check_temps=True):
+        """
+        Args:
+            tpp: A numpy array of temperature, use dew point temperature
+                if available [degree C].
+            _precipitation: A numpy array of precipitation in millimeters.
+                Currently unused.
+            t_max: Max temperature that snow density is modeled.
+                   Default is 0.0 Degrees C.
+            t_min: Minimum temperature that snow density is changing.
+                   Default is -10.0 Degrees C.
+            check_temps: A boolean value check to apply special temperature
+                constraints, this is done using
+                :func:`~smrf.envphys.Snow.check_temperature`. Default is True.
+
+        Returns:
+            Dictionary:
+               - **pcs** (*numpy.array*) - Percent of the precipitation that is
+                    snow in values 0.0-1.0.
+               - **rho_s** (*numpy.array*) - Density of the fresh snow
+                    in kg/m^3.
+        """
+        # again, this shouldn't be needed in this context
+        if check_temps:
+            tpp, tsnow = check_temperature(
+                tpp, t_max=t_max, t_min=t_min
+            )
+
+        pcs = calc_percent_snow(tpp, t_max=t_max)
+
+        # New snow density - no compaction
+        t_range = t_max - t_min
+        ex = PiecewiseSusong1999.EX_MIN + (
+            ((t_range + (tsnow - t_max)) / t_range) * PiecewiseSusong1999.EXR
+        )
+
+        ex[ex > PiecewiseSusong1999.EX_MAX] = PiecewiseSusong1999.EX_MAX
+
+        rho_ns = 50.0 + (1.7 * ((tpp - PiecewiseSusong1999.TZ) + 15.0)**ex)
+
+        return {
+            'pcs': pcs,
+            'rho_s': rho_ns
+        }