# -*- coding: utf-8 -*-
"""
Created on Tue Jan 16 15:47:46 2024

@author: chuffard
"""

# -*- coding: utf-8 -*-
"""
Created on Wed Dec 13 09:34:42 2023

@author: chuffard
"""

import numpy as np

import pandas as pd

#import spectrum #doesn't want to install...tried everything I can find

import matplotlib.pyplot as plt

#from spectrum import Periodogram, data_cosine

from datetime import datetime

import scipy as SP

import scipy.signal as signal

import seaborn as sns

#import endaq as endaq

import statsmodels.tsa.stattools as sts

from statsmodels.tsa.stattools import acf

import pyflakes as pyf

#import lzip as lz

from pandas import read_csv

from statsmodels.graphics.tsaplots import plot_acf

from statsmodels.graphics.tsaplots import plot_pacf


import ephem

import pygam

import re

from matplotlib.backends.backend_pdf import PdfPages

mako15 = sns.color_palette("rocket", 15)
mako60 = sns.color_palette("rocket", 60)
mako14 = sns.color_palette("rocket", 14)
sdiverging_colors = sns.color_palette("RdBu", 24)
sdiverging_colors80 = sns.color_palette("RdBu", 80)
#This is the time period without the crab

# read in data; use first line as header
Main_Hourly_SES_MARS = pd.read_csv('G:/Analyses/Main_Hourly_SES_MARS.csv', parse_dates=True, index_col="datetime")

Main_Hourly_SES_MARS["atn"].plot()
#remove infinite values
Main_Hourly_SES_MARS= Main_Hourly_SES_MARS[Main_Hourly_SES_MARS['atn']!= np.inf]
list(Main_Hourly_SES_MARS.columns)

# ['atn',
#  'mooninterp',
#  'hour',
#  'year',
#  'month',
#  'day',
#  'Year-Week',
#  'justdate',
#  'ABSdayssincesolst',
#  'moonbins',
#  'Pred',
#  'tidal_diff',
#  'dayssincesolstice',
#  'NUMdayssincesolst',
#  'sunrisedecimalhour',
#  'sunsetdecimalhour',
#  'sunrise',
#  'sunset']

#plot time series of eac column

Main_Hourly_SES_MARS.plot(subplots=True, layout=(18,1), figsize=(15,15))
plt.savefig('all data hourly.pdf')




#Try autocorrelation
#from https://machinelearningmastery.com/gentle-introduction-autocorrelation-partial-autocorrelation/
#Trim to period with no big gaps
trimhourly = Main_Hourly_SES_MARS.loc['2023-06-27':'2023-08-14']  
trimhourly.plot(subplots=True, layout=(15,1), figsize=(15,15))
plt.savefig('TRIMMED for autocorrelation all data hourly.pdf')

plot_acf(trimhourly['atn'], lags=60,missing = 'conservative')
plt.title('hourly, autocorrelation atn')
plt.xlabel('lag (hours)')
plt.savefig('atn acf hourly.pdf')

#try spectray density welch
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html
fs = 1
f, Pxx_den = signal.welch(trimhourly['atn'])
plt.subplot(2, 1, 1)
plt.semilogy(f, Pxx_den)
plt.ylim([0.5e-3, 1])
plt.xlabel('frequency [units]')
plt.ylabel('PSD [V**2/units]')
plt.title('hourly, trimmed to no crab period, welch')
plt.savefig('welch acf hourly.pdf')

#Now try this with the tides 
#Try autocorrelation
#from https://machinelearningmastery.com/gentle-introduction-autocorrelation-partial-autocorrelation/
#Trim to period with no big gaps

plot_acf(trimhourly['Pred'], lags=60,missing = 'conservative')
plt.xlabel('lag (hours)')
plt.title('hourly, autocorrelation tidal height')
plt.savefig('tidal height acf hourly.pdf')

#try spectray density welch
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html
fs = 1
f, Pxx_den = signal.welch(trimhourly['Pred'])
plt.subplot(2, 1, 1)
plt.semilogy(f, Pxx_den)
plt.ylim([0.5e-3, 1])
plt.xlabel('frequency [units]')
plt.ylabel('PSD [V**2/units]')
plt.title('hourly, tidal height, welch')
plt.savefig('welch tidal height hourly.pdf')

#Now try this with the tidal change
#Try autocorrelation
#from https://machinelearningmastery.com/gentle-introduction-autocorrelation-partial-autocorrelation/
#Trim to period with no big gaps

