import pandas as pd
from io import BytesIO
import pycurl

site_number = '01589440'
date_from = '1996-06-06'
date_to = '2006-06-06'

url = 'https://nwis.waterdata.usgs.gov/nwis/uv/?parameterCd=00060,00065&format=rdb&site_no=%s&period=&begin_date=%s\
&end_date=%s&siteStatus=all' % (site_number, date_from, date_to)

def download_data_from_url( url, savename='test.csv' ):
    ''' 
    One can download data with a different method that supports resume. If data is two large then it takes
    lot of time and the connection to the server might be interrupted.
    '''
    
    buffer = BytesIO()
    c = pycurl.Curl()
    c.setopt(c.URL, url)
    c.setopt(c.WRITEDATA, buffer)
    c.perform()
    c.close()

    body = buffer.getvalue() # Body is a byte string.
    with open( savename, 'w' ) as /assets/images/20190215/output:
        /assets/images/20190215/output.write( body.decode('utf-8')) # We have to know the encoding in order to print it to a text file
        
download_data_from_url(url, '%s.csv' % site_number)

data = pd.read_csv('%s.csv' % site_number, sep='\t', comment='#', skiprows=[29], parse_dates=[2])
data.head(n=10)

data['timestamp'] = data.datetime.map(lambda x: x.timestamp())
data.head(n=10)

data = data.rename(index=str, columns={"69744_00060": "Discharge"})
data.head(n=10)

from scipy.optimize import curve_fit

hist, bins = np.histogram(data['Discharge'], bins=500)

def quadratic(x, a, b, c):
    return a*x**2 + b*x + c

N = 30
# to fit the peak and a small portion of the tail after it
x = np.log(bins[:-1][:N])
y = np.log(hist[:N])

popt, pcov = curve_fit(quadratic, x, y)
perr = np.sqrt(np.diag(pcov))

print('a, b, c = ', popt)
print('err_a, err_b, err_c = ', perr)

x_test = np.linspace(0, 5, 20)
y_test = quadratic(x_test, *popt)

plt.plot(x_test, quadratic(x_test, *popt), 'r--')
plt.plot(x, y, 'o--')

plt.show()

data.plot('timestamp', 'Discharge')

from scipy.signal import find_peaks
import matplotlib.pyplot as plt
import numpy as np

x = data['Discharge']
ts = data['timestamp']
peaks, _ = find_peaks(x, height=60, width=15)

plt.figure(figsize=(15, 8))
plt.plot(ts, x)
plt.plot(ts[peaks], x[peaks], "x")
plt.xticks(np.linspace(0.85e9, 1.15e9, 5))
plt.show()

y_peak = x[peaks]

hist, bins = np.histogram(y_peak, bins=200)
plt.plot(bins[:-1], hist, 'x')

from scipy.optimize import curve_fit

N = 19
# to exclude zeros and ones
x = np.log(bins[:-1][:N])
y = np.log(hist[:N])

def lin(x, a, b):
    return a*x + b

popt, pcov = curve_fit(lin, x, y)
perr = np.sqrt(np.diag(pcov))
print('a, b = ', popt)
print('err_a, err_b = ', perr)

plt.plot(x, y, 'x')

x_test = np.linspace(2, 8, 20)
y_test = lin(x_test, *popt)

plt.plot(x_test, y_test)

event_distance = [end_date - start_date for start_date, end_date in zip(ts[peaks], ts[peaks][1:])]

# event distance in hours instead of seconds
hist, bins = np.histogram(np.array(event_distance) / 360., bins=200)

# to exclude zeros and ones
N = 20
x = np.log(bins[:-1][:N])
y = np.log(hist[:N])

popt, pcov = curve_fit(quadratic, x, y)
perr = np.sqrt(np.diag(pcov))

print('a, b, c = ', popt)
print('err_a, err_b, err_c = ', perr)

x_test = np.linspace(0, 9, 20)
y_test = quadratic(x_test, *popt)

plt.plot(x_test, quadratic(x_test, *popt), 'r--')

plt.plot(x, y, 'o')
plt.show()

### fifteen minute samples one year is approx. 24*4*365 = 35040

year = 35040
maximum_dc = []
years = data['Discharge'].shape[0] // year + 1

for yr in range(0, years):
    maximum_dc.append(np.max(data['Discharge'][yr*year:(yr+1)*year]))

plt.plot(maximum_dc, 'o--')
plt.xticks(np.linspace(0, 10, 11))

def sin(x, a, b, c, d):
    return a*np.sin(b*x + c) + d

popt, pcov = curve_fit(sin, np.linspace(0, 10, 10), maximum_dc, p0=[1000, 2.*np.pi/7., 5., 1200])
perr = np.sqrt(np.diag(pcov)) 0, 100)
y_test = sin(x_test, *popt)

plt.plot(x_test, y_test, 'r--')

plt.show()

from datetime import datetime
import collections

month = 2880
rainy_days = []
months = data['Discharge'].shape[0] // month  + 1
months_dict = dict()

for time_stamp in ts[peaks]:
    months_dict[datetime.utcfromtimestamp(int(time_stamp)).month] = 0

for time_stamp in ts[peaks]:
    months_dict[datetime.utcfromtimestamp(int(time_stamp)).month] += 1

D = collections.OrderedDict(sorted(months_dict.items()))
    
plt.bar(D.keys(), D.values(), align='center')
plt.show()

import pandas as pd

site_number = '04240010' # approx. 220 miles north - Syracuse

url = 'https://nwis.waterdata.usgs.gov/nwis/uv/?parameterCd=00060,00065&format=rdb&site_no=%s&period=&begin_date=%s\
&end_date=%s&siteStatus=all' % (site_number, date_from, date_to)

download_data_from_url(url, "%s.csv" % site_number)

data2 = pd.read_csv('%s.csv' % site_number, sep='\t', comment='#', skiprows=[32], parse_dates=[2])
data2.head(n=10)

data2['timestamp'] = data2.datetime.map(lambda x: x.timestamp())
data2 = data2.rename(index=str, columns={"107697_00060": "Discharge"})
data2.head(n=10)

# showing it on top of each other

plt.figure(figsize=(15, 8))
plt.plot(data["timestamp"], data["Discharge"] - 200, 'b--')
plt.plot(data2["timestamp"], data2["Discharge"], 'r--')
plt.show()

# weekly rolling window 24*4*7 = 672

### fifteen minute samples one year is approx. 24*4*365 = 35040

week = 672
maximum_dc = []
maximum_dc2 = []

# assuming same number of weeks
weeks = data['Discharge'].shape[0] // week + 1

for wk in range(0, weeks):
    maximum_dc.append(np.max(data['Discharge'][wk*week:(wk+1)*week]))
    maximum_dc2.append(np.max(data2['Discharge'][wk*week:(wk+1)*week]))

plt.figure(figsize=(15, 8))

from scipy.signal import savgol_filter
    
plt.plot(savgol_filter(maximum_dc, 31, 2), 'b--')
plt.plot(savgol_filter(maximum_dc2, 31, 2), 'r--')
plt.xticks(np.linspace(0, 500, 26))

plt.show()

shift = 20

plt.figure(figsize=(15, 8))
plt.plot(np.linspace(shift, 520, len(maximum_dc)), savgol_filter(maximum_dc, 31, 2), 'b--')
plt.plot(savgol_filter(maximum_dc2, 31, 2), 'r--')
plt.xticks(np.linspace(0, 500, 26))
plt.show()