""" Created on March 14, 2017 Originally written by Scott Havens in 2015 @author: Micah Johnson """ import os from datetime import datetime import numpy as np import pandas as pd import pytz [docs]def storms(precipitation, perc_snow, mass=1, time=4, stormDays=None, stormPrecip=None, ps_thresh=0.5): """ Calculate the decimal days since the last storm given a precip time series, percent snow, mass threshold, and time threshold - Will look for pixels where perc_snow > 50% as storm locations - A new storm will start if the mass at the pixel has exceeded the mass limit, this ensures that the enough has accumulated Args: precipitation: Precipitation values perc_snow: Precent of precipitation that was snow mass: Threshold for the mass to start a new storm time: Threshold for the time to start a new storm stormDays: If specified, this is the output from a previous run of storms stormPrecip: Keeps track of the total storm precip Returns: tuple: - **stormDays** - Array representing the days since the last storm at a pixel - **stormPrecip** - Array representing the precip accumulated during the most recent storm Created April 17, 2015 @author: Scott Havens """ # either preallocate or use the input if stormDays is None: stormDays = np.zeros(precipitation.shape) if stormPrecip is None: stormPrecip = np.zeros(precipitation.shape) # if there is no snow, don't reset the counter # This ensures that the albedo won't be reset stormDays += 1 if np.sum(perc_snow) == 0: # stormDays = np.add(stormDays, 1) stormPrecip = np.zeros(precipitation.shape) return stormDays, stormPrecip # determine locations where it has snowed idx = perc_snow >= ps_thresh # determine locations where the time threshold has passed # these areas, the stormPrecip will be set back to zero idx_time = stormDays >= time stormPrecip[idx_time] = 0 # add the values to the stormPrecip stormPrecip[idx] =+ precipitation[idx] # see if the mass threshold has been passed idx_mass = stormPrecip >= mass # reset the stormDays to zero where the storm is present stormDays[idx_mass] = 0 return stormDays, stormPrecip [docs]def time_since_storm(precipitation, perc_snow, time_step=1/24, mass=1.0, time=4, stormDays=None, stormPrecip=None, ps_thresh=0.5): """ Calculate the decimal days since the last storm given a precip time series, percent snow, mass threshold, and time threshold - Will look for pixels where perc_snow > 50% as storm locations - A new storm will start if the mass at the pixel has exceeded the mass limit, this ensures that the enough has accumulated Args: precipitation: Precipitation values perc_snow: Percent of precipitation that was snow time_step: Step in days of the model run mass: Threshold for the mass to start a new storm time: Threshold for the time to start a new storm stormDays: If specified, this is the output from a previous run of storms else it will be set to the date_time value stormPrecip: Keeps track of the total storm precip Returns: tuple: - **stormDays** - Array representing the days since the last storm at a pixel - **stormPrecip** - Array representing the precip accumulated during the most recent storm Created Janurary 5, 2016 @author: Scott Havens """ # either preallocate or use the input if stormDays is None: stormDays = np.zeros(precipitation.shape) if stormPrecip is None: stormPrecip = np.zeros(precipitation.shape) # if there is no snow, don't reset the counter # This ensures that the albedo won't be reset stormDays += time_step if np.sum(perc_snow) == 0: # stormDays = np.add(stormDays, 1) stormPrecip = np.zeros(precipitation.shape) return stormDays, stormPrecip # determine locations where it has snowed idx = perc_snow >= ps_thresh # determine locations where the time threshold has passed # these areas, the stormPrecip will be set back to zero idx_time = stormDays >= time stormPrecip[idx_time] = 0 # add the values to the stormPrecip stormPrecip[idx] =+ precipitation[idx] idx_mass = stormPrecip >= mass # reset the stormDays to zero where the storm is present stormDays[idx_mass] = 0 return stormDays, stormPrecip [docs]def time_since_storm_basin(precipitation, storm, stormid, storming, time, time_step=1/24, stormDays=None): """ Calculate the decimal days since the last storm given a precip time series, days since last storm in basin, and if it is currently storming - Will assign uniform decimal days since last storm to every pixel Args: precipitation: Precipitation values storm: current or most recent storm time_step: step in days of the model run last_storm_day_basin: time since last storm for the basin stormid: ID of current storm storming: if it is currently storming time: current time stormDays: unifrom days since last storm on pixel basis Returns: stormDays: unifrom days since last storm on pixel basis Created May 9, 2017 @author: Scott Havens modified by Micah Sandusky """ # either preallocate or use the input if stormDays is None: stormDays = np.zeros(precipitation.shape) #if before first storm, just add timestep if not storming and time <= storm['start'] and stormid == 0: stormDays += time_step # if during a storm than reset to zero elif time <= storm['end']: stormDays = np.zeros(precipitation.shape) # else assign uniform to days from last storm for the basin else: stormDays += time_step return stormDays [docs]def time_since_storm_pixel(precipitation, dpt, perc_snow, storming, time_step=1/24, stormDays=None, mass=1.0, ps_thresh=0.5): """ Calculate the decimal days since the last storm given a precip time series - Will assign decimal days since last storm to every pixel Args: precipitation: Precipitation values dpt: dew point values perc_snow: percent_snow values storming: if it is stomring time_step: step in days of the model run stormDays: unifrom days since last storm on pixel basis mass: Threshold for the mass to start a new storm ps_thresh: Threshold for percent_snow Returns: stormDays: days since last storm on pixel basis Created October 16, 2017 @author: Micah Sandusky """ # either preallocate or use the input if stormDays is None: stormDays = np.zeros(precipitation.shape) # add timestep stormDays += time_step # only reset if stomring and not overly warm if storming and dpt.min() < 2.0: # determine location where there is enough mass idx_mass = precipitation >= mass # determine locations where it has snowed idx = perc_snow >= ps_thresh # reset the stormDays to zero where the storm is present stormDays[(idx_mass & idx)] = 0 return stormDays [docs]def tracking_by_station(precip, mass_thresh = 0.01, steps_thresh = 3): """ Processes the vector station data prior to the data being distributed Args: precipitation: precipitation values time: Time step that smrf is on time_steps_since_precip: time steps since the last precipitation storm_lst: list that store the storm cycles in order. A storm is recorded by its start and its end. The list is passed by reference and modified internally. Each storm entry should be in the format of: [{start:Storm Start, end:Storm End}] e.g. [ {start:date_time1,end:date_time2,'BOG1':100, 'ATL1':85}, {start:date_time3,end:date_time4,'BOG1':50, 'ATL1':45}, ] #would be a two storms at stations BOG1 and ATL1 mass_thresh: mass amount that constitutes a real precip event, default = 0.01. steps_thresh: Number of time steps that constitutes the end of a precip event, default = 2 steps (typically 2 hours) Returns: tuple: - **storms** - A list of dictionaries containing storm start,stop, mass accumulated, of given storm. - **storm_count** - A total number of storms found Created April 24, 2017 @author: Micah Johnson """ storm_columns = ['start','end'] stations = list(precip) storm_columns+=stations storms = [] stations = list(precip) is_storming = False time_steps_since_precip= 0 tzinfo = pytz.timezone('UTC') for i,row in precip.iterrows(): time = pd.Timestamp(i) #Storm Idenificiation if row.max() > mass_thresh: #Start a new storm if not is_storming: new_storm = {} new_storm['start'] = time for sta,p in row.iteritems(): new_storm[sta] = 0 #Create a new row is_storming = True #print "=="*10 + "> New storm!" time_steps_since_precip = 0 #Always add the latest end date to avoid unclosed storms new_storm['end'] = time #Accumulate precip for storm total for sta,mass in row.iteritems(): new_storm[sta] += mass elif is_storming and time_steps_since_precip < steps_thresh: #storm_lst[-1]['end'] = time new_storm['end'] = time time_steps_since_precip+=1 #print "=="*10 +"> Hours since precip = {0}".format(time_steps_since_precip) #print "=="*10 + "> still storming but no precip!" if time_steps_since_precip >= steps_thresh and is_storming: is_storming = False storms.append(new_storm) #print "=="*10 + "> not storming!" #Append the last storm if we ended during a storm if is_storming: storms.append(new_storm) storm_count = len(storms) #Make sure we have storms if storm_count == 0: empty_data = {} for col in storm_columns: empty_data[col] = [] storms = pd.DataFrame(empty_data) else: storms = pd.DataFrame(storms) return storms,storm_count [docs]def tracking_by_basin(precipitation, time, storm_lst, time_steps_since_precip, is_storming, mass_thresh = 0.01, steps_thresh=2): """ Args: precipitation: precipitation values time: Time step that smrf is on time_steps_since_precip: time steps since the last precipitation storm_lst: list that store the storm cycles in order. A storm is recorded by its start and its end. The list is passed by reference and modified internally. Each storm entry should be in the format of: [{start:Storm Start, end:Storm End}] e.g. [ {start:date_time1,end:date_time2}, {start:date_time3,end:date_time4}, ] #would be a two storms mass_thresh: mass amount that constitutes a real precip event, default = 0.0. steps_thresh: Number of time steps that constitutes the end of a precip event, default = 2 steps (typically 2 hours) Returns: tuple: storm_lst - updated storm_lst time_steps_since_precip - updated time_steps_since_precip is_storming - True or False whether the storm is ongoing or not Created March 3, 2017 @author: Micah Johnson """ # print "--"*10 +"> Max precip = {0}".format(precipitation.max()) if precipitation.max() > mass_thresh: #Start a new storm if len(storm_lst)== 0 or not is_storming : storm_lst.append({'start':time,'end':None}) is_storming = True # print "--"*10 + "> New storm!" #always append the most recent timestep to avoid unended storms storm_lst[-1]['end'] = time time_steps_since_precip = 0 #print "--"*10 + "> still storming!" elif is_storming and time_steps_since_precip < steps_thresh: #storm_lst[-1]['end'] = time time_steps_since_precip+=1 # print "--"*10 +"> Hours since precip = {0}".format(time_steps_since_precip) # print "--"*10 + "> still storming but no precip!" if time_steps_since_precip >= steps_thresh: is_storming = False # print "--"*10 + "> not storming!" return storm_lst, time_steps_since_precip, is_storming [docs]def clip_and_correct(precip,storms,stations = []): """ Meant to go along with the storm tracking, we correct the data here by adding in the precip we would miss by ignoring it. This is mostly because will get rain on snow events when there is snow because of the storm definitions and still try to distribute precip data. Args: precip: Vector station data representing the measured precipitation storms: Storm list with dictionaries as defined in :func:`~smrf.envphys.storms.tracking_by_station` stations: Desired stations that are being used for clipping. If stations is not passed, then use all in the dataframe Returns: The correct precip that ensures there is no precip outside of the defined storms with the clipped amount of precip proportionally added back to storms. Created May 3, 2017 @author: Micah Johnson """ #Specify zeros where were not storming precip_clipped = precip.copy() precip_clipped[:] = 0 for j,storm in storms.iterrows(): storm_start = storm['start'] storm_end = storm['end'] my_slice = precip.loc[storm_start:storm_end] precip_clipped.loc[storm_start:storm_end] = my_slice correction = {} #print "{0:>10} {1:>10} {2:>10}".format("Original","Clipped", "C") if len(stations) == 0: stations = precip.columns #Correct the precip to be equal to the sum. for station in stations: original = precip[station].sum() clipped = precip_clipped[station].sum() if original == 0: c = 1.0 elif clipped == 0: c = 0 else: c = original/clipped #print "{0:>10} {1:>10} {2:>10}".format(original,clipped,c) correction[station] = c #print "Max precip from the correction:{0}".format(precip_clipped.max()) #return pd.Series(correction) return precip_clipped.mul(pd.Series(correction), axis=1)