import matplotlib.pyplot as plt
from salishsea_tools.nowcast import figures, residuals
import numpy as np
import datetime
from dateutil import tz

from salishsea_tools import stormtools
import arrow
import pandas as pd
import requests
from cStringIO import StringIO
import matplotlib

%matplotlib inline

start = '01-Jan-2014'; end = '31-Dec-2014'

wlev_2014 = figures.load_archived_observations('Tofino',start,end )
mean_2014 =  wlev_2014.wlev.mean()

start = '01-Jan-2005'; end = '31-Dec-2014'

wlev_10yr = figures.load_archived_observations('Tofino',start,end )

mean_10yr = wlev_10yr.wlev.mean()

fig,ax=plt.subplots(1,1,figsize=(10,5))


s=wlev_2014['time'][wlev_2014.index[0]]
e=wlev_2014['time'][wlev_2014.index[-1]]

ax.plot(wlev_2014.time,wlev_2014.wlev,'g', label='2014 water levels')
ax.plot([s,e],[mean_2014,mean_2014],'-b',lw=2, label='2014 mean')
ax.plot([s,e],[mean_10yr,mean_10yr],'--k',lw=2,label='10yr mean')
ax.set_xlabel('metres')
plt.legend()

print 'mean 2014', mean_2014
print 'mean 10yr', mean_10yr

means_ann={}
means_winter={}
means_summer={}

means=[]
for yr in np.arange(2005,2015,1):
    
    st = datetime.datetime(yr, 1,1); st=st.replace(tzinfo=tz.tzutc())
    ed = datetime.datetime(yr,12,31); ed=ed.replace(tzinfo=tz.tzutc())
    mean = wlev_10yr[(wlev_10yr.time <= ed) & (wlev_10yr.time >= st)]['wlev'].mean()
    means.append(mean)
    print '{} Mean: {}'.format(yr, mean)
print 'Cummulative mean:', np.mean(np.mean(means))
means_ann['Tofino']=np.array(means)

means = [] 
for yr in np.arange(2005,2015,1):
    
    st = datetime.datetime(yr, 11,1); st=st.replace(tzinfo=tz.tzutc())
    ed = datetime.datetime(yr+1,1,31); ed=ed.replace(tzinfo=tz.tzutc())
    
    year = wlev_10yr[(wlev_10yr.time <= ed) & (wlev_10yr.time >= st)]
    
    winter_mean = year['wlev'].mean()
    means.append(winter_mean)
    print '{} to {}  Mean: {}'.format(st.strftime('%d-%b-%Y'), ed.strftime('%d-%b-%Y'),winter_mean)
    
print 'Cummulative Mean:', np.mean(np.array(means))

means_winter['Tofino']=np.array(means)

means = [] 
for yr in np.arange(2005,2015,1):
    
    st = datetime.datetime(yr, 6,1); st=st.replace(tzinfo=tz.tzutc())
    ed = datetime.datetime(yr,8,31); ed=ed.replace(tzinfo=tz.tzutc())
    
    year = wlev_10yr[(wlev_10yr.time <= ed) & (wlev_10yr.time >= st)]
    
    summer_mean = year['wlev'].mean()
    means.append(summer_mean)
    print '{} to {}  Mean: {}'.format(st.strftime('%d-%b-%Y'), ed.strftime('%d-%b-%Y'),summer_mean)
print '10yr Mean:', np.mean(np.array(means))

means_summer['Tofino']=np.array(means)

def get_NOAA(station_no, start_date, end_date, product):
    """Retrieves recent NOAA water levels from a station in a given date range.

    NOAA water levels are at 6 minute intervals and are relative to
    mean sea level.
    See: http://tidesandcurrents.noaa.gov/stations.html?type=Water+Levels.

    :arg station_no: NOAA station number.
    :type station_no: int

    :arg start_date: The start of the date range; e.g. 01-Jan-2014.
    :type start_date: str

    :arg end_date: The end of the date range; e.g. 02-Jan-2014.
    :type end_date: str

    :returns: DataFrame object (obs) with time and wlev columns,
              among others that are irrelevant.
    """

    # Time range
    st_ar = arrow.Arrow.strptime(start_date, '%d-%b-%Y')
    end_ar = arrow.Arrow.strptime(end_date, '%d-%b-%Y')

    base_url = (
        'http://tidesandcurrents.noaa.gov/api/datagetter?')
    params = {
        'product': product,
        'application': 'NOS.COOPS.TAC.WL',
        'begin_date': st_ar.format('YYYYMMDD'),
        'end_date': end_ar.format('YYYYMMDD'),
        'datum': 'STND',
        'station': str(station_no),
        'time_zone': 'GMT',
        'units': 'metric',
        'format': 'csv',
    }
    response = requests.get(base_url, params=params)

    fakefile = StringIO(response.content)
    try:
        obs = pd.read_csv(
            fakefile, parse_dates=[0], date_parser=figures.dateparse_NOAA)
    except ValueError:
        data = {'Date Time': st_ar.datetime, ' Water Level': float('NaN')}
        obs = pd.DataFrame(data=data, index=[0])
    obs = obs.rename(columns={'Date Time': 'time', ' Water Level': 'wlev'})
    return obs

