USDA-ARS-NWRC / smrf / 328

pending completion

Build # 328

Build Type

Pull #170

travis-ci-com

Committed by

web-flow

Commit Message

<a href="https://github.com/USDA-ARS-NWRC/smrf/commit/<a class=hub.com/USDA-ARS-NWRC/smrf/commit/d68ba8972be73946c9b6a30a35b935d930ab07d8">d68ba8972<a href="https://github.com/USDA-ARS-NWRC/smrf/commit/d68ba8972be73946c9b6a30a35b935d930ab07d8">">Merge </a><a class="double-link" href="https://github.com/USDA-ARS-NWRC/smrf/commit/<a class="double-link" href="https://github.com/USDA-ARS-NWRC/smrf/commit/893f30cce1f1b0a8178bf90984c4e333552160b8">893f30cce</a>">893f30cce</a><a href="https://github.com/USDA-ARS-NWRC/smrf/commit/d68ba8972be73946c9b6a30a35b935d930ab07d8"> into 4e1ebc5df">4e1ebc5df</a>

Pull Request Pull Request #170: SMRF package

Run Details

261 of 261 new or added lines in 51 files covered. (100.0%)

3414 of 4319 relevant lines covered (79.05%)

0.79 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

76.47

/smrf/distribute/precipitation.py


import logging

import netCDF4 as nc
import numpy as np

from smrf.distribute import image_data
from smrf.envphys import precip, snow, storms
from smrf.utils import utils