plot_acf(trimhourly['tidal_diff'], lags=60,missing = 'conservative')
plt.title('hourly, autocorrelation tidal change')
plt.xlabel('lag (hours)')
plt.savefig('tidal change acf hourly.pdf')
#try spectray density welch
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html
fs = 1
f, Pxx_den = signal.welch(trimhourly['tidal_diff'])
plt.subplot(2, 1, 1)
plt.semilogy(f, Pxx_den)
plt.ylim([0.5e-3, 1])
plt.xlabel('frequency [units]')
plt.ylabel('PSD [V**2/units]')
plt.title('hourly, tidal change, welch')
plt.savefig('welch tidal change hourly.pdf')


#COMPLETE CASES FROM HERE ON OUT#########
#rfilter to complete cases
Main_Hourly_SES_MARSComplete= Main_Hourly_SES_MARS.loc[Main_Hourly_SES_MARS.notna().all(axis='columns')]

#plot by hour of day
hourly_atn = Main_Hourly_SES_MARSComplete[['hour', 'atn']] 
hourly_atn= hourly_atn.loc[hourly_atn.notna().all(axis='columns')]
boxplot = hourly_atn.boxplot(by='hour')
plt.savefig('bot plots atn by hour.pdf')

#plot by hour of day
hourly_tideheight = Main_Hourly_SES_MARSComplete[['hour', 'Pred']] 
hourly_tideheight= hourly_tideheight.loc[hourly_tideheight.notna().all(axis='columns')]
boxplot = hourly_tideheight.boxplot(by='hour')
#why is this slightly different in the midnight-2 am hours from the last one?
plt.savefig('bot plots tidal height by hour.pdf')

#plot by hour of day
hourly_tidal_diff = Main_Hourly_SES_MARSComplete[['hour', 'tidal_diff']] 
hourly_tidal_diff= hourly_tidal_diff.loc[hourly_tidal_diff.notna().all(axis='columns')]
boxplot = hourly_tidal_diff.boxplot(by='hour')
#why is this slightly different in the midnight-2 am hours from the last one?
plt.savefig('bot plots hourly_tidal_diff by hour.pdf')


sns.relplot(data=Main_Hourly_SES_MARSComplete, x="hour", y="atn", kind="line")
plt.savefig('atn by 24 hour ribbon.pdf')

sns.relplot(data=Main_Hourly_SES_MARSComplete, x="hour", y="Pred", kind="line")
plt.savefig('tidal height by 24 hour ribbon.pdf')

sns.relplot(data=Main_Hourly_SES_MARSComplete, x="hour", y="tidal_diff", kind="line")
plt.savefig('tidal_diff by 24 hour ribbon.pdf')

#technicaolly works but looks terrible
# fig, (ax1, ax2, ax3) = plt.subplots(1, 3, subplot_kw=dict(projection='polar'))
# ax1.plot(Main_Hourly_SES_MARSComplete['hour'], Main_Hourly_SES_MARSComplete['atn'])
# ax2.plot(Main_Hourly_SES_MARSComplete['hour'], Main_Hourly_SES_MARSComplete['Pred'])
# ax3.plot(Main_Hourly_SES_MARSComplete['hour'],Main_Hourly_SES_MARSComplete['tidal_diff'])
# plt.show()

minsunrise = Main_Hourly_SES_MARSComplete['sunrisedecimalhour'].min()
maxsunrise = Main_Hourly_SES_MARSComplete['sunrisedecimalhour'].max()
minsunset = Main_Hourly_SES_MARSComplete['sunsetdecimalhour'].min()
maxsunset = Main_Hourly_SES_MARSComplete['sunsetdecimalhour'].max()

fig, axes = plt.subplots(3, 1, figsize=(10, 10))
for ax in axes:
    ax.axvspan(minsunrise, maxsunrise, color=sns.xkcd_rgb['grey'], alpha=0.5)
    ax.axvspan(minsunset, maxsunset,  color=sns.xkcd_rgb['grey'], alpha=0.5)