means = []
for yr in np.arange(2005,2015,1):
    
    st = datetime.datetime(yr, 1,1); st=st.replace(tzinfo=tz.tzutc())
    ed = datetime.datetime(yr,12,31); ed=ed.replace(tzinfo=tz.tzutc())
    
    wlev_NB =get_NOAA(figures.SITES['Neah Bay']['stn_no'],st.strftime('%d-%b-%Y'),ed.strftime('%d-%b-%Y'),'hourly_height')
    mean = wlev_NB.wlev.mean()
    means.append(mean)
    print '{} Mean: {}'.format(yr, mean)

print '10yr mean', np.mean(np.array(means))

means_ann['Neah Bay']=np.array(means)

means = []
for yr in np.arange(2005,2015,1):
    
    st = datetime.datetime(yr, 11,1); st=st.replace(tzinfo=tz.tzutc())
    ed = datetime.datetime(yr+1,1,31); ed=ed.replace(tzinfo=tz.tzutc())
    
    wlev_NB =get_NOAA(figures.SITES['Neah Bay']['stn_no'],st.strftime('%d-%b-%Y'),ed.strftime('%d-%b-%Y'),'hourly_height')
    mean = wlev_NB.wlev.mean()
    means.append(mean)
    print '{} to {}  Mean: {}'.format(st.strftime('%d-%b-%Y'), ed.strftime('%d-%b-%Y'),mean)
print 'Cummulative mean', np.mean(np.array(means))

means_winter['Neah Bay'] = np.array(means)

means = []
for yr in np.arange(2005,2015,1):
    
    st = datetime.datetime(yr, 6,1); st=st.replace(tzinfo=tz.tzutc())
    ed = datetime.datetime(yr,8,31); ed=ed.replace(tzinfo=tz.tzutc())
    
    wlev_NB =get_NOAA(figures.SITES['Neah Bay']['stn_no'],st.strftime('%d-%b-%Y'),ed.strftime('%d-%b-%Y'),'hourly_height')
    mean = wlev_NB.wlev.mean()
    means.append(mean)
    print '{} to {}  Mean: {}'.format(st.strftime('%d-%b-%Y'), ed.strftime('%d-%b-%Y'),mean)
print 'Cummulative mean', np.mean(np.array(means))

means_summer['Neah Bay']=np.array(means)

fig,axs=plt.subplots(3,1,figsize=(10,10))


yrs = np.arange(2005,2015,1)
for key, c in zip(['Tofino','Neah Bay'], ['b','g']):
    
    ax=axs[0]
    ax.plot(yrs,means_ann[key], c=c,label =key)
    ax.plot([yrs[0],yrs[-1]],[np.mean(means_ann[key]), np.mean(means_ann[key])], '--',c=c)
    ax.set_title('Annual Mean Water Level')
 
    ax=axs[1]
    ax.plot(yrs,means_winter[key], c=c,label =key)
    ax.plot([yrs[0],yrs[-1]],[np.mean(means_winter[key]), np.mean(means_winter[key])], '--',c=c)
    ax.set_title('Nov-Jan Mean water level')

    ax=axs[2]
    ax.plot(yrs,means_summer[key], c=c,label =key)
    ax.plot([yrs[0],yrs[-1]],[np.mean(means_summer[key]), np.mean(means_summer[key])], '--',c=c)
    ax.set_title('Jun-Aug Mean water level')


x_formatter = matplotlib.ticker.ScalarFormatter(useOffset=False)

for ax in axs:
    ax.xaxis.set_major_formatter(x_formatter)    
    ax.set_ylim([1.7,2.4])
    ax.legend()

from salishsea_tools.nowcast import residuals

t_orig=datetime.datetime(2014,1,1)
t_final = datetime.datetime(2014,12,31)


start_date = t_orig.strftime('%d-%b-%Y')
end_date = t_final.strftime('%d-%b-%Y')
stn_no = figures.SITES['Neah Bay']['stn_no']
obs = get_NOAA(stn_no, start_date, end_date,'hourly_height')
tides = get_NOAA(stn_no, start_date, end_date, 'predictions')
res_obs_NB = residuals.calculate_residual(obs.wlev, obs.time, tides[' Prediction'], tides.time)