class ppt(image_data.image_data):
    """
    The :mod:`~smrf.distribute.precip.ppt` class allows for variable specific
    distributions that go beyond the base class.

    The instantaneous precipitation typically has a positive trend with
    elevation due to orographic effects. However, the precipitation
    distribution can be further complicated for storms that have isolated
    impact at only a few measurement locations, for example thunderstorms
    or small precipitation events. Some distribution methods may be better
    suited than others for capturing the trend of these small events with
    multiple stations that record no precipitation may have a negative impact
    on the distribution.

    The precipitation phase, or the amount of precipitation falling as rain or
    snow, can significantly alter the energy and mass balance of the snowpack,
    either leading to snow accumulation or inducing melt :cite:`Marks&al:1998`
    :cite:`Kormos&al:2014`. The precipitation phase and initial snow density
    estimated using a variety of models that can be set in the configuration
    file.

    For more information on the available models, checkout
    :mod:`~smrf.envphys.snow`.

    After the precipitation phase is calculated, the storm information can be
    determined. The spatial resolution for which storm definitions are applied
    is based on the snow model thats selected.

    The time since last storm is based on an accumulated precipitation mass
    threshold, the time elapsed since it last snowed, and the precipitation
    phase.  These factors determine the start and end time of a storm that
    has produced enough precipitation as snow to change the surface albedo.

    Args:
        pptConfig: The [precip] section of the configuration file
        time_step: The time step in minutes of the data, defaults to 60

    Attributes:
        config: configuration from [precip] section
        precip: numpy array of the precipitation
        percent_snow: numpy array of the percent of time step that was snow
        snow_density: numpy array of the snow density
        storm_days: numpy array of the days since last storm
        storm_total: numpy array of the precipitation mass for the storm
        last_storm_day: numpy array of the day of the last storm (decimal day)
        last_storm_day_basin: maximum value of last_storm day within the mask
            if specified
        min: minimum value of precipitation is 0
        max: maximum value of precipitation is infinite
        stations: stations to be used in alphabetical order

    """

    variable = 'precip'

    # these are variables that can be output
    output_variables = {
        'precip': {
            'units': 'mm',
            'standard_name': 'precipitation_mass',
            'long_name': 'Precipitation mass'
        },
        'percent_snow': {
            'units': '%',
            'standard_name': 'percent_snow',
            'long_name': 'Percent of precipitation as snow'
        },
        'snow_density': {
            'units': 'kg/m3',
            'standard_name': 'snow_density',
            'long_name': 'Precipitation snow density'
        },
        'storm_days': {
            'units': 'day',
            'standard_name': 'days_since_last_storm',
            'long_name': 'Days since the last storm'
        },
        'storm_total': {
            'units': 'mm',
            'standard_name': 'precipitation_mass_storm',
            'long_name': 'Precipitation mass for the storm period'
        },
        'last_storm_day': {
            'units': 'day',
            'standard_name': 'day_of_last_storm',
            'long_name': 'Decimal day of the last storm since Oct 1'
        }
    }

    # these are variables that are operate at the end only and do not need to
    # be written during main distribute loop
    post_process_variables = {}

    def __init__(self, pptConfig, start_date, time_step=60):

        # extend the base class
        image_data.image_data.__init__(self, self.variable)
        self._logger = logging.getLogger(__name__)

        # check and assign the configuration
        self.getConfig(pptConfig)
        self.time_step = float(time_step)
        self.start_date = start_date

    def initialize(self, topo, data):

        self._logger.debug('Initializing distribute.precip')

        self._initialize(topo, data.metadata)
        self.percent_snow = np.zeros((topo.ny, topo.nx))
        self.snow_density = np.zeros((topo.ny, topo.nx))
        self.storm_days = np.zeros((topo.ny, topo.nx))
        self.storm_total = np.zeros((topo.ny, topo.nx))
        self.last_storm_day = np.zeros((topo.ny, topo.nx))
        self.dem = topo.dem

        # Assign storm_days array if given
        if self.config["storm_days_restart"] is not None:
            self._logger.debug('Reading {} from {}'.format(
                'storm_days', self.config['storm_days_restart']))
            f = nc.Dataset(self.config['storm_days_restart'], 'r')

            if 'storm_days' in f.variables:
                time = f.variables['time'][:]
                t = f.variables['time']
                time = nc.num2date(t[:], t.getncattr(
                    'units'), t.getncattr('calendar'))

                # start at index of storm_days - 1
                time_ind = np.where(time == self.start_date)[0] - 1

                if not time_ind:
                    self._logger.warning(
                        'Invalid storm_days input! Setting to 0.0')
                    self.storm_days = np.zeros((topo.ny, topo.nx))

                else:
                    self.storm_days = f.variables['storm_days'][time_ind, :, :][0]  # noqa
            else:
                self._logger.warning(
                    'Variable {} not in {}, setting to 0.0'.format(
                        'storm_days', self.config['storm_days_restart']))
                self.storm_days = np.zeros((topo.ny, topo.nx))

            f.close()

        self.ppt_threshold = self.config['storm_mass_threshold']

        # Time steps needed to end a storm definition
        self.time_to_end_storm = self.config['marks2017_timesteps_to_end_storms']  # noqa

        self.nasde_model = self.config['new_snow_density_model']

        self._logger.info(
            """Using {0} for the new accumulated snow density"""
            """ model:  """.format(self.nasde_model))

        if self.nasde_model == 'marks2017':

            self.storm_total = np.zeros((topo.ny, topo.nx))

            self.storms = []
            self.time_steps_since_precip = 0
            self.storming = False

            # Clip and adjust the precip data so that there is only precip
            # during the storm and ad back in the missing data to conserve mass
            data.precip = data.precip[self.stations]
            self.storms, storm_count = storms.tracking_by_station(
                data.precip,
                mass_thresh=self.ppt_threshold,
                steps_thresh=self.time_to_end_storm)
            self.corrected_precip = storms.clip_and_correct(
                data.precip,
                self.storms,
                stations=self.stations)

            if storm_count != 0:
                self._logger.info(
                    "Identified Storms:\n{0}".format(self.storms))
                self.storm_id = 0
                self._logger.info(
                    "Estimated number of storms: {0}".format(storm_count))

            else:
                if (data.precip.sum() > 0).any():
                    self.storm_id = np.nan
                    self._logger.warning("Zero events triggered a storm "
                                         " definition, None of the precip will"
                                         " be used in this run.")

        # if redistributing due to wind
        if self.config['precip_rescaling_model'] == 'winstral':
            self._tbreak_file = nc.Dataset(
                self.config['winstral_tbreak_netcdf'], 'r')
            self.tbreak = self._tbreak_file.variables['tbreak'][:]
            self.tbreak_direction = self._tbreak_file.variables['direction'][:]
            self._tbreak_file.close()
            self._logger.debug('Read data from {}'
                               .format(self.config['winstral_tbreak_netcdf']))

            # get the veg values
            self.veg_type = topo.veg_type
            matching = [s for s in self.config.keys() if "winstral_veg_" in s]
            v = {}
            for m in matching:
                if m != 'winstral_veg_default':
                    ms = m.split('_')
                    # v[ms[1]] = float(self.config[m])
                    if type(self.config[m]) == list:
                        v[ms[1]] = float(self.config[m][0])
                    else:
                        v[ms[1]] = float(self.config[m])
            self.veg = v

        self.mask = topo.mask

    def distribute(self, data, dpt, precip_temp, ta, time, wind, temp, az,
                   dir_round_cell, wind_speed, cell_maxus, mask=None):
        """
        Distribute given a Panda's dataframe for a single time step. Calls
        :mod:`smrf.distribute.image_data.image_data._distribute`.

        The following steps are taken when distributing precip, if there is
        precipitation measured:

        1. Distribute the instantaneous precipitation from the measurement data
        2. Determine the distributed precipitation phase based on the
            precipitation temperature
        3. Calculate the storms based on the accumulated mass, time since last
            storm, and precipitation phase threshold


        Args:
            data:           Pandas dataframe for a single time step from precip
            dpt:            dew point numpy array that will be used for
            precip_temp:    numpy array of the precipitation temperature
            ta:             air temp numpy array
            time:           pass in the time were are currently on
            wind:           station wind speed at time step
            temp:           station air temperature at time step
            az:             numpy array for simulated wind direction
            dir_round_cell: numpy array for wind direction in discreet
                            increments for referencing maxus at a specific
                            direction
            wind_speed:     numpy array of wind speed
            cell_maxus:     numpy array for maxus at correct wind directions
            mask:           basin mask to apply to the storm days for
                            calculating the last storm day for the basin
        """

        self._logger.debug('%s Distributing all precip' % data.name)
        data = data[self.stations]

        if self.config['distribution'] != 'grid':
            if self.nasde_model == 'marks2017':
                # Adjust the precip for undercatchment
                if self.config['station_adjust_for_undercatch']:
                    self._logger.debug(
                        '%s Adjusting precip for undercatch...' % data.name)
                    self.corrected_precip.loc[time] = \
                        precip.adjust_for_undercatch(
                            self.corrected_precip.loc[time],
                        wind,
                        temp,
                        self.config,
                        self.metadata)

                # Use the clipped and corrected precip
                self.distribute_for_marks2017(
                    self.corrected_precip.loc[time],
                    precip_temp,
                    ta,
                    time,
                    mask=mask)

            else:
                # Adjust the precip for undercatchment
                if self.config['station_adjust_for_undercatch']:
                    self._logger.debug(
                        '%s Adjusting precip for undercatch...' % data.name)
                    data = precip.adjust_for_undercatch(
                        data,
                        wind,
                        temp,
                        self.config,
                        self.metadata)

                    self.distribute_for_susong1999(
                        data, precip_temp, time, mask=mask)
        else:
            self.distribute_for_susong1999(data, precip_temp, time, mask=mask)

        # redistribute due to wind to account for driftin
        if self.config['precip_rescaling_model'] == 'winstral':
            self._logger.debug('%s Redistributing due to wind' % data.name)
            if np.any(dpt < 0.5):
                self.precip = precip.dist_precip_wind(
                    self.precip,
                    dpt,
                    az,
                    dir_round_cell,
                    wind_speed,
                    cell_maxus,
                    self.tbreak,
                    self.tbreak_direction,
                    self.veg_type,
                    self.veg,
                    self.config)

    def distribute_for_marks2017(self, data, precip_temp, ta, time, mask=None):
        """
        Specialized distribute function for working with the new accumulated
        snow density model Marks2017 requires storm total and a corrected
        precipitation as to avoid precip between storms.
        """
        # self.corrected_precip # = data.mul(self.storm_correction)

        if data.sum() > 0.0:
            # Check for time in every storm
            for i, s in self.storms.iterrows():
                storm_start = s['start']
                storm_end = s['end']

                if time >= storm_start and time <= storm_end:
                    # establish storm info
                    self.storm_id = i
                    storm = self.storms.iloc[self.storm_id]
                    self.storming = True
                    break
                else:
                    self.storming = False

            self._logger.debug("Storming? {0}".format(self.storming))
            self._logger.debug("Current Storm ID = {0}".format(self.storm_id))

            # distribute data and set the min/max
            self._distribute(data, zeros=None)
            self.precip = utils.set_min_max(self.precip, self.min, self.max)

            if time == storm_start:
                # Entered into a new storm period distribute the storm total
                self._logger.debug('{0} Entering storm #{1}'
                                   .format(data.name, self.storm_id+1))
                if precip_temp.min() < 2.0:
                    self._logger.debug('''Distributing Total Precip
                                        for Storm #{0}'''
                                       .format(self.storm_id+1))
                    self._distribute(storm[self.stations].astype(float),
                                     other_attribute='storm_total')
                    self.storm_total = utils.set_min_max(self.storm_total,
                                                         self.min,
                                                         self.max)

            if self.storming and precip_temp.min() < 2.0:
                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(
                    precip_temp,
                    self.precip,
                    nasde_model=self.nasde_model)

            else:
                snow_den = np.zeros(self.precip.shape)
                perc_snow = np.zeros(self.precip.shape)

            # calculate decimal days since last storm
            self.storm_days = storms.time_since_storm_pixel(
                self.precip,
                precip_temp,
                perc_snow,
                storming=self.storming,
                time_step=self.time_step/60.0/24.0,
                stormDays=self.storm_days,
                mass=self.ppt_threshold)

        else:
            self.storm_days += self.time_step/60.0/24.0
            self.precip = np.zeros(self.storm_days.shape)
            perc_snow = np.zeros(self.storm_days.shape)
            snow_den = np.zeros(self.storm_days.shape)

        # save the model state
        self.percent_snow = perc_snow
        self.snow_density = snow_den

        # day of last storm, this will be used in albedo
        self.last_storm_day = utils.water_day(data.name)[0] - \
            self.storm_days - 0.001

        # get the time since most recent storm
        if mask is not None:
            self.last_storm_day_basin = np.max(mask * self.last_storm_day)
        else:
            self.last_storm_day_basin = np.max(self.last_storm_day)

    def distribute_for_susong1999(self, data, ppt_temp, time, mask=None):
        """Susong 1999 estimates percent snow and snow density based on
        Susong et al, (1999) :cite:`Susong&al:1999`.

        Args:
            data (pd.DataFrame): Precipitation mass data
            ppt_temp (pd.DataFrame): Precipitation temperature data
            time : Unused
            mask (np.array, optional): Mask the output. Defaults to None.
        """

        if data.sum() > 0:

            self._distribute(data)
            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(
                ppt_temp,
                self.precip,
                nasde_model=self.nasde_model)

            # determine the time since last storm
            stormDays, stormPrecip = storms.time_since_storm(
                self.precip,
                perc_snow,
                time_step=self.time_step/60/24,
                mass=self.ppt_threshold,
                time=self.time_to_end_storm,
                stormDays=self.storm_days,
                stormPrecip=self.storm_total)

            # save the model state
            self.percent_snow = perc_snow
            self.snow_density = snow_den
            self.storm_days = stormDays
            self.storm_total = stormPrecip

        else:

            self.storm_days += self.time_step/60/24

            # make everything else zeros
            self.precip = np.zeros(self.storm_days.shape)
            self.percent_snow = np.zeros(self.storm_days.shape)
            self.snow_density = np.zeros(self.storm_days.shape)

        # day of last storm, this will be used in albedo
        self.last_storm_day = utils.water_day(data.name)[0] - \
            self.storm_days - 0.001

        # get the time since most recent storm
        if mask is not None:
            self.last_storm_day_basin = np.max(mask * self.last_storm_day)
        else:
            self.last_storm_day_basin = np.max(self.last_storm_day)

    def distribute_thread(self, queue, data, date, mask=None):
        """
        Distribute the data using threading and queue. All data is provided and
        ``distribute_thread`` will go through each time step and call
        :mod:`smrf.distribute.precip.ppt.distribute` then puts the distributed
        data into the queue for:

        * :py:attr:`percent_snow`

        * :py:attr:`snow_density`
        * :py:attr:`storm_days`
        * :py:attr:`last_storm_day_basin`

        Args:
            queue: queue dictionary for all variables
            data: pandas dataframe for all data, indexed by date time
        """
        self._logger.info("Distributing {}".format(self.variable))

        for t in data.precip.index:

            dpt = queue['dew_point'].get(t)
            precip_temp = queue['precip_temp'].get(t)
            ta = queue['air_temp'].get(t)

            # variables for wind redistribution
            az = queue['wind_direction'].get(t)
            flatwind = queue['flatwind'].get(t)
            dir_round_cell = queue['dir_round_cell'].get(t)
            cell_maxus = queue['cellmaxus'].get(t)

            self.distribute(data.precip.loc[t], dpt, precip_temp, ta, t,
                            data.wind_speed.loc[t], data.air_temp.loc[t],
                            az, dir_round_cell, flatwind, cell_maxus,
                            mask=mask)

            queue[self.variable].put([t, self.precip])
            queue['percent_snow'].put([t, self.percent_snow])
            queue['snow_density'].put([t, self.snow_density])
            queue['last_storm_day_basin'].put([t, self.last_storm_day_basin])
            queue['storm_days'].put([t, self.storm_days])
            queue['storm_total'].put([t, self.storm_total])
            if self.nasde_model == "marks2017":
                queue['storm_id'].put([t, self.storm_id])

    def post_processor(self, main_obj, threaded=False):

        pass

