<jupyter_start><jupyter_text>ESCI 895 Final Project<NAME>  Project Summary:    The goal of this project is to determine the discharge associated with the 2, 5, 10, 25, 50, and 100-year floods on the Lamprey and Sugar Rivers in New Hampshire. Two methods will be compared, one using annual maximum discharge records and the other using instaneous discharge data. The first method will fit the data to a Gumbel distribution. The second with use the General Pareto distribtuion (GPD) with a Peaks-Over-Threshold (POT) method.        The annual maxiumum data has a longer record period for both rivers but only takes into account one large flood per year. While unlikely, extreme floods that happen within the same year will not be accounted for in the analysis. The instaneous data, using a peaks over threshold method to select only large floods, does not have this issue. The record length for this data set is much shorter, which may not provide accurate estimates for return periods of floods longer than the record itself, like the 25-year, 50-year, and 100-year floods.         The instananoeus data also requires more processing to select the maximum values from individual floods events, rather than analyzing multiple values measured duing the same flood event that may be above the threshold. A threshold value will also need to be determined. Other studies have recommend using a threshold that is equal to the magnitude of the discharge asscoiated 1-year to 2-year flood.     Study Sites: UPDATE TO INCLUDE BOTH INSTANEOUS AND ANNUAL DATA SETS AND RECORD LENGTHS Lamprey River:	The Lamprey River is located in southeastern New Hampshire. It 47 miles long, flowing from Northwood, NH to Newmarket, NH where it empties into the Great Bay. The Lamprey River watershed is 212 square miles. Daily discharge data was recorded at the USGS 01060003 station located in Newmarket, NH. This station has 86 years of consecutive daily discharge data (USGS). The gauge station is located directly upstream of Packer’s Falls. The drainage area upstream of the flow gauge is 185 square miles. The Lamprey River has some regulation by Pawtuckaway and Mendums Ponds. Occasionally, water from the river is diverted for the municipal supply of Durham. Newmarket, NH has a humid continental climate with hot summers and cold winters. Newmarket, NH has an average annual rainfall of 48 inches per year and snowfall of 53 inches per year.Add a figure Sugar River:	The Sugar River is located in western New Hampshire. It originates as an outlet of Lake Sunapee in Sunapee, NH, and heads west through Newport and Claremont, NH before emptying into the Connecticut River. The River is 27 miles long. The USGS gauge used in this project is located in West Claremont, NH upstream of a small waterfall. The gauge station is a few miles downstream of a large waterfall located in Claremont, NH. The gauge has a station identification number of USGS 01152500. The drainage basin draining to the location of this gauge is 269 square miles (USGS). The Sugar River is regulated by Sunapee Lake and occasionally diverted by mills upstream. This station has recorded 93 years of daily discharge data. Claremont, NH has a humid continental climate with a yearly average of 41.3 in of rain and 74.7 in of snowfall.Add a figure  Methods:     Data was analyzed using python ____ (say edition used?). ____ were imported to analyze the data in following steps.<jupyter_code>#%% Import libraries
import numpy as np
from matplotlib import pyplot as plt
import math
from scipy.special import gamma, factorial
import pandas as pd
import datetime as datetime<jupyter_output><empty_output><jupyter_text>Annual maximum discharge and instaneous discharge data sets were retrieved from USGS databases for each gage and loaded into python jupyter. Excess data was removed, leaving only the timestamps and discharge values. Missing or corrupt data was replaced using a linear interpolation method.<jupyter_code>#%% Specify inputs
filenames = ['hourlyqlamprey.txt', 'SugarInstant.txt', 'annualn.txt']
rivers = ["Lamprey River", "Sugar River"]<jupyter_output><empty_output><jupyter_text>First, data from the annual maximum discharge values were analyzed.     Second, instaneous discharge data was analyzed using the POT method and General Pareto Distribution (GPD). First, A threshold value needs to be selected for the POT method. Then values below the threshold can be removed for the rest of this analysis. ***Data is both in 1 hour and 15 min increments so need to remove data not on hour increments<jupyter_code>#%% Part 1
def GEV(file, threshold):
    #Load discharge file into dfpeak and format
    dfpeak = pd.read_csv(filenames[0], delimiter="\t", comment='#', header=1, parse_dates=['20d'], na_values = (9999, 999, 997, "Eqp"))
    #Rename columns
    dfpeak = dfpeak.rename(columns={"20d": "DATE","14n": "discharge_cfs"})
    #Set date column as index
    dfpeak =dfpeak.set_index('DATE')
    #Remove not needed columns
    dfpeak = dfpeak[['discharge_cfs']]
    #Fill nan values through linear interpolation
    dfpeak.interpolate(method = 'linear',inplace = True)
    
    
    #POT 
    #select discharges greater than threshold
    #threshold = 250
    #Drop all values below the threshold which are not needed for the extreme events analysis.
    dfpeak.drop(dfpeak[dfpeak['discharge_cfs'] < threshold].index, inplace = True)
    #Resample data to daily max daily data.
    #This accounts for uneven sampling intervals for data sets with hourly, 30 minute, and 15 minute measurement intervals.
    #This also helps simply analysis by lessening the amount of data values. This can be done as multiple extreme flood events wouldn't likely happen in the same day).
    dfpeak = dfpeak.resample('D').max()
    #Drop values below threshold.
    dfpeak = dfpeak.dropna()
    
    #Still need to select ind. floods events
    #something like if the row before is nan, and row below
    #is not then select max of all values until next nan?
    #But what if there are two major floods and discharge
    #never drops below threhold between them?
    
    #Running frequency analysis
    #Make a datframe to fill in later
    interp = np.array([2,5,10,25,50,100,200,500,1000])
    dfinterp = pd.DataFrame(interp, columns=['Return Period (yrs)'])
    dfinterp['EP'] = 1/dfinterp['Return Period (yrs)']
    dfinterp['1 - EP'] = 1 - dfinterp['EP']
    #Take mean and std of peak annual discharge data
    peak_mean = dfpeak['discharge_cfs'].mean()
    peak_std = dfpeak['discharge_cfs'].std()
    
    #Rank and order data
    count = dfpeak['discharge_cfs'].count()
    dfpeak = dfpeak.sort_values('discharge_cfs', ascending = True)
    dfpeak['rank'] = dfpeak['discharge_cfs'].rank(ascending= False)
    dfpeak = dfpeak.sort_values(by = 'rank', ascending = True)
    #Calculate L moments
    num = len(dfpeak)
    dfpeak['b1'] = ((num-dfpeak['rank'])/(num*(num-1)))*dfpeak['discharge_cfs']
    dfpeak['b2'] = (((num-dfpeak['rank'])*(num-dfpeak['rank']-1))/(num*(num-1)*(num-2)))*dfpeak['discharge_cfs']
    B1 = sum(dfpeak['b1'])
    B2 = sum(dfpeak['b2'])
    
    lamda1 = np.mean(dfpeak['discharge_cfs'])
    lamda2 = 2*B1-lamda1
    lamda3 = 6*B2-6*B1+lamda1
    skew = lamda3/lamda2
    
    #GEV - need to change to be GPD
    c = 2*lamda2/(lamda3+3*lamda2)-np.log(2)/np.log(3)
    k = 7.859*c+2.9554*c**2
    alpha = (k*lamda2)/(gamma(1+k)*(1-2**(-k)))
    squiggle = lamda1 + (alpha/k)*(math.gamma(1+k)-1)
    print(c , k , alpha, squiggle)
    dfinterp['GEV'] = squiggle+(alpha/k)*(1-(-np.log(dfinterp['1 - EP']))**k)
    #Create figures to show this
    print(dfinterp)<jupyter_output><empty_output><jupyter_text>GEV analysis for the Lamprey River:Selecting a threshold for the Lamprey River... Selected 250 cfs, which is the 1.5-year flood estimate from the previous analysis using annual maximum data and the GEV??? method.<jupyter_code>GEV(filenames[0], 250)<jupyter_output>-0.05500071989807398 -0.4233103386428873 179.62098848316572 410.9113332392692
   Return Period (yrs)     EP  1 - EP          GEV