fig,ax = plt.subplots(1,1,figsize=(20,5))
ax.plot(obs.time,res_obs_NB)
mean = res_obs_NB.mean()

#ax.plot([obs.time.index[0], obs.time.index[-1]], [mean,mean],'--b')
ax.set_title('Time series for forcing and observed residuals')
print 'Annual mean residual', mean

from salishsea_tools.nowcast import residuals

t_orig=datetime.datetime(2013,1,1)
t_final = datetime.datetime(2013,12,31)


start_date = t_orig.strftime('%d-%b-%Y')
end_date = t_final.strftime('%d-%b-%Y')
stn_no = figures.SITES['Neah Bay']['stn_no']
obs = get_NOAA(stn_no, start_date, end_date,'hourly_height')
tides = get_NOAA(stn_no, start_date, end_date, 'predictions')
res_obs_NB = residuals.calculate_residual(obs.wlev, obs.time, tides[' Prediction'], tides.time)

fig,ax = plt.subplots(1,1,figsize=(20,5))
ax.plot(obs.time,res_obs_NB)
mean = res_obs_NB.mean()
#ax.plot([obs.time.index[0], obs.time.index[-1]], [mean,mean],'--b')
ax.set_title('Time series for forcing and observed residuals')
print 'Annual mean residual', mean

res={}

means = []
for yr in np.arange(2005,2015,1):
    
    st = datetime.datetime(yr, 1,1); st=st.replace(tzinfo=tz.tzutc())
    ed = datetime.datetime(yr,12,31); ed=ed.replace(tzinfo=tz.tzutc())
    
    
    start_date = st.strftime('%d-%b-%Y')
    end_date = ed.strftime('%d-%b-%Y')
    stn_no = figures.SITES['Neah Bay']['stn_no']
    obs = get_NOAA(stn_no, start_date, end_date,'hourly_height')
    tides = get_NOAA(stn_no, start_date, end_date, 'predictions')
    res_obs_NB = residuals.calculate_residual(obs.wlev, obs.time, tides[' Prediction'], tides.time)

    mean = np.mean(res_obs_NB)
    means.append(mean)
    print '{} to {}  Mean: {}'.format(st.strftime('%d-%b-%Y'), ed.strftime('%d-%b-%Y'),mean)
print 'Cummulative mean', np.mean(np.array(means))

res['annual'] = means

means = []
for yr in np.arange(2005,2015,1):
    
    st = datetime.datetime(yr, 11,1); st=st.replace(tzinfo=tz.tzutc())
    ed = datetime.datetime(yr+1,1,31); ed=ed.replace(tzinfo=tz.tzutc())
    
    
    start_date = st.strftime('%d-%b-%Y')
    end_date = ed.strftime('%d-%b-%Y')
    stn_no = figures.SITES['Neah Bay']['stn_no']
    obs = get_NOAA(stn_no, start_date, end_date,'hourly_height')
    tides = get_NOAA(stn_no, start_date, end_date, 'predictions')
    res_obs_NB = residuals.calculate_residual(obs.wlev, obs.time, tides[' Prediction'], tides.time)

    mean = np.mean(res_obs_NB)
    means.append(mean)
    print '{} to {}  Mean: {}'.format(st.strftime('%d-%b-%Y'), ed.strftime('%d-%b-%Y'),mean)
print 'Cummulative mean', np.mean(np.array(means))

res['winter'] = means

means = []
for yr in np.arange(2005,2015,1):
    
    st = datetime.datetime(yr, 6,1); st=st.replace(tzinfo=tz.tzutc())
    ed = datetime.datetime(yr,8,31); ed=ed.replace(tzinfo=tz.tzutc())
    
    
    start_date = st.strftime('%d-%b-%Y')
    end_date = ed.strftime('%d-%b-%Y')
    stn_no = figures.SITES['Neah Bay']['stn_no']
    obs = get_NOAA(stn_no, start_date, end_date,'hourly_height')
    tides = get_NOAA(stn_no, start_date, end_date, 'predictions')
    res_obs_NB = residuals.calculate_residual(obs.wlev, obs.time, tides[' Prediction'], tides.time)

    mean = np.mean(res_obs_NB)
    means.append(mean)
    print '{} to {}  Mean: {}'.format(st.strftime('%d-%b-%Y'), ed.strftime('%d-%b-%Y'),mean)
print 'Cummulative mean', np.mean(np.array(means))

res['summer'] = means

fig,axs=plt.subplots(3,1,figsize=(10,10))


yrs = np.arange(2005,2015,1)
key='Neah Bay'

ax=axs[0]
ax.plot(yrs,means_ann[key], c='b',label ='mean wlev')
ax.plot(yrs,res['annual'], c='g',label ='residual')
ax.set_title('Annual Means')
 
