USDA-ARS-NWRC / smrf / 324

pending completion

Build # 324

Build Type

Pull #170

travis-ci-com

Committed by

web-flow

Commit Message

Merge 610a59fd5 into 4e1ebc5df

Pull Request Pull Request #170: SMRF package

Run Details

234 of 234 new or added lines in 50 files covered. (100.0%)

3418 of 4338 relevant lines covered (78.79%)

0.79 hits per line

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

34.67

/smrf/envphys/precip.py

'''
Created on Apr 15, 2015

@author: scott
'''

import numpy as np

from smrf.utils import utils


def catchment_ratios(ws, gauge_type, snowing):
    """
    Point models for catchment ratios of the
    """

    if gauge_type == "us_nws_8_shielded":
        if snowing:
            CR = np.exp(4.61 - 0.04*ws**1.75)
        else:
            CR = 101.04 - 5.62*ws

    elif gauge_type == "us_nws_8_unshielded":
        if snowing:
            CR = np.exp(4.61 - 0.16*ws**1.28)
        else:
            CR = 100.77 - 8.34*ws
    else:
        raise ValueError(
            "Unknown catchement adjustment model ----> {0}".format(gauge_type))

    #    elif type == "Hellmann unshielded":
    #        if snowing:
    #             CR = CR = 100.00 + 1.13*Ws - 19.45*Ws
    #        else:
    #             CR = 96.63 + 0.41*Ws - 9.84*Ws + 5.95 * Tmean

    # elif type == "nipher":
    #     if snowing:
    #          CR= 100.00 - 0.44*Ws - 1.98*Ws
    #     else:
    #          CR = 97.29 - 3.18*Ws+ 0.58* Tmax - 0.67*Tmin
    #
    # elif type == "tretyakov":
    #     if snowing:
    #          CR = 103.11 - 8.67 * Ws + 0.30 * Tmax
    #     else:
    #          CR =  96.99 - 4.46 *Ws + 0.88 * Tmax + 0.22*Tmin
    # Avoid corrupting data
    if np.isnan(CR):
        CR = 100.0
    # CR is in percent.
    CR = CR/100.0
    return CR


def adjust_for_undercatch(p_vec, wind, temp, sta_type, metadata):
    """
    Adjusts the vector precip station data for undercatchment. Relationships
    should be added to :func:`~smrf.envphys.precip.catchment_ratio`.

    Args:
        p_vec - The station vector data in pandas series
        wind -  The vector wind data
        temp - The vector air_temp data
        sta_type - A dictionary of station names and the type of correction
            to apply
        station_metadata - station metadata TODO merge in the station_dict
            info to metadata

    Returns:
        adj_precip - Adjust precip accoding to the corrections applied.
    """
    adj_precip = p_vec.copy()
    for sta in p_vec.index:
        if sta in temp.keys():
            T = temp[sta]
            if sta in wind.keys():
                ws = wind[sta]
                if ws > 6.0:
                    ws = 6.0

                if T < -0.5:
                    snowing = True
                else:
                    snowing = False

                if sta in sta_type.keys():
                    gauge_type = sta_type[sta]
                else:
                    gauge_type = sta_type['station_undercatch_model_default']

                cr = catchment_ratios(ws, gauge_type, snowing)
                adj_precip[sta] = p_vec[sta]/cr

    return adj_precip