0                    2  0.500   0.500   482.126415
1                    5  0.200   0.800   787.248687
2                   10  0.100   0.900  1086.627455
3                   25  0.040   0.960  1629.902589
4                   50  0.020   0.980  2199.864163
5                  100  0.010   0.990  2960.962434
6                  200  0.005   0.995  3979.477505
7                  500  0.002   0.998  5875.241837
8                 1000  0.001   0.999  7884.952900<jupyter_text>GEV analysis for the Sugar River:Selecting a threshold for the Sugar River... The selected threshold values was 400?? cfs, which is the 1.5-year flood estimate from the previous analysis using annual maximum data and the GEV??? method.<jupyter_code>GEV(filenames[1], 400)
#Show a figure of the distribution... 
#maybe highlight/show all the extreme values in different color<jupyter_output><empty_output><jupyter_text>***Ask Anne about how to talk about steps done within the function (like it as a function because its easy to run for different data sets but isn't conducive to juptyer format for explanation then code blocks.  Results:<jupyter_code>#display dfinterp data frames... make one that shows values for both distributions.

#For GPD can show CDF and PDF graphs<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07615787784789817, 'fraction_numerical': 0.047598673654936356, 'mean_word_length': 3.839544513457557, 'pattern_counts': {'":': 2, '<': 26, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '403602', 'license': 'MIT', 'n_tokens_mistral': 3057, 'n_tokens_neox': 2636, 'n_words': 1191, 'path': 'ESCI_895_Project_LK.ipynb', 'repo_name': 'lak1039/ESCI895', 'size': '9351'}