ax=axs[1]
ax.plot(yrs,means_winter[key], c='b',label ='mean wlev')
ax.plot(yrs,res['winter'], c='g',label ='residual')
ax.set_title('Winter Means')

ax=axs[2]
ax.plot(yrs,means_summer[key], c='b',label ='mean wlev')
ax.plot(yrs,res['summer'], c='g',label ='residual')
ax.set_title('Summer Means')

x_formatter = matplotlib.ticker.ScalarFormatter(useOffset=False)

for ax in axs:
    ax.xaxis.set_major_formatter(x_formatter)    
    #ax.set_ylim([1.7,2.4])
    ax.legend()

fig,axs=plt.subplots(3,1,figsize=(10,10))


yrs = np.arange(2005,2015,1)
key='Neah Bay'

ax=axs[0]
ax.plot(yrs,means_ann[key] - np.mean(means_ann[key]), c='b',label ='mean wlev')
ax.plot(yrs,res['annual']-np.mean(res['annual']), c='g',label ='residual')
ax.set_title('Annual Means')
 
ax=axs[1]
ax.plot(yrs,means_winter[key] - np.mean(means_winter[key]), c='b',label ='mean wlev')
ax.plot(yrs,res['winter']-  np.mean(res['winter']), c='g',label ='residual')
ax.set_title('Winter Means')

ax=axs[2]
ax.plot(yrs,means_summer[key]- np.mean(means_summer[key]), c='b',label ='mean wlev')
ax.plot(yrs,res['summer'] - np.mean(res['summer']), c='g',label ='residual')
ax.set_title('Summer Means')

x_formatter = matplotlib.ticker.ScalarFormatter(useOffset=False)

for ax in axs:
    ax.xaxis.set_major_formatter(x_formatter)    
    #ax.set_ylim([1.7,2.4])
    ax.legend(loc=0)
    ax.grid()


start = '01-Jan-2005'; end = '31-Dec-2014'

wlev_10yr = figures.load_archived_observations('Point Atkinson',start,end )

tides,msl = stormtools.load_tidal_predictions('Point Atkinson_atide_compare8_01-Jan-2005_31-Dec-2018.csv')

mean_10yr = wlev_10yr.wlev.mean()

means=[]
for yr in np.arange(2005,2015,1):
    
    st = datetime.datetime(yr, 1,1); st=st.replace(tzinfo=tz.tzutc())
    ed = datetime.datetime(yr,12,31); ed=ed.replace(tzinfo=tz.tzutc())
    mean = wlev_10yr[(wlev_10yr.time <= ed) & (wlev_10yr.time >= st)]['wlev'].mean()
    means.append(mean)
    print '{} Mean: {}'.format(yr, mean)
print 'Cummulative mean:', np.mean(np.mean(means))
means_ann['Point Atkinson']=np.array(means)

means=[]
means_res={}
msl=figures.SITES['Point Atkinson']['msl']
for yr in np.arange(2005,2015,1):
    
    st = datetime.datetime(yr, 1,1); st=st.replace(tzinfo=tz.tzutc())
    ed = datetime.datetime(yr,12,31); ed=ed.replace(tzinfo=tz.tzutc())
    
    wlev = wlev_10yr[(wlev_10yr.time <= ed) & (wlev_10yr.time >= st)]['wlev']
    time = wlev_10yr[(wlev_10yr.time <= ed) & (wlev_10yr.time >= st)]['time']
    
    residual = residuals.calculate_residual(np.array(wlev),np.array(time), np.array(tides.pred_all)+msl, np.array(tides.time))
    mean = np.mean(residual)

    means.append(mean)
    print '{} Mean: {}'.format(yr, mean)
print 'Cummulative mean:', np.mean(np.mean(means))
means_res['Point Atkinson']=np.array(means)

fig,axs=plt.subplots(1,1,figsize=(10,5))


yrs = np.arange(2005,2015,1)
key='Point Atkinson'

ax=axs
ax.plot(yrs,means_ann[key] - np.mean(means_ann[key]), c='b',label ='mean wlev PA')
ax.plot(yrs,means_res[key]-np.mean(means_res[key]), c='g',label ='residual PA')

key = 'Neah Bay'
ax.plot(yrs,means_ann[key] - np.mean(means_ann[key]), '--r',label ='mean wlev NB')
ax.plot(yrs,res['annual']-np.mean(res['annual']), '--k',label ='residual NB')
ax.set_title('Annual Means')
 

x_formatter = matplotlib.ticker.ScalarFormatter(useOffset=False)
ax.xaxis.set_major_formatter(x_formatter)    
    #ax.set_ylim([1.7,2.4])
ax.legend(loc=0)
ax.grid()
ax.set_ylabel('[m]')