# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # GE 254 - Intro to Geomorph
# ## HW 1 - Using Python to calculate Flood Frequency Analysis
#
# #### Due: Thursday Sept. 24th by 11 am
#
# For this assignment, we will be taking the work we did in excel on flood frequency and instead creating a python script to do the same thing. You will be given a partially done code. You need to finish the code. Answer any questions shown as comments in the code. And most importantly, comment the code extensively. The end result should be a useful python script that takes any USGS flow data and calculates predictive flood frequency analysis. 
#
# This exercise will provide some experience with methods used for predicting flood frequency and magnitude. We will be using the US Geological Survey (USGS) website to retrieve historical stream gauge data of the sort used to predict the likelihood of flood events of particular magnitudes during a given time interval. Such predictions are the basis for numerous engineering, restoration and development projects in and around rivers.

# #### Objectives:
#     1. Practice using Python Pandas DataFrames to import, manipulate, and visualize data
#     2. Learn some powerful tools for subsetting your dataframes
#     3. Practice visualizing data 
#     
# * Dr <NAME>, Sept. 2020
# * Homework assignment adapted from T. Perron Lab on flood frequency
# * Used in a 200 level undergraduate geomorphology course, after doing the same work in lab using excel (or google sheets) to calculate flood frequency and rating curves. 
# * This homework assumes a basic introduction to python, pandas, and dataframes.

# +
#import libraries

#let's use pandas dataframes again
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt 

#so I've imported 3 libraries using common 
#abbreviations to reference these libraries in my code. This is personal choice.
#If I just wrote "import numpy", it would still work but when I call a numpy command
#I would have to write numpy.matrix() vs. np.matrix() - it's all about laziness
#and minimizing typing.
# -

# First make sure you've downloaded the correct data. I want you to download peak flow and daily flow measurements for 1 of these four sites (you will be assigned in class): 
# * 02086849 
# * 02081000
# * 02081747  
# * 02083000 
#
# I also want you to download daily data for your site, on the waterdata.usgs.gov/nwis/sw site click on daily data instead of peak-flow data. 
#
# Remember to download all available data for discharge and gage height (if available). Remember to select tab-separated file then download that file as a .txt, by right clicking on the page and selecting "save as".  Make sure you upload these files into JupyterHub. 

# Let's start by reading in the peak-flow measurements - this is similar to what we did in Lab 2

# +
#ok let's break down this command - can you tell me what header does?
#what is the delimiter mean?
# why did I write skiprows?
#what is usecols do?
#for help with this use google - or go to the help menu above and select pandas

#don't forget to add your file name to ADD_FILE_NAME_HERE.txt part of the '' in the line below
peak_all = pd.read_csv('ADD_FILE_NAME_HERE.txt',header=72,delimiter="\t",\
                       skiprows=[73],usecols=range(0,8),parse_dates=True)
peak_all.head()

#if for some reason when peak_all is displayed below and you see weird values
#in the first couple of rows, you need to adjust the values given to header
#and skiprows. 
# -

#now let's rename these columns to something a bit more useful
new_column_names = ['Agency', 'SiteNo', 'Date', 'Time', \
                    'Discharge (cfs)', 'Discharge_quality','Gage_ht (ft)',\
                    'Gage_quality']
#ok now what does the "\" do in the above line of code?

peak_all.columns = new_column_names

peak_all.head()

# +
#now for some data cleanup
peak_all['Discharge (m3/sec)'] = peak_all['Discharge (cfs)'] * 0.028316847
peak_all['Gage_ht (m)'] = peak_all['Gage_ht (ft)'] * 3.28084
new_station_name = "0" + str(peak_all['SiteNo'].unique()[0])
peak_all['SiteNo'] = new_station_name

#ok what happens in this cell?
peak_all.head()
# -

#the fun date-time work
peak_all['Date'] = pd.to_datetime(peak_all.Date)
peak_all['Year'] = peak_all['Date'].dt.year
peak_all.head()
#explain what the above lines of code do please

# +
#ok now that we've done some data management, let's pull out only the data we need
peak = peak_all[['Year','Discharge (m3/sec)','Gage_ht (m)']]

#now let's remove NaN  measurements
peak = peak.dropna()
print(peak.head())
# -

peak = peak.sort_values('Discharge (m3/sec)',ascending=False) #what does this do?
peak.head()

peak['Rank'] = range(1,peak.shape[0]+1) #what is the range command?
n = peak.shape[0] #what value(s) does shape give us? why do I index at first position (0)?
print('The total number of observations:', n)
print(peak.head())

# #### helpful description of indexing dataframes in python & pandas
# *https://www.shanelynn.ie/select-pandas-dataframe-rows-and-columns-using-iloc-loc-and-ix/*

# +
#ok now create a new column in peak dataframe that is the Recurrence Interval (RI = (1+n)/rank)


# -

# The usual procedure is to assume that the frequency distribution of floods in our record conforms to a known distribution, and plot the observed floods in such a way that they will fall along a straight line if the data conform to the assumed distribution. Distributions in common use include the Gumbel, Weibull, and Pearson Type III distributions, all of which are designed to characterize extreme value phenomena. To keep our analysis simple, we will assume that flood frequency is lognormally distributed. This implies that there should be a semilogarithmic relationship between flood magnitude and RI.

#ok now to plot data
peak.plot(x='Recurrence Interval (years)',y='Discharge (m3/sec)',\
          title='Flood Frequency of Station ' + peak_all['SiteNo'][0], \
          kind='scatter',logx=True)

#calculate the trendline to our floods
x = peak['Recurrence Interval (years)']
y = peak['Discharge (m3/sec)']
f = np.polyfit(np.log10(x),y,1,w=np.sqrt(y))
print(f)
xf = [min(x),max(x)]
yf = f[0]*np.log10(xf) + f[1]
print(xf,yf)