1
2	import logging	1✔
3
4	import netCDF4 as nc	1✔
5	import numpy as np	1✔
6
7	from smrf.distribute import image_data	1✔
8	from smrf.envphys import precip, snow, storms	1✔
9	from smrf.utils import utils	1✔
10
11
12	class ppt(image_data.image_data):	1✔
13	"""
14	The :mod:`~smrf.distribute.precip.ppt` class allows for variable specific
15	distributions that go beyond the base class.
16
17	The instantaneous precipitation typically has a positive trend with
18	elevation due to orographic effects. However, the precipitation
19	distribution can be further complicated for storms that have isolated
20	impact at only a few measurement locations, for example thunderstorms
21	or small precipitation events. Some distribution methods may be better
22	suited than others for capturing the trend of these small events with
23	multiple stations that record no precipitation may have a negative impact
24	on the distribution.
25
26	The precipitation phase, or the amount of precipitation falling as rain or
27	snow, can significantly alter the energy and mass balance of the snowpack,
28	either leading to snow accumulation or inducing melt :cite:`Marks&al:1998`
29	:cite:`Kormos&al:2014`. The precipitation phase and initial snow density
30	estimated using a variety of models that can be set in the configuration
31	file.
32
33	For more information on the available models, checkout
34	:mod:`~smrf.envphys.snow`.
35
36	After the precipitation phase is calculated, the storm information can be
37	determined. The spatial resolution for which storm definitions are applied
38	is based on the snow model thats selected.
39
40	The time since last storm is based on an accumulated precipitation mass
41	threshold, the time elapsed since it last snowed, and the precipitation
42	phase. These factors determine the start and end time of a storm that
43	has produced enough precipitation as snow to change the surface albedo.
44
45	Args:
46	pptConfig: The [precip] section of the configuration file
47	time_step: The time step in minutes of the data, defaults to 60
48
49	Attributes:
50	config: configuration from [precip] section
51	precip: numpy array of the precipitation
52	percent_snow: numpy array of the percent of time step that was snow
53	snow_density: numpy array of the snow density
54	storm_days: numpy array of the days since last storm
55	storm_total: numpy array of the precipitation mass for the storm
56	last_storm_day: numpy array of the day of the last storm (decimal day)
57	last_storm_day_basin: maximum value of last_storm day within the mask
58	if specified
59	min: minimum value of precipitation is 0
60	max: maximum value of precipitation is infinite
61	stations: stations to be used in alphabetical order
62
63	"""
64
65	variable = 'precip'	1✔
66
67	# these are variables that can be output
68	output_variables = {	1✔
69	'precip': {
70	'units': 'mm',
71	'standard_name': 'precipitation_mass',
72	'long_name': 'Precipitation mass'
73	},
74	'percent_snow': {
75	'units': '%',
76	'standard_name': 'percent_snow',
77	'long_name': 'Percent of precipitation as snow'
78	},
79	'snow_density': {
80	'units': 'kg/m3',
81	'standard_name': 'snow_density',
82	'long_name': 'Precipitation snow density'
83	},
84	'storm_days': {
85	'units': 'day',
86	'standard_name': 'days_since_last_storm',
87	'long_name': 'Days since the last storm'
88	},
89	'storm_total': {
90	'units': 'mm',
91	'standard_name': 'precipitation_mass_storm',
92	'long_name': 'Precipitation mass for the storm period'
93	},
94	'last_storm_day': {
95	'units': 'day',
96	'standard_name': 'day_of_last_storm',
97	'long_name': 'Decimal day of the last storm since Oct 1'
98	}
99	}
100
101	# these are variables that are operate at the end only and do not need to
102	# be written during main distribute loop
103	post_process_variables = {}	1✔
104
105	def __init__(self, pptConfig, start_date, time_step=60):	1✔
106
107	# extend the base class
108	image_data.image_data.__init__(self, self.variable)	1✔
109	self._logger = logging.getLogger(__name__)	1✔
110
111	# check and assign the configuration
112	self.getConfig(pptConfig)	1✔
113	self.time_step = float(time_step)	1✔
114	self.start_date = start_date	1✔
115
116	def initialize(self, topo, data):	1✔
117
118	self._logger.debug('Initializing distribute.precip')	1✔
119
120	self._initialize(topo, data.metadata)	1✔
121	self.percent_snow = np.zeros((topo.ny, topo.nx))	1✔
122	self.snow_density = np.zeros((topo.ny, topo.nx))	1✔
123	self.storm_days = np.zeros((topo.ny, topo.nx))	1✔
124	self.storm_total = np.zeros((topo.ny, topo.nx))	1✔
125	self.last_storm_day = np.zeros((topo.ny, topo.nx))	1✔
126	self.dem = topo.dem	1✔
127
128	# Assign storm_days array if given
129	if self.config["storm_days_restart"] is not None:	1✔
130	self._logger.debug('Reading {} from {}'.format(	×
131	'storm_days', self.config['storm_days_restart']))
132	f = nc.Dataset(self.config['storm_days_restart'], 'r')	×
133
134	if 'storm_days' in f.variables:	×
135	time = f.variables['time'][:]	×
136	t = f.variables['time']	×
137	time = nc.num2date(t[:], t.getncattr(	×
138	'units'), t.getncattr('calendar'))
139
140	# start at index of storm_days - 1
141	time_ind = np.where(time == self.start_date)[0] - 1	×
142
143	if not time_ind:	×
144	self._logger.warning(	×
145	'Invalid storm_days input! Setting to 0.0')
146	self.storm_days = np.zeros((topo.ny, topo.nx))	×
147
148	else:
149	self.storm_days = f.variables['storm_days'][time_ind, :, :][0] # noqa	×
150	else:
151	self._logger.warning(	×
152	'Variable {} not in {}, setting to 0.0'.format(
153	'storm_days', self.config['storm_days_restart']))
154	self.storm_days = np.zeros((topo.ny, topo.nx))	×
155
156	f.close()	×
157
158	self.ppt_threshold = self.config['storm_mass_threshold']	1✔
159
160	# Time steps needed to end a storm definition
161	self.time_to_end_storm = self.config['marks2017_timesteps_to_end_storms'] # noqa	1✔
162
163	self.nasde_model = self.config['new_snow_density_model']	1✔
164
165	self._logger.info(	1✔
166	"""Using {0} for the new accumulated snow density"""
167	""" model: """.format(self.nasde_model))
168
169	if self.nasde_model == 'marks2017':	1✔
170
171	self.storm_total = np.zeros((topo.ny, topo.nx))	1✔
172
173	self.storms = []	1✔
174	self.time_steps_since_precip = 0	1✔
175	self.storming = False	1✔
176
177	# Clip and adjust the precip data so that there is only precip
178	# during the storm and ad back in the missing data to conserve mass
179	data.precip = data.precip[self.stations]	1✔
180	self.storms, storm_count = storms.tracking_by_station(	1✔
181	data.precip,
182	mass_thresh=self.ppt_threshold,
183	steps_thresh=self.time_to_end_storm)
184	self.corrected_precip = storms.clip_and_correct(	1✔
185	data.precip,
186	self.storms,
187	stations=self.stations)
188
189	if storm_count != 0:	1✔
190	self._logger.info(	1✔
191	"Identified Storms:\n{0}".format(self.storms))
192	self.storm_id = 0	1✔
193	self._logger.info(	1✔
194	"Estimated number of storms: {0}".format(storm_count))
195
196	else:
197	if (data.precip.sum() > 0).any():	1✔
198	self.storm_id = np.nan	1✔
199	self._logger.warning("Zero events triggered a storm "	1✔
200	" definition, None of the precip will"
201	" be used in this run.")
202
203	# if redistributing due to wind
204	if self.config['precip_rescaling_model'] == 'winstral':	1✔
205	self._tbreak_file = nc.Dataset(	×
206	self.config['winstral_tbreak_netcdf'], 'r')
207	self.tbreak = self._tbreak_file.variables['tbreak'][:]	×
208	self.tbreak_direction = self._tbreak_file.variables['direction'][:]	×
209	self._tbreak_file.close()	×
210	self._logger.debug('Read data from {}'	×
211	.format(self.config['winstral_tbreak_netcdf']))
212
213	# get the veg values
214	self.veg_type = topo.veg_type	×
215	matching = [s for s in self.config.keys() if "winstral_veg_" in s]	×
216	v = {}	×
217	for m in matching:	×
218	if m != 'winstral_veg_default':	×
219	ms = m.split('_')	×
220	# v[ms[1]] = float(self.config[m])
221	if type(self.config[m]) == list:	×
222	v[ms[1]] = float(self.config[m][0])	×
223	else:
224	v[ms[1]] = float(self.config[m])	×
225	self.veg = v	×
226
227	self.mask = topo.mask	1✔
228
229	def distribute(self, data, dpt, precip_temp, ta, time, wind, temp, az,	1✔
230	dir_round_cell, wind_speed, cell_maxus, mask=None):
231	"""
232	Distribute given a Panda's dataframe for a single time step. Calls
233	:mod:`smrf.distribute.image_data.image_data._distribute`.
234
235	The following steps are taken when distributing precip, if there is
236	precipitation measured:
237
238	1. Distribute the instantaneous precipitation from the measurement data
239	2. Determine the distributed precipitation phase based on the
240	precipitation temperature
241	3. Calculate the storms based on the accumulated mass, time since last
242	storm, and precipitation phase threshold
243
244
245	Args:
246	data: Pandas dataframe for a single time step from precip
247	dpt: dew point numpy array that will be used for
248	precip_temp: numpy array of the precipitation temperature
249	ta: air temp numpy array
250	time: pass in the time were are currently on
251	wind: station wind speed at time step
252	temp: station air temperature at time step
253	az: numpy array for simulated wind direction
254	dir_round_cell: numpy array for wind direction in discreet
255	increments for referencing maxus at a specific
256	direction
257	wind_speed: numpy array of wind speed
258	cell_maxus: numpy array for maxus at correct wind directions
259	mask: basin mask to apply to the storm days for
260	calculating the last storm day for the basin
261	"""
262
263	self._logger.debug('%s Distributing all precip' % data.name)	1✔
264	data = data[self.stations]	1✔
265
266	if self.config['distribution'] != 'grid':	1✔
267	if self.nasde_model == 'marks2017':	1✔
268	# Adjust the precip for undercatchment
269	if self.config['station_adjust_for_undercatch']:	1✔
270	self._logger.debug(	1✔
271	'%s Adjusting precip for undercatch...' % data.name)
272	self.corrected_precip.loc[time] = \	1✔
273	precip.adjust_for_undercatch(
274	self.corrected_precip.loc[time],
275	wind,
276	temp,
277	self.config,
278	self.metadata)
279
280	# Use the clipped and corrected precip
281	self.distribute_for_marks2017(	1✔
282	self.corrected_precip.loc[time],
283	precip_temp,
284	ta,
285	time,
286	mask=mask)
287
288	else:
289	# Adjust the precip for undercatchment
290	if self.config['station_adjust_for_undercatch']:	×
291	self._logger.debug(	×
292	'%s Adjusting precip for undercatch...' % data.name)
293	data = precip.adjust_for_undercatch(	×
294	data,
295	wind,
296	temp,
297	self.config,
298	self.metadata)
299
300	self.distribute_for_susong1999(	×
301	data, precip_temp, time, mask=mask)
302	else:
303	self.distribute_for_susong1999(data, precip_temp, time, mask=mask)	1✔
304
305	# redistribute due to wind to account for driftin
306	if self.config['precip_rescaling_model'] == 'winstral':	1✔
307	self._logger.debug('%s Redistributing due to wind' % data.name)	×
308	if np.any(dpt < 0.5):	×
309	self.precip = precip.dist_precip_wind(	×
310	self.precip,
311	dpt,
312	az,
313	dir_round_cell,
314	wind_speed,
315	cell_maxus,
316	self.tbreak,
317	self.tbreak_direction,
318	self.veg_type,
319	self.veg,
320	self.config)
321
322	def distribute_for_marks2017(self, data, precip_temp, ta, time, mask=None):	1✔
323	"""
324	Specialized distribute function for working with the new accumulated
325	snow density model Marks2017 requires storm total and a corrected
326	precipitation as to avoid precip between storms.
327	"""
328	# self.corrected_precip # = data.mul(self.storm_correction)
329
330	if data.sum() > 0.0:	1✔
331	# Check for time in every storm
332	for i, s in self.storms.iterrows():	1✔
333	storm_start = s['start']	1✔
334	storm_end = s['end']	1✔
335
336	if time >= storm_start and time <= storm_end:	1✔
337	# establish storm info
338	self.storm_id = i	1✔
339	storm = self.storms.iloc[self.storm_id]	1✔
340	self.storming = True	1✔
341	break	1✔
342	else:
343	self.storming = False	×
344
345	self._logger.debug("Storming? {0}".format(self.storming))	1✔
346	self._logger.debug("Current Storm ID = {0}".format(self.storm_id))	1✔
347
348	# distribute data and set the min/max
349	self._distribute(data, zeros=None)	1✔
350	self.precip = utils.set_min_max(self.precip, self.min, self.max)	1✔
351
352	if time == storm_start:	1✔
353	# Entered into a new storm period distribute the storm total
354	self._logger.debug('{0} Entering storm #{1}'	1✔
355	.format(data.name, self.storm_id+1))
356	if precip_temp.min() < 2.0:	1✔
357	self._logger.debug('''Distributing Total Precip	1✔
358	for Storm #{0}'''
359	.format(self.storm_id+1))
360	self._distribute(storm[self.stations].astype(float),	1✔
361	other_attribute='storm_total')
362	self.storm_total = utils.set_min_max(self.storm_total,	1✔
363	self.min,
364	self.max)
365
366	if self.storming and precip_temp.min() < 2.0:	1✔
367	self._logger.debug('''Calculating new snow density for	1✔
368	storm #{0}'''.format(self.storm_id+1))
369	# determine the precip phase and den
370	snow_den, perc_snow = snow.calc_phase_and_density(	1✔
371	precip_temp,
372	self.precip,
373	nasde_model=self.nasde_model)
374
375	else:
376	snow_den = np.zeros(self.precip.shape)	×
377	perc_snow = np.zeros(self.precip.shape)	×
378
379	# calculate decimal days since last storm
380	self.storm_days = storms.time_since_storm_pixel(	1✔
381	self.precip,
382	precip_temp,
383	perc_snow,
384	storming=self.storming,
385	time_step=self.time_step/60.0/24.0,
386	stormDays=self.storm_days,
387	mass=self.ppt_threshold)
388
389	else:
390	self.storm_days += self.time_step/60.0/24.0	1✔
391	self.precip = np.zeros(self.storm_days.shape)	1✔
392	perc_snow = np.zeros(self.storm_days.shape)	1✔
393	snow_den = np.zeros(self.storm_days.shape)	1✔
394
395	# save the model state
396	self.percent_snow = perc_snow	1✔
397	self.snow_density = snow_den	1✔
398
399	# day of last storm, this will be used in albedo
400	self.last_storm_day = utils.water_day(data.name)[0] - \	1✔
401	self.storm_days - 0.001
402
403	# get the time since most recent storm
404	if mask is not None:	1✔
405	self.last_storm_day_basin = np.max(mask * self.last_storm_day)	1✔
406	else:
407	self.last_storm_day_basin = np.max(self.last_storm_day)	1✔
408
409	def distribute_for_susong1999(self, data, ppt_temp, time, mask=None):	1✔
410	"""Susong 1999 estimates percent snow and snow density based on
411	Susong et al, (1999) :cite:`Susong&al:1999`.
412
413	Args:
414	data (pd.DataFrame): Precipitation mass data
415	ppt_temp (pd.DataFrame): Precipitation temperature data
416	time : Unused
417	mask (np.array, optional): Mask the output. Defaults to None.
418	"""
419
420	if data.sum() > 0:	1✔
421
422	self._distribute(data)	1✔
423	self.precip = utils.set_min_max(self.precip, self.min, self.max)	1✔
424
425	# determine the precip phase and den
426	snow_den, perc_snow = snow.calc_phase_and_density(	1✔
427	ppt_temp,
428	self.precip,
429	nasde_model=self.nasde_model)
430
431	# determine the time since last storm
432	stormDays, stormPrecip = storms.time_since_storm(	1✔
433	self.precip,
434	perc_snow,
435	time_step=self.time_step/60/24,
436	mass=self.ppt_threshold,
437	time=self.time_to_end_storm,
438	stormDays=self.storm_days,
439	stormPrecip=self.storm_total)
440
441	# save the model state
442	self.percent_snow = perc_snow	1✔
443	self.snow_density = snow_den	1✔
444	self.storm_days = stormDays	1✔
445	self.storm_total = stormPrecip	1✔
446
447	else:
448
449	self.storm_days += self.time_step/60/24	1✔
450
451	# make everything else zeros
452	self.precip = np.zeros(self.storm_days.shape)	1✔
453	self.percent_snow = np.zeros(self.storm_days.shape)	1✔
454	self.snow_density = np.zeros(self.storm_days.shape)	1✔
455
456	# day of last storm, this will be used in albedo
457	self.last_storm_day = utils.water_day(data.name)[0] - \	1✔
458	self.storm_days - 0.001
459
460	# get the time since most recent storm
461	if mask is not None:	1✔
462	self.last_storm_day_basin = np.max(mask * self.last_storm_day)	×
463	else:
464	self.last_storm_day_basin = np.max(self.last_storm_day)	1✔
465
466	def distribute_thread(self, queue, data, date, mask=None):	1✔
467	"""
468	Distribute the data using threading and queue. All data is provided and
469	``distribute_thread`` will go through each time step and call
470	:mod:`smrf.distribute.precip.ppt.distribute` then puts the distributed
471	data into the queue for:
472
473	* :py:attr:`percent_snow`
474
475	* :py:attr:`snow_density`
476	* :py:attr:`storm_days`
477	* :py:attr:`last_storm_day_basin`
478
479	Args:
480	queue: queue dictionary for all variables
481	data: pandas dataframe for all data, indexed by date time
482	"""
483	self._logger.info("Distributing {}".format(self.variable))	1✔
484
485	for t in data.precip.index:	1✔
486
487	dpt = queue['dew_point'].get(t)	1✔
488	precip_temp = queue['precip_temp'].get(t)	1✔
489	ta = queue['air_temp'].get(t)	1✔
490
491	# variables for wind redistribution
492	az = queue['wind_direction'].get(t)	1✔
493	flatwind = queue['flatwind'].get(t)	1✔
494	dir_round_cell = queue['dir_round_cell'].get(t)	1✔
495	cell_maxus = queue['cellmaxus'].get(t)	1✔
496
497	self.distribute(data.precip.loc[t], dpt, precip_temp, ta, t,	1✔
498	data.wind_speed.loc[t], data.air_temp.loc[t],
499	az, dir_round_cell, flatwind, cell_maxus,
500	mask=mask)
501
502	queue[self.variable].put([t, self.precip])	1✔
503	queue['percent_snow'].put([t, self.percent_snow])	1✔
504	queue['snow_density'].put([t, self.snow_density])	1✔
505	queue['last_storm_day_basin'].put([t, self.last_storm_day_basin])	1✔
506	queue['storm_days'].put([t, self.storm_days])	1✔
507	queue['storm_total'].put([t, self.storm_total])	1✔
508	if self.nasde_model == "marks2017":	1✔
509	queue['storm_id'].put([t, self.storm_id])	1✔
510
511	def post_processor(self, main_obj, threaded=False):	1✔
512
513	pass	1✔

USDA-ARS-NWRC / smrf / 328

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous