diff --git a/core/logic/__init__.py b/core/logic/__init__.py index 2143d2e..38280bc 100644 --- a/core/logic/__init__.py +++ b/core/logic/__init__.py @@ -250,15 +250,12 @@ def main( # ----------------------------------------------------- # CALCULATE PEAK FLOW FOR ALL PRECIP PERIODS - - - # with everything generate peak flow estimates for the catchment + # generate peak flow estimates for the catchment for all storm frequencies in QP_HEADER peak_flow_ests = calculate_peak_flow( catchment_area_sqkm=each_catchment['area_up'], tc_hr=time_of_concentration, avg_cn=each_catchment['avg_cn'], precip_table=precip_tab_1d, - uid=each_catchment['id'], qp_header=QP_HEADER ) diff --git a/core/logic/calc.py b/core/logic/calc.py index 11e3779..c0afbd4 100644 --- a/core/logic/calc.py +++ b/core/logic/calc.py @@ -13,22 +13,23 @@ from .utils import msg -QP_HEADER=['Y1','Y2','Y5','Y10','Y25','Y50','Y100','Y200','Y500','Y1000'] +QP_HEADER=['Y1','Y2','Y5','Y10','Y25','Y50','Y100','Y200','Y500','Y1000'] def calculate_tc( max_flow_length, #units of meters mean_slope, # percent slope - const_a=0.000325, - const_b=0.77, + const_a=0.000325, + const_b=0.77, const_c=-0.385 -): + ): """ calculate time of concentration (hourly) Inputs: - max_flow_length: maximum flow length of a catchment area, derived from the DEM for the catchment area. - - mean_slope: average slope, from the DEM *for just the catchment area* + - mean_slope: average slope, from the DEM *for just the catchment area*. This must be + percent slope, provided as an integer (e.g., 23, not 0.23) Outputs: tc_hr: time of concentration (hourly) @@ -39,30 +40,33 @@ def calculate_tc( return tc_hr def calculate_peak_flow( - catchment_area_sqkm, - tc_hr, - avg_cn, + catchment_area_sqkm, + tc_hr, + avg_cn, precip_table, - uid=None, qp_header=QP_HEADER ): - """ - calculate peak runoff statistics at a "pour point" (e.g., a stormwater + """Calculate peak runoff statistics at a "pour point" (e.g., a stormwater inlet, a culvert, or otherwise a basin's outlet of some sort) using parameters dervied from prior analysis of that pour point's catchment area - (i.e., it's watershed or contributing area) and precipitation estimates. + (i.e., it's watershed or contributing area) and *24-hour* precipitation estimates. + + Note that the TR-55 methodology is designed around a 24-hour storm *duration*. YMMV + if providing rainfall estimates (via the precip_table parameter) for other storm durations. + + This calculator by default returns peak flow for storm *frequencies* ranging from 1 to 1000 year events. Inputs: - catchment_area_sqkm: area measurement of catchment in *square kilometers* - tc_hr: hourly time of concentration number for the catchment area - avg_cn: average curve number for the catchment area - - precip_table: a precipitation frequency estimates "table" as a Numpy - Array, derived from standard NOAA Preciptation Frequency Estimates + - precip_table: precipitation estimates as a 1D array (a list) + derived from standard NOAA Preciptation Frequency Estimates. Values in centimeters tables (the `precip_table_etl()` function can automatically prep this) Outputs: - runoff: a dictionary indicating peak runoff at the pour point for - various storm events by duration and frequency + storm events by frequency """ # reference some variables: @@ -75,53 +79,55 @@ def calculate_peak_flow( # (this indicates invalid data) if cn in [0,'',None] or tc in [0,'',None]: qp_data = [0 for i in range(0,len(qp_header))] - return OrderedDict(zip(qp_header,qp_data)) + return OrderedDict(zip(qp_header, qp_data)) # array for storing peak flows Qp = [] # calculate storage, S in cm # NOTE: THIS ASSUMES THE CURVE NUMBER RASTER IS IN METERS - Storage = 0.1*((25400.0/cn)-254.0) #cn is the average curve number of the catchment area - #msg("\tStorage: {0}".format(Storage)) - Ia = 0.2*Storage #inital abstraction, amount of precip that never has a chance to become runoff - #msg("\tIa: {0}".format(Ia)) - #setup precip list for the correct watershed from dictionary - P = numpy.array(precip_table) - #msg("\tP: {0}".format(P)) + Storage = 0.1 * ((25400.0 / cn) - 254.0) #cn is the average curve number of the catchment area + #msg("Storage: {0}".format(Storage)) + Ia = 0.2 * Storage #inital abstraction, amount of precip that never has a chance to become runoff + #msg("Ia: {0}".format(Ia)) + # setup precip list for the correct watershed from dictionary + P = numpy.array(precip_table) #P in cm + #msg("P: {0}".format(P)) + #calculate depth of runoff from each storm #if P < Ia NO runoff is produced - #P in cm Pe = (P - Ia) Pe = numpy.array([0 if i < 0 else i for i in Pe]) # get rid of negative Pe's - #msg("\tPe: {0}".format(Pe)) + #msg("Pe: {0}".format(Pe)) Q = (Pe**2)/(P+(Storage-Ia)) - #msg("\tQ: {0}".format(Q)) + #msg("Q: {0}".format(Q)) - #calculate q_peak, cubic meters per second + # calculate q_peak, cubic meters per second # q_u is an adjustment because these watersheds are very small. It is a function of tc, - # and constants Const0, Const1, and Const2 which are in turn functions of Ia/P (rain_ratio) and rainfall type - # We are using rainfall Type II because that is applicable to most of New York State - # rain_ratio is a vector with one element per input return period + # and constants Const0, Const1, and Const2 which are in turn functions of Ia/P (rain_ratio) and rainfall type + # We are using rainfall Type II because that is applicable to most of New York State + # rain_ratio is a vector with one element per input return period rain_ratio = Ia/P rain_ratio = numpy.array([.1 if i < .1 else .5 if i > .5 else i for i in rain_ratio]) # keep rain ratio within limits set by TR55 - #msg("\tRain Ratio: {0}".format(rain_ratio)) + #msg("Rain Ratio: {0}".format(rain_ratio)) - Const0 = (rain_ratio ** 2) * -2.2349 + (rain_ratio * 0.4759) + 2.5273 - Const1 = (rain_ratio ** 2) * 1.5555 - (rain_ratio * 0.7081) - 0.5584 - Const2 = (rain_ratio ** 2)* 0.6041 + (rain_ratio * 0.0437) - 0.1761 + # TODO: expose these as parameters; document here what they are (referencing the TR-55 documentation) + # TODO: some of these are geographically-derived; use geodata to pull the correct/suggested ones in (possibly + # in a function that precedes this) + Const0 = (rain_ratio**2) * -2.2349 + (rain_ratio * 0.4759) + 2.5273 + Const1 = (rain_ratio**2) * 1.5555 - (rain_ratio * 0.7081) - 0.5584 + Const2 = (rain_ratio**2) * 0.6041 + (rain_ratio * 0.0437) - 0.1761 #qu has weird units which take care of the difference between Q in cm and area in km2 (m^3 s^-1 km^-2 cm^-1) - qu = 10 ** (Const0+Const1*numpy.log10(tc)+Const2*(numpy.log10(tc))**2-2.366) - #msg("\tqu: {0}".format(qu)) - q_peak = Q*qu*catchment_area_sqkm - Qp = q_peak.tolist() # m^3 s^-1 - #msg("\tq_peak: {0}".format(Qp)) + qu = 10 ** (Const0 + Const1 * numpy.log10(tc) + Const2 * (numpy.log10(tc))**2 - 2.366) + #msg("qu: {0}".format(qu)) + q_peak = Q * qu * catchment_area_sqkm # m^3 s^-1 + #msg("q_peak: {0}".format(q_peak.tolist())) # TODO: better parameterize the header here (goes all the way back to how NOAA csv is ingested) results = OrderedDict(zip(qp_header,q_peak)) #msg("Results:") # for i in results.items(): - #msg("\t%-5s: %s" % (i[0], i[1])) + #msg("%-5s: %s" % (i[0], i[1])) return results