# +
#now plot the data and the trendline

#hint to add the trendline - use plt.plot(x,y) where you need to figure out what x and y should be!

# -

#What is the 100 year flood discharge? 
dis100 = 
print('the 100 year discharge is:', dis100, 'm3/s')

# +
#now lets look at the rating curve - aka stage vs. discharge


#can you fit a line to this plot? For predicting discharge for a given stage?

#Hint use polyfit again

#now plot the linear trendline on the rating curve
# -

#what is the predicted stage for the 100 year flood?
#stg100 = 
print('the 100 year flood stage is:', stg100, 'm')

# 1. Ok now that we've recreated flood recurrence interval analysis, what would you have to change to run this code on a different USGS station? 
# 2. What parameters might be different? What would you expect to be the same? 
# 3. How easy would it be to change these things?
#
# ## visualizing and analyzing daily flow measurements. 

#alright nice job - now let's look at daily measurements. Make sure to fill in the correct file name in the ''
#make sure you've uploaded the .txt file to jupyter hub
daily_all = pd.read_csv('',header=32,delimiter="\t",\
                        skiprows=[33],usecols=[0,1,2,3,9],parse_dates=True)
daily_all.head()
#again - note the header and skip-row number, these might need to change
#for your file...

new_column_names = ['Agency', 'SiteNo', 'OldDateTime', \
                    'Discharge (cfs)', 'Gage_ht (ft)']
daily_all.columns = new_column_names
daily_all.head()

# +
#now for some data cleanup

#first change all units for gage height & discharge to mks


#then add in the corrected station number

# -

# ## Questions about the above code:
# 1. Why do I keep printing out/displaying daily_all or peak or peak_all? 
# 2. What is the use or reason for this in the code? 
# 3. What are other methods you can use to check your code validity?

# +
#now calculate the datetime objects


#now add in separate columns for day, month, and year - year is shown below, do the other two
daily_all['Year'] = daily_all['DateTime'].dt.year



print(daily_all.head())

# +
#ok how can you now find the average flow for YOUR birthday? change the code as needed
avgdis_bday = daily_all['Discharge (m3/sec)'][((daily_all.Month==7)&(daily_all.Day==7))].mean()

print('The average discharge for 7/7 is: ', avgdis_bday, 'm3/sec')
#go ahead and redo this to find the minimum, maximum, and mean flow for YOUR birthday


# +
#now go ahead and plot the daily average flow (aka average per day)
#this is how we use the very useful function group by

avg_daily = daily_all.groupby(['Month','Day'],as_index=False).mean()
#ok how does the above line of code work? Look up groupby and explain what is done.

avg_daily.Year = 2000 #random year chosen - must be a leap year for datetime to work!

# -

avg_daily['Date'] = pd.to_datetime(avg_daily[['Year','Month','Day']])


# +
#now go ahead and plot your avg daily value
avg_daily.plot(x='Date',y='Discharge (m3/sec)',\
               title='Averaged Daily Flow of Station ' + peak_all['SiteNo'][0], \
              kind='scatter')

#now can you plot a monthly mean? Let's use the groupby again
avg_m = daily_all.groupby(['Month'],as_index=False).mean()
avg_m['Day'] = 15 #add in a fake day
avg_m['Year'] = 2000 #add in a fake year that fits the data - make sure it matches above
avg_m['Date'] = pd.to_datetime(avg_m[['Year','Month','Day']])
print(avg_m)
plt.plot(avg_m.Date,avg_m['Discharge (m3/sec)'],'--m',linewidth=4)

#does the monthly average match the daily averages values - are the trends holding?
# -

print(daily_all.head())

# ## Resampling Data
# Now let's calculate the per month average flow over our timeseries (aka resample our data)
# look at this or some helpful information *https://sergilehkyi.com/tips-on-working-with-datetime-index-in-pandas/* 

daily_all2 = daily_all.set_index('DateTime')
print(daily_all2.head())

ma = daily_all2.resample('M').mean() #look up what this command does - this can be very powerful

print(ma.head())
print(daily_all2.head()) 
#explain the two new dataframes created (ma and daily_all2). What information is in these?
#why did I create these? 

daily_all.plot(x='DateTime', y='Discharge (m3/sec)',kind='scatter')
plt.plot(ma.index,ma['Discharge (m3/sec)'],"--k",linewidth=3)

# ## Questions
# 1. What do you see in the above plots? 
# 2. What error might be found by the order in which we analyzed data?
# 3. Looking at this plot, do you think you might switch the order of operations of this homework? 
# 4. Would you maybe analyze or even visualize one piece of data before another? Why or why not?

# ## Ok last task:
#     1. Plot the rating curve based on the daily data (see daily_all)
#     2. On top of that, plot the rating curve based on the monthly averaged data (see ma dataframe)
#     3. Finally plot on top of that, the rating curve based on peak yearly flows (see peak)
#     4. What differences occur in these different datasets? Make sure you plot each dataset with a different color and/or marker to easily highlight the differences. Make sure to add a legend! 
#
# Ok now fit a linear trendlines to each dataset (you should have 3 in total, and remember, you've already calculated one of these trendlines earlier on). 
# * How do these trendlines differ? What might be driving these differences? 
# * What errors might be caused by using one dataset over the other? 
# * Which dataset would you use to most accurately predict discharge from the stage?

# Ok now send me a copy of your .ipynb file and the two .txt files you've analyzed in this code. Please remember to give it a unique name to the file name (aka add your initials). **Nice job with Homework 1!**



StarcoderdataJupyterScriptsDedupFiltered