sns.lineplot(ax=axes[2], x="hour", y="tidal_diff", data=Main_Hourly_SES_MARSComplete)
sns.lineplot(ax=axes[0], x="hour", y="atn", data=Main_Hourly_SES_MARSComplete)
sns.lineplot(ax=axes[1], x="hour", y="Pred", data=Main_Hourly_SES_MARSComplete)
axes[0].set(xlabel=None)
axes[1].set(xlabel=None)
axes[0].set(ylabel="attenuance")
axes[1].set(ylabel="tidal height (m)")
axes[2].set(ylabel='change in tidal height (m/h)')
plt.savefig('24 hour ribbon.pdf')
# a = sns.lineplot(data=useatnmoonphaseCRAB, 
#              x="hour", 
#              y="atn", 
#              hue = 'Year-Week',
#              palette = mako14, ci = None, ax=ax[0])
fig, axes = plt.subplots(3, 1, figsize=(10, 10))
for ax in axes:
    ax.axvspan(minsunrise, maxsunrise, color=sns.xkcd_rgb['grey'], alpha=0.5)
    ax.axvspan(minsunset, maxsunset,  color=sns.xkcd_rgb['grey'], alpha=0.5)
sns.lineplot(ax=axes[0], x="hour", y="atn", data=Main_Hourly_SES_MARSComplete, hue = 'NUMdayssincesolst',palette = 'RdBu')
sns.lineplot(ax=axes[1], x="hour", y="Pred", data=Main_Hourly_SES_MARSComplete, hue = 'NUMdayssincesolst',palette = 'RdBu')
sns.lineplot(ax=axes[2], x="hour", y="tidal_diff", data=Main_Hourly_SES_MARSComplete, hue = 'NUMdayssincesolst',palette = 'RdBu')
axes[0].get_legend().remove()
axes[1].get_legend().remove()
axes[0].set(xlabel=None)
axes[1].set(xlabel=None)
#sns.move_legend(axes[2], "lower center")
plt.legend(bbox_to_anchor =(0.5,-0.3), loc='lower center',mode = "expand", ncol = 4)
plt.tight_layout()
plt.savefig('24 hour by days from solstice.pdf')

#moon phase bins for coloring that line plot
bins = [0, .25, .5, .75, 1]
Main_Hourly_SES_MARSComplete['moonbins'] = pd.cut(Main_Hourly_SES_MARSComplete['mooninterp'], bins)

mako4 = sns.color_palette("rocket", 4)
sns.lineplot(data=Main_Hourly_SES_MARSComplete, 
             x="hour", 
             y="atn", 
             hue = 'moonbins',
             palette = mako4)
plt.legend(loc=(1.04, 0))
plt.savefig('atn by 24 hour no crab one panel moon bins CI.pdf')

mako4 = sns.color_palette("rocket", 4)
sns.lineplot(data=Main_Hourly_SES_MARSComplete, 
             x="hour", 
             y="atn", 
             hue = 'moonbins',
             palette = mako4,
             ci = None)
plt.legend(loc=(1.04, 0))
plt.savefig('atn by 24 hour no crab one panel moon bins CI none.pdf')


sns.lineplot(data=Main_Hourly_SES_MARSComplete, 
             x="hour", 
             y="atn", 
             hue = 'Year-Week',
             palette = mako14, ci = None)
plt.legend(loc=(1.04, 0))
plt.savefig('atn by 24 hour no crab one panel.pdf', bbox_inches='tight')

#CONTOUR MAPS
plt.close('all')
plt.tricontourf(Main_Hourly_SES_MARSComplete["hour"], Main_Hourly_SES_MARSComplete["mooninterp"], Main_Hourly_SES_MARSComplete["atn"],palette = 'rocket')
plt.title('Attenuance by hour and lunar illumination ')
plt.xlabel('hour')
plt.ylabel('lunar illumination')
plt.savefig('atn by 24 hour.pdf')

plt.close('all')
plt.tricontourf(Main_Hourly_SES_MARSComplete["hour"], Main_Hourly_SES_MARSComplete["Pred"], Main_Hourly_SES_MARSComplete["atn"],palette = 'rocket')
plt.title('Attenuance by hour and tidal height')
plt.xlabel('hour')
plt.ylabel('tidal height (m)')
plt.savefig('atn by 24 hour tidal height.pdf')


plt.close('all')
plt.tricontourf(Main_Hourly_SES_MARSComplete["hour"], Main_Hourly_SES_MARSComplete["tidal_diff"], Main_Hourly_SES_MARSComplete["atn"],palette = 'rocket')
plt.title('Attenuance by hour and tidal_diff (m/h)')
plt.xlabel('hour')
plt.ylabel('tidal_diff (m/h)')
plt.savefig('atn by 24 hour tidal_diff .pdf')

##Try gam
#https://gist.github.com/josef-pkt/453de603b019143e831fbdd4dfb6aa30