def dist_precip_wind(precip, precip_temp, az, dir_round_cell, wind_speed,
                     cell_maxus, tbreak, tbreak_direction, veg_type, veg_fact,
                     cfg, mask=None):
    """
    Redistribute the precip based on wind speed and direciton
    to account for drifting.

    Args:
        precip:             distributed precip
        precip_temp:        precip_temp array
        az:                 wind direction
        dir_round_cell:     from wind module
        wind_speed:         wind speed array
        cell_maxus:         max upwind slope from maxus file
        tbreak:             relative local slope from tbreak file
        tbreak_direction:   direction array from tbreak file
        veg_type:           user input veg types to correct
        veg_factor:         interception correction for veg types. ppt mult is
                            multiplied by 1/veg_factor

    Returns:
        precip_drift:       numpy array of precip redistributed for wind

    """
    # thresholds
    tbreak_threshold = cfg['tbreak_threshold']
    min_scour = cfg['winstral_min_scour']
    max_scour = cfg['winstral_max_scour']
    min_drift = cfg['winstral_min_drift']
    max_drift = cfg['winstral_max_drift']
    precip_temp_threshold = 0.5

    # polynomial factors
    drift_poly = {}
    drift_poly['a'] = 0.0289
    drift_poly['b'] = -0.0956
    drift_poly['c'] = 1.000761
    ppt_poly = {}
    ppt_poly['a'] = 0.0001737
    ppt_poly['b'] = 0.002549
    ppt_poly['c'] = 0.03265
    ppt_poly['d'] = 0.5929

    # initialize arrays
    celltbreak = np.ones(dir_round_cell.shape)
    drift_factor = np.ones(dir_round_cell.shape)
    precip_drift = np.zeros(dir_round_cell.shape)
    pptmult = np.ones(dir_round_cell.shape)

    # classify tbreak
    dir_unique = np.unique(dir_round_cell)

    for d in dir_unique:
        # find all values for matching direction
        ind = dir_round_cell == d
        i = np.argwhere(tbreak_direction == d)[0][0]
        celltbreak[ind] = tbreak[i][ind]

    # ################################################### #
    # Routine for drift cells
    # ################################################### #
    # for tbreak cells >  threshold
    idx = ((celltbreak > tbreak_threshold) & (
        precip_temp < precip_temp_threshold))

    # calculate drift factor
    drift_factor[idx] = drift_poly['c'] * np.exp(
        drift_poly['b'] * wind_speed[idx] +
        drift_poly['a'] * wind_speed[idx]**2
    )

    # cap drift factor
    drift_factor[idx] = utils.set_min_max(
        drift_factor[idx], min_drift, max_drift)

    # multiply precip by drift factor for drift cells
    precip_drift[idx] = drift_factor[idx]*precip[idx]

    # ################################################### #
    # Now lets deal with non drift cells
    # ################################################### #

    # for tbreak cells <= threshold (i.e. the rest of them)
    idx = ((celltbreak <= tbreak_threshold) & (
        precip_temp < precip_temp_threshold))

    # reset pptmult for exposed pixels
    pptmult = np.ones(dir_round_cell.shape)

    # original from manuscript
    pptmult[idx] = ppt_poly['a'] + ppt_poly['c'] * cell_maxus[idx] - \
        ppt_poly['b'] * cell_maxus[idx]**2 + ppt_poly['a'] * cell_maxus[idx]**3

    # veg effects at indices that we are working on where veg type matches
    for i, v in enumerate(veg_fact):
        idv = ((veg_type == int(v)) & (idx))
        pptmult[idv] = pptmult[idv] * 1.0 / veg_fact[v]

    # hardcode for pine
    pptmult[(idx) & ((veg_type == 3055) | (veg_type == 3053)
                     | (veg_type == 42))] = 1.0  # 0.92
    # cap ppt mult
    pptmult[idx] = utils.set_min_max(pptmult[idx], min_scour, max_scour)

    # multiply precip by scour factor
    precip_drift[idx] = pptmult[idx] * precip[idx]

    # ############################## #
    # no precip redistribution where dew point >= threshold
    idx = precip_temp >= precip_temp_threshold
    precip_drift[idx] = precip[idx]

    return precip_drift

1	'''
2	Created on Apr 15, 2015
3
4	@author: scott
5	'''
6
7	import numpy as np	1✔
8
9	from smrf.utils import utils	1✔
10
11
12	def catchment_ratios(ws, gauge_type, snowing):	1✔
13	"""
14	Point models for catchment ratios of the
15	"""
16
17	if gauge_type == "us_nws_8_shielded":	1✔
18	if snowing:	1✔
19	CR = np.exp(4.61 - 0.04ws*1.75)	1✔
20	else:
21	CR = 101.04 - 5.62*ws	×
22
23	elif gauge_type == "us_nws_8_unshielded":	×
24	if snowing:	×
25	CR = np.exp(4.61 - 0.16ws*1.28)	×
26	else:
27	CR = 100.77 - 8.34*ws	×
28	else:
29	raise ValueError(	×
30	"Unknown catchement adjustment model ----> {0}".format(gauge_type))
31
32	# elif type == "Hellmann unshielded":
33	# if snowing:
34	# CR = CR = 100.00 + 1.13Ws - 19.45Ws
35	# else:
36	# CR = 96.63 + 0.41Ws - 9.84Ws + 5.95 * Tmean
37
38	# elif type == "nipher":
39	# if snowing:
40	# CR= 100.00 - 0.44Ws - 1.98Ws
41	# else:
42	# CR = 97.29 - 3.18Ws+ 0.58 Tmax - 0.67*Tmin
43	#
44	# elif type == "tretyakov":
45	# if snowing:
46	# CR = 103.11 - 8.67 * Ws + 0.30 * Tmax
47	# else:
48	# CR = 96.99 - 4.46 Ws + 0.88 Tmax + 0.22*Tmin
49	# Avoid corrupting data
50	if np.isnan(CR):	1✔
51	CR = 100.0	×
52	# CR is in percent.
53	CR = CR/100.0	1✔
54	return CR	1✔
55
56
57	def adjust_for_undercatch(p_vec, wind, temp, sta_type, metadata):	1✔
58	"""
59	Adjusts the vector precip station data for undercatchment. Relationships
60	should be added to :func:`~smrf.envphys.precip.catchment_ratio`.
61
62	Args:
63	p_vec - The station vector data in pandas series
64	wind - The vector wind data
65	temp - The vector air_temp data
66	sta_type - A dictionary of station names and the type of correction
67	to apply
68	station_metadata - station metadata TODO merge in the station_dict
69	info to metadata
70
71	Returns:
72	adj_precip - Adjust precip accoding to the corrections applied.
73	"""
74	adj_precip = p_vec.copy()	1✔
75	for sta in p_vec.index:	1✔
76	if sta in temp.keys():	1✔
77	T = temp[sta]	1✔
78	if sta in wind.keys():	1✔
79	ws = wind[sta]	1✔
80	if ws > 6.0:	1✔
81	ws = 6.0	1✔
82
83	if T < -0.5:	1✔
84	snowing = True	1✔
85	else:
86	snowing = False	×
87
88	if sta in sta_type.keys():	1✔
89	gauge_type = sta_type[sta]	×
90	else:
91	gauge_type = sta_type['station_undercatch_model_default']	1✔
92
93	cr = catchment_ratios(ws, gauge_type, snowing)	1✔
94	adj_precip[sta] = p_vec[sta]/cr	1✔
95
96	return adj_precip	1✔
97
98
99	def dist_precip_wind(precip, precip_temp, az, dir_round_cell, wind_speed,	1✔
100	cell_maxus, tbreak, tbreak_direction, veg_type, veg_fact,
101	cfg, mask=None):
102	"""
103	Redistribute the precip based on wind speed and direciton
104	to account for drifting.
105
106	Args:
107	precip: distributed precip
108	precip_temp: precip_temp array
109	az: wind direction
110	dir_round_cell: from wind module
111	wind_speed: wind speed array
112	cell_maxus: max upwind slope from maxus file
113	tbreak: relative local slope from tbreak file
114	tbreak_direction: direction array from tbreak file
115	veg_type: user input veg types to correct
116	veg_factor: interception correction for veg types. ppt mult is
117	multiplied by 1/veg_factor
118
119	Returns:
120	precip_drift: numpy array of precip redistributed for wind
121
122	"""
123	# thresholds
124	tbreak_threshold = cfg['tbreak_threshold']	×
125	min_scour = cfg['winstral_min_scour']	×
126	max_scour = cfg['winstral_max_scour']	×
127	min_drift = cfg['winstral_min_drift']	×
128	max_drift = cfg['winstral_max_drift']	×
129	precip_temp_threshold = 0.5	×
130
131	# polynomial factors
132	drift_poly = {}	×
133	drift_poly['a'] = 0.0289	×
134	drift_poly['b'] = -0.0956	×
135	drift_poly['c'] = 1.000761	×
136	ppt_poly = {}	×
137	ppt_poly['a'] = 0.0001737	×
138	ppt_poly['b'] = 0.002549	×
139	ppt_poly['c'] = 0.03265	×
140	ppt_poly['d'] = 0.5929	×
141
142	# initialize arrays
143	celltbreak = np.ones(dir_round_cell.shape)	×
144	drift_factor = np.ones(dir_round_cell.shape)	×
145	precip_drift = np.zeros(dir_round_cell.shape)	×
146	pptmult = np.ones(dir_round_cell.shape)	×
147
148	# classify tbreak
149	dir_unique = np.unique(dir_round_cell)	×
150
151	for d in dir_unique:	×
152	# find all values for matching direction
153	ind = dir_round_cell == d	×
154	i = np.argwhere(tbreak_direction == d)[0][0]	×
155	celltbreak[ind] = tbreak[i][ind]	×
156
157	# ################################################### #
158	# Routine for drift cells
159	# ################################################### #
160	# for tbreak cells > threshold
161	idx = ((celltbreak > tbreak_threshold) & (	×
162	precip_temp < precip_temp_threshold))
163
164	# calculate drift factor
165	drift_factor[idx] = drift_poly['c'] * np.exp(	×
166	drift_poly['b'] * wind_speed[idx] +
167	drift_poly['a'] * wind_speed[idx]**2
168	)
169
170	# cap drift factor
171	drift_factor[idx] = utils.set_min_max(	×
172	drift_factor[idx], min_drift, max_drift)
173
174	# multiply precip by drift factor for drift cells
175	precip_drift[idx] = drift_factor[idx]*precip[idx]	×
176
177	# ################################################### #
178	# Now lets deal with non drift cells
179	# ################################################### #
180
181	# for tbreak cells <= threshold (i.e. the rest of them)
182	idx = ((celltbreak <= tbreak_threshold) & (	×
183	precip_temp < precip_temp_threshold))
184
185	# reset pptmult for exposed pixels
186	pptmult = np.ones(dir_round_cell.shape)	×
187
188	# original from manuscript
189	pptmult[idx] = ppt_poly['a'] + ppt_poly['c'] * cell_maxus[idx] - \	×
190	ppt_poly['b'] * cell_maxus[idx]*2 + ppt_poly['a'] cell_maxus[idx]**3
191
192	# veg effects at indices that we are working on where veg type matches
193	for i, v in enumerate(veg_fact):	×
194	idv = ((veg_type == int(v)) & (idx))	×
195	pptmult[idv] = pptmult[idv] * 1.0 / veg_fact[v]	×
196
197	# hardcode for pine
198	pptmult[(idx) & ((veg_type == 3055) \| (veg_type == 3053)	×
199	\| (veg_type == 42))] = 1.0 # 0.92
200	# cap ppt mult
201	pptmult[idx] = utils.set_min_max(pptmult[idx], min_scour, max_scour)	×
202
203	# multiply precip by scour factor
204	precip_drift[idx] = pptmult[idx] * precip[idx]	×
205
206	# ############################## #
207	# no precip redistribution where dew point >= threshold
208	idx = precip_temp >= precip_temp_threshold	×
209	precip_drift[idx] = precip[idx]	×
210
211	return precip_drift	×

USDA-ARS-NWRC / smrf / 324

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