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

@author: chuffard
"""

#Data prep and exploration for SES 

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



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

#remove infinite values
SESdata= SESdata[SESdata['atn']!= np.inf]

SESdata.dtypes
SEStrim = SESdata.drop(['Date', 'Time'], axis=1)

SEStrim["atn"].plot()

#find hourly average of SES time series
hourly_SESatn = SEStrim['atn'].resample('1H').mean()

#transpose
hourly_SESatndf = pd.DataFrame([hourly_SESatn]).T

#fix weird time zone stuff
for col in hourly_SESatndf.select_dtypes(include=['datetime64[ns, UTC]']).columns:
    hourly_SESatndf[col] = hourly_SESatndf[col].apply(lambda x: x.tz_localize(None))
    
#Trim to period w prior to crab clog
trimhourly = hourly_SESatndf.loc['2023-06-27':'2023-08-14']  
trimhourly.plot()
plt.title('hourly, trimmed to no crab period')
plt.ylabel('attenuance')

#fill na with interpolation
timehourlyinterpolate = trimhourly.interpolate(method='linear')



#Trim to period w prior to crab clog
nocrabhourly = hourly_SESatndf.loc[:'2023-08-14']  
nocrabhourly.plot()
plt.title('hourly, trimmed to no crab period')
plt.ylabel('attenuance')

#can't inperpolate the no crab period because there are big gaps

#try spectray density welch
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html
fs = 1
f, Pxx_den = signal.welch(timehourlyinterpolate['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.show()

#Try autocorrelation
#from https://machinelearningmastery.com/gentle-introduction-autocorrelation-partial-autocorrelation/

plot_acf(timehourlyinterpolate['atn'], lags=60)
plt.title('hourly, trimmed to no crab period, autocorrelation')

#adjusting N doesn't make a difference
#fft foesn't make a difference
#drop na no difference

plot_acf(timehourlyinterpolate['atn'], lags=60, missing = 'drop')
plt.title('hourly, trimmed to no crab period, autocorrelation drop na')

#adjusting N doesn't make a difference

#These are the acf results for hourly averages
acfarray = acf(timehourlyinterpolate['atn'], adjusted=False, 
                   nlags=60, qstat=True, fft=False, alpha=0.05, 
                   bartlett_confint=True, missing='none')

results = pd.DataFrame()
results["acf"] = acf(timehourlyinterpolate['atn'], nlags=60)

results.to_csv("C:/Python_temp/acfhourly.csv")


plot_pacf(timehourlyinterpolate['atn'], lags=60)
plt.title('hourly, trimmed to no crab period, partial autocorrelation')



#plot by hours of the day
timehourlyinterpolate.index.dtype #dtype('<M8[ns]')
timehourlyinterpolate['datetime'] = timehourlyinterpolate.index.strftime('%Y-%m-%d %H:%M:%S')
timehourlyinterpolate['datetime']  = pd.to_datetime(timehourlyinterpolate['datetime'] , format='%Y-%m-%d %H:%M:%S')
timehourlyinterpolate.datetime.dtype #dtype('O')
timehourlyinterpolate['hour'] = timehourlyinterpolate['datetime'].dt.hour


#nocrabhourly
nocrabhourly.index.dtype #dtype('<M8[ns]')
nocrabhourly['datetime'] = nocrabhourly.index.strftime('%Y-%m-%d %H:%M:%S')
nocrabhourly['datetime']  = pd.to_datetime(nocrabhourly['datetime'] , format='%Y-%m-%d %H:%M:%S')
nocrabhourly.datetime.dtype #dtype('O')
nocrabhourly['hour'] = nocrabhourly['datetime'].dt.hour

boxplot = timehourlyinterpolate.boxplot(by='hour')
plt.savefig('bot plots atn by hour.eps',format='eps')

#nocrabhourly
boxplot = nocrabhourly.boxplot(by='hour')
plt.title('no crab bot plots atn by hour')
plt.savefig('no crab bot plots atn by hour.pdf',)


sns.relplot(data=timehourlyinterpolate, x="hour", y="atn", kind="line")
plt.savefig('atn by 24 hour.eps',format='eps')

#nocrabhourly
sns.relplot(data=nocrabhourly, x="hour", y="atn", kind="line")
plt.title('no crab bot plots atn by hour')
plt.savefig('no crab atn by 24 hour.pdf')

#plot all against attenuation time series

#moon phase illumination
nocrabhourly['year'] = nocrabhourly['datetime'].dt.year
nocrabhourly['month'] = nocrabhourly['datetime'].dt.month
nocrabhourly['day'] = nocrabhourly['datetime'].dt.day
nocrabhourly['simpleEdatetime'] = nocrabhourly['year'].map(str) +'/'+ nocrabhourly['month'].map(str) + '/' + nocrabhourly['day'].map(str)



#how many days in sample 8/14/2023 - 6/27/2023
n = nocrabhourly.nunique(axis=0)
print(n) #simpleEdatetime

m = ephem.Moon()
m.compute('2023/05/12')
j = ephem.Date('2023/05/12')
for add in range(0, 96):
    print (ephem.date(j + add))
    print(m.moon_phase)



b = pd.DataFrame()
c = []
e=[]
l=[]
#needs padding of the date for interpolating the moon phase
j = ephem.Date('2023/05/12')
for add in range(0, 96):
    d = (ephem.date(j + add))
    m = ephem.Moon(d)
    moonphase = m.moon_phase
    c.append(moonphase)
    print(ephem.date(d + ephem.hour))
    local = ephem.localtime(d)
    e.append(d)
    l.append(local)
    
moonphasedata = pd.DataFrame(
    {'local_time': l,
     'UTCtime': e,
     'moonphase': c
    })


CRABtimehourly = nocrabhourly.copy()

#Needs to get into same time zone. this says UTC but it's actually local
CRABtimehourly['local_time'] = CRABtimehourly['datetime']
CRABtimehourly['UTC'] = pd.to_datetime(CRABtimehourly['local_time'], utc = True)
moonphasedata['UTC'] = pd.to_datetime(moonphasedata['local_time'], utc = True)
atnmoonphase = pd.concat([CRABtimehourly, moonphasedata])
#this does funky things to the dates and times, so aggregate by datetime
useatnmoonphaseCRAB = pd.DataFrame(atnmoonphase, columns=['local_time', 'atn','moonphase'])

useatnmoonphaseCRAB = useatnmoonphaseCRAB.resample('1H', on='local_time').mean()

#now turn index to column
#turn index into a column
useatnmoonphaseCRAB['datetime'] = useatnmoonphaseCRAB.index

#now interpolate moonphase
useatnmoonphaseCRAB['mooninterp'] = useatnmoonphaseCRAB['moonphase'].interpolate(method='linear')


#plot moon phase against interpolated moon phase to make sure it's working

ax1 = useatnmoonphaseCRAB.plot(x = 'datetime', y = ["mooninterp"])
ax2 = useatnmoonphaseCRAB.plot.scatter(x = 'datetime', y = ['moonphase'], marker = '^', ax = ax1)
plt.xticks(rotation=45);

#now remove the original moonphase because of all the NAs
useatnmoonphaseCRAB= useatnmoonphaseCRAB.drop(columns=['moonphase'])

#now bring back the hour of the day
from datetime import date
useatnmoonphaseCRAB['hour'] = useatnmoonphaseCRAB['datetime'].dt.hour
useatnmoonphaseCRAB['year'] = useatnmoonphaseCRAB['datetime'].dt.year
useatnmoonphaseCRAB['month'] = useatnmoonphaseCRAB['datetime'].dt.month
useatnmoonphaseCRAB['day'] = useatnmoonphaseCRAB['datetime'].dt.day
useatnmoonphaseCRAB['Year-Week'] = useatnmoonphaseCRAB['datetime'].dt.strftime('%Y-%U')
useatnmoonphaseCRAB['justdate'] = useatnmoonphaseCRAB['datetime'].dt.date
useatnmoonphaseCRAB['dayssincesolstice'] = useatnmoonphaseCRAB['justdate'] - date(2023,6,21)
#turn this into a string
useatnmoonphaseCRAB['daysstr'] = useatnmoonphaseCRAB['dayssincesolstice'].astype(str)
#now extract the values 
useatnmoonphaseCRAB['result'] = useatnmoonphaseCRAB['daysstr'].str.extract(r'(\d+)', expand=False)
useatnmoonphaseCRAB.dtypes
useatnmoonphaseCRAB['NUMdayssincesolstice'] = pd.to_numeric(useatnmoonphaseCRAB['result'])


mako4 = sns.color_palette("rocket", 4)
mako9 = sns.color_palette("rocket", 9)
sdiverging_colors = sns.color_palette("RdBu", 24)
#rfilter to complete cases
useatnmoonphaseCRAB= useatnmoonphaseCRAB.loc[useatnmoonphaseCRAB.notna().all(axis='columns')]

# #now to get the two weeks column
# len(useatnmoonphaseCRAB) #1176
# len(useatnmoonphaseCRAB) /336
# twoweeks = 336 # number of times you repeat each item




sns.scatterplot(data=useatnmoonphaseCRAB, x='mooninterp', y='atn', hue='month')
sns.scatterplot(data=useatnmoonphaseCRAB, x='hour', y='atn', hue='mooninterp')

plt.tricontourf(useatnmoonphaseCRAB["hour"], useatnmoonphaseCRAB["mooninterp"], useatnmoonphaseCRAB["atn"],palette = 'rocket')
plt.title('Attenuance by hour and lunar illumination no crab')
plt.xlabel('hour')
plt.ylabel('lunar illumination')
plt.savefig('atn by 24 hour no crab.pdf')


plt.tricontourf(useatnmoonphaseCRAB["hour"], useatnmoonphaseCRAB["NUMdayssincesolstice"], useatnmoonphaseCRAB["atn"],palette = 'rocket')
plt.title('Attenuance by hour and lunar illumination no crab')
plt.xlabel('hour')
plt.ylabel('days since solstice')
plt.savefig('atn by 24 hour solstice no crab.pdf')





mako60 = sns.color_palette("rocket", 60)
mako14 = sns.color_palette("rocket", 14)
sdiverging_colors = sns.color_palette("RdBu", 24)
# instanciate figure and axes instances
fig,ax = plt.subplots(1,2,figsize=(6, 4)) 
a = sns.lineplot(data=useatnmoonphaseCRAB, 
             x="hour", 
             y="atn", 
             hue = 'Year-Week',
             palette = mako14, ci = None, ax=ax[0])
ax[0].legend_.remove()
b = sns.scatterplot(data = useatnmoonphaseCRAB,
               x = 'hour', 
               y = 'atn',
               hue = 'Year-Week',
               palette = mako14,ax=ax[1])
b.set(ylabel=None)
# sns.move_legend(
#     ax[1],
#     loc="right",
#     bbox_to_anchor=(-0.1, 1.1),
#     ncol=2,
#     title=None,
#     frameon=False,
# )
sns.move_legend(ax[1], "upper left", bbox_to_anchor=(1, 1))
ax[0].autoscale()
ax[1].autoscale()
plt.savefig('atn by 24 hour no crab two panels tight.pdf', bbox_inches='tight')




sns.lineplot(data=useatnmoonphaseCRAB, 
             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')



sns.lineplot(data=useatnmoonphaseCRAB, 
             x="hour", 
             y="atn", 
             hue = 'NUMdayssincesolstice',
             palette =mako60, ci = None)
plt.legend(loc=(1.04, 0), ncol=4)
plt.savefig('atn by 24 hour no crab days since solstice.pdf', bbox_inches='tight')




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

mako4 = sns.color_palette("rocket", 4)
sns.lineplot(data=useatnmoonphaseCRAB, 
             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=useatnmoonphaseCRAB, 
             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')


#### Has peak of atn changed with the timing of sunrise and sunset, before and after the solstice?
df= useatnmoonphaseCRAB.copy(deep=True)

MARS      = ephem.Observer()

#Location of Fredericton, Canada
MARS.lon  = str(-122.18689833333333) #Note that lon should be in string format
MARS.lat  = str(36.712468333333334)      #Note that lat should be in string format
sun = ephem.Sun()
#122 degrees 11.2139 minutes
#36 degrees 42.7481
#Elevation of sea level
MARS.elev = 0
#To get U.S. Naval Astronomical Almanac values, use these settings
MARS.pressure= 0
MARS.horizon = '-12' #-6=civil twilight, -12=nautical, -18=astronomical


#### Has peak of atn changed with the timing of sunrise and sunset, before and after the solstice?
MARS      = ephem.Observer()

#Location of Fredericton, Canada
MARS.lon  = str(-122.18689833333333) #Note that lon should be in string format
MARS.lat  = str(36.712468333333334)      #Note that lat should be in string format
sun = ephem.Sun()

def get_time(obs, obj, func):
    func = getattr(obs, func)
    def inner(date):
        obs.date = date
        return ephem.localtime(func(obj))
    return inner

df['rise'] = pd.Series(df.justdate).apply(get_time(MARS, sun, 'next_rising'))
df['set'] = pd.Series(df.justdate).apply(get_time(MARS, sun, 'next_setting'))

#ok that works, now use for main fataframe
useatnmoonphaseCRAB['sunrise'] = pd.Series(useatnmoonphaseCRAB.justdate).apply(get_time(MARS, sun, 'next_rising'))
useatnmoonphaseCRAB['sunset'] = pd.Series(useatnmoonphaseCRAB.justdate).apply(get_time(MARS, sun, 'next_setting'))


useatnmoonphaseCRAB.to_csv("SESMARS_hourly.csv")

#trim to morning and evening datasets- midnight to noon, noon to midnight

nightSESMARS = useatnmoonphaseCRAB.query('hour > 11')
morningSESMARS = useatnmoonphaseCRAB.query('hour < 12')


#now what's the time of the peak each day
idx = nightSESMARS.groupby('justdate')['atn'].idxmax()
max_scores_night = nightSESMARS.loc[idx]
max_scores_night['sunsettime']= max_scores_night['sunset'].dt.time
max_scores_night['atnpeaktime']= max_scores_night['datetime'].dt.time
max_scores_night.to_csv("max_scores_night.csv")







#now what's the time of the peak each day
idx = morningSESMARS.groupby('justdate')['atn'].idxmax()
max_scores_morning = morningSESMARS.loc[idx]
max_scores_morning.to_csv("max_scores_morning.csv")



sns.scatterplot(data=max_scores_morning, x='sunrise', y='hour', hue='mooninterp')
plt.xticks(rotation=45);
plt.legend(loc=(1.04, 0))

sns.scatterplot(data=max_scores_night, x='sunset', y='hour', hue='mooninterp')
plt.xticks(rotation=45);
plt.legend(loc=(1.04, 0))


sns.scatterplot(data=max_scores_night, x='sunset', y='hour', hue='NUMdayssincesolstice')
plt.xticks(rotation=45);
plt.legend(loc=(1.04, 0))

sns.scatterplot(data=max_scores_morning, x='sunset', y='hour', hue='NUMdayssincesolstice')
plt.xticks(rotation=45);
plt.legend(loc=(1.04, 0))


plt.tricontourf(max_scores_morning["hour"], max_scores_morning["NUMdayssincesolstice"], max_scores_morning["atn"],palette = 'rocket')
plt.title('morning Attenuance by hour and lunar illumination no crab')
plt.xlabel('hour')
plt.ylabel('NUMdayssincesolstice')
plt.savefig('morning NUMdayssincesolstice atn by 24 hour no crab.pdf')



plt.tricontourf(useatnmoonphaseCRAB["hour"], useatnmoonphaseCRAB["NUMdayssincesolstice"], useatnmoonphaseCRAB["atn"],palette = 'rocket')
plt.title(' Attenuance by hour and lunar illumination no crab')
plt.xlabel('hour')
plt.ylabel('NUMdayssincesolstice')
plt.savefig(' NUMdayssincesolstice atn by 24 hour no crab.pdf')


sns.scatterplot(data=useatnmoonphaseCRAB, x='hour', y='atn', hue='NUMdayssincesolstice')
plt.legend(loc=(1.04, 0))

#PyGam is useless- gives wrong p values and lambda

#import tides
MossLandingTides = pd.read_csv('G:/Analyses/ASLO Poster/Moss_landing_tides_CSV.csv', parse_dates=[['Date', 'Time']], index_col="Date_Time")
#Find hourly averages for merging with full dataset
hourly_MossLandingTides = MossLandingTides['Pred'].resample('1H').mean()
#fill in tides for every hour
hourly_MossLandingTidesinterp = hourly_MossLandingTides.interpolate(method='linear')
tidesDF = pd.DataFrame(hourly_MossLandingTidesinterp)

sns.lineplot(data=tidesDF, x=tidesDF.index, y=tidesDF['Pred'])
plt.xticks(rotation=45);




#find speed of change between tidal predictions as a proxy for current speed
tidesDF['tidal_diff'] = tidesDF.diff(axis = 0, periods = 1) 
sns.lineplot(data=tidesDF, x=tidesDF.index, y=tidesDF['tidal_diff'])
plt.xticks(rotation=45);
tidesDF = tidesDF.reset_index()
#rename Datae_time to datetime
tidesDF = tidesDF.rename(columns={'Date_Time': 'datetime'})



SESdataHourly = pd.read_csv('G:/Analyses/ASLO poster/SES_Hourly_USE.csv', parse_dates=True, index_col="datetime")
#SESdataHourly['datetime'] = SESdataHourly.index
SESdataHourly = SESdataHourly.reset_index()

#now merge with SES hourly data
Main_Hourly_SES_MARS = SESdataHourly.merge(tidesDF, how='outer',on='datetime')
#now set the index back to the datetime
Main_Hourly_SES_MARS = Main_Hourly_SES_MARS.set_index('datetime')
Main_Hourly_SES_MARS = Main_Hourly_SES_MARS.loc['2023-05-13':'2023-08-14']  

Main_Hourly_SES_MARS.index = pd.to_datetime(Main_Hourly_SES_MARS.index)

#now bring back the hour of the day
from datetime import date
Main_Hourly_SES_MARS['hour'] = Main_Hourly_SES_MARS.index.hour
Main_Hourly_SES_MARS['year'] = Main_Hourly_SES_MARS.index.year
Main_Hourly_SES_MARS['month'] = Main_Hourly_SES_MARS.index.month
Main_Hourly_SES_MARS['day'] = Main_Hourly_SES_MARS.index.day
Main_Hourly_SES_MARS['Year-Week'] = Main_Hourly_SES_MARS.index.strftime('%Y-%U')
Main_Hourly_SES_MARS['justdate'] = Main_Hourly_SES_MARS.index.date
Main_Hourly_SES_MARS['dayssincesolstice'] = Main_Hourly_SES_MARS['justdate'] - date(2023,6,21)
#turn this into a string
Main_Hourly_SES_MARS['daysstr'] = Main_Hourly_SES_MARS['dayssincesolstice'].astype(str)
#now extract the values 
Main_Hourly_SES_MARS['result'] = Main_Hourly_SES_MARS['daysstr'].str.extract(r'(\d+)', expand=False)
Main_Hourly_SES_MARS.dtypes
Main_Hourly_SES_MARS['NUMdayssincesolstice'] = pd.to_numeric(Main_Hourly_SES_MARS['result'])

#delete columns not going to use
Main_Hourly_SES_MARS = Main_Hourly_SES_MARS.drop(columns=['result', 'daysstr'])

Main_Hourly_SES_MARS.to_csv("Main_Hourly_SES_MARS.csv")


# ##Gam
# useatnmoonphaseCRABNOGAPS= useatnmoonphaseCRAB.loc[useatnmoonphaseCRAB.notna().all(axis='columns')]
# useatnmoonphaseCRABNOGAPS = useatnmoonphaseCRABNOGAPS.reset_index()

# list(useatnmoonphaseCRAB.columns)
# # ['atn',
# #  'datetime',
# #  'mooninterp',
# #  'hour',
# #  'year',
# #  'month',
# #  'day',
# #  'Year-Week',
# #  'justdate',
# #  'dayssincesolstice',
# #  'moonbins',
# #  'sunrise',
# #  'sunset']

# #https://www.youtube.com/watch?v=XQ1vk7wEI7c&ab_channel=PyData
# useatnmoonphaseCRABx = useatnmoonphaseCRABNOGAPS[['mooninterp', 'hour','NUMdayssincesolstice']].values
# useatnmoonphaseCRABy = useatnmoonphaseCRABNOGAPS['atn']


# from pygam import LinearGAM
# gam = LinearGAM().gridsearch(useatnmoonphaseCRABx, useatnmoonphaseCRABy)

# plt.plot(gam.predict(useatnmoonphaseCRABx))
# useatnmoonphaseCRABNOGAPS['atn'].plot(figsize =(12,8))

# #pdeps, confi = gam.partial_dependence(width = 0.95)
# #this doesn't work. try resetting index

# #https://pygam.readthedocs.io/en/latest/notebooks/quick_start.html
# for i, term in enumerate(gam.terms):
#     if term.isintercept:
#         continue

#     XX = gam.generate_X_grid(term=i)
#     pdep, confi = gam.partial_dependence(term=i, X=XX, width=0.95)

#     plt.figure()
#     plt.plot(XX[:, term.feature], pdep)
#     plt.plot(XX[:, term.feature], confi, c='r', ls='--')
#     plt.title(repr(term))
#     plt.show()

# gam.summary()

# # LinearGAM                                                                                                 
# # =============================================== ==========================================================
# # Distribution:                        NormalDist Effective DoF:                                     25.4885
# # Link Function:                     IdentityLink Log Likelihood:                                -13702.2544
# # Number of Samples:                         1693 AIC:                                            27457.4858
# #                                                 AICc:                                           27458.3602
# #                                                 GCV:                                                0.0496
# #                                                 Scale:                                              0.0483
# #                                                 Pseudo R-Squared:                                   0.1972
# # ==========================================================================================================
# # Feature Function                  Lambda               Rank         EDoF         P > x        Sig. Code   
# # ================================= ==================== ============ ============ ============ ============
# # s(0)                              [15.8489]            20           9.3          1.11e-16     ***         
# # s(1)                              [15.8489]            20           8.4          2.03e-14     ***         
# # s(2)                              [15.8489]            20           7.7          1.11e-16     ***         
# # intercept                                              1            0.0          1.11e-16     ***         
# # ==========================================================================================================
# # Significance codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1







# list(max_scores_night.columns)


# max_scores_nightx = max_scores_night[['mooninterp', 'hour','NUMdayssincesolstice']].values
# max_scores_nighty = max_scores_night['atn']








# from pygam import LinearGAM
# gam = LinearGAM().gridsearch(max_scores_nightx, max_scores_nighty)

# plt.plot(gam.predict(max_scores_nightx))
# useatnmoonphaseCRABNOGAPS['atn'].plot(figsize =(12,8))

# #pdeps, confi = gam.partial_dependence(width = 0.95)
# #this doesn't work. try resetting index
# for i, term in enumerate(gam.terms):
#     if term.isintercept:
#         continue

#     XX = gam.generate_X_grid(term=i)
#     pdep, confi = gam.partial_dependence(term=i, X=XX, width=0.95)

#     plt.figure()
#     plt.plot(XX[:, term.feature], pdep)
#     plt.plot(XX[:, term.feature], confi, c='r', ls='--')
#     plt.title(repr(term))
#     plt.show()


# gam.summary()
# # gam.summary()
# # LinearGAM                                                                                                 
# # =============================================== ==========================================================
# # Distribution:                        NormalDist Effective DoF:                                      5.6974
# # Link Function:                     IdentityLink Log Likelihood:                                  -348.6154
# # Number of Samples:                           75 AIC:                                              710.6255
# #                                                 AICc:                                             712.1574
# #                                                 GCV:                                                0.0842
# #                                                 Scale:                                              0.0728
# #                                                 Pseudo R-Squared:                                   0.2007
# # ==========================================================================================================
# # Feature Function                  Lambda               Rank         EDoF         P > x        Sig. Code   
# # ================================= ==================== ============ ============ ============ ============
# # s(0)                              [1000.]              20           3.1          1.25e-02     *           
# # s(1)                              [1000.]              20           1.7          5.90e-01                 
# # s(2)                              [1000.]              20           0.9          1.75e-01                 
# # intercept                                              1            0.0          1.11e-16     ***         
# # ==========================================================================================================
# # Significance codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# # WARNING: Fitting splines and a linear function to a feature introduces a model identifiability problem
# #          which can cause p-values to appear significant when they are not.

# # WARNING: p-values calculated in this manner behave correctly for un-penalized models or models with
# #          known smoothing parameters, but when smoothing parameters have been estimated, the p-values
# #          are typically lower than they should be, meaning that the tests reject the null too readily.
# # C:\Users\chuffard\AppData\Local\Temp\ipykernel_16868\3358381670.py:1: UserWarning: KNOWN BUG: p-values computed in this summary are likely much smaller than they should be. 
 
# # Please do not make inferences based on these values! 

# # Collaborate on a solution, and stay up to date at: 
# # github.com/dswah/pyGAM/issues/163 

# #   gam.summary()


# max_scores_morningx = max_scores_morning[['mooninterp', 'hour','NUMdayssincesolstice']].values
# max_scores_morningy = max_scores_morning['atn']


# from pygam import LinearGAM
# gam = LinearGAM().fit(max_scores_morningx, max_scores_morningy)

# plt.plot(gam.predict(max_scores_morningx))
# useatnmoonphaseCRABNOGAPS['atn'].plot(figsize =(12,8))

# #pdeps, confi = gam.partial_dependence(width = 0.95)
# #this doesn't work. try resetting index
# for i, term in enumerate(gam.terms):
#     if term.isintercept:
#         continue

#     XX = gam.generate_X_grid(term=i)
#     pdep, confi = gam.partial_dependence(term=i, X=XX, width=0.95)

#     plt.figure()
#     plt.plot(XX[:, term.feature], pdep)
#     plt.plot(XX[:, term.feature], confi, c='r', ls='--')
#     plt.title(repr(term))
#     plt.show()


# gam.summary()
# # LinearGAM                                                                                                 
# # =============================================== ==========================================================
# # Distribution:                        NormalDist Effective DoF:                                     24.6611
# # Link Function:                     IdentityLink Log Likelihood:                                  -628.6854
# # Number of Samples:                           72 AIC:                                             1308.6929
# #                                                 AICc:                                            1338.8723
# #                                                 GCV:                                                0.0703
# #                                                 Scale:                                               0.029
# #                                                 Pseudo R-Squared:                                   0.6076
# # ==========================================================================================================
# # Feature Function                  Lambda               Rank         EDoF         P > x        Sig. Code   
# # ================================= ==================== ============ ============ ============ ============
# # s(0)                              [0.6]                20           11.1         7.41e-02     .           
# # s(1)                              [0.6]                20           8.6          1.61e-01                 
# # s(2)                              [0.6]                20           5.0          2.56e-03     **          
# # intercept                                              1            0.0          1.11e-16     ***         
# # ==========================================================================================================
# # Significance codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# # WARNING: Fitting splines and a linear function to a feature introduces a model identifiability problem
# #          which can cause p-values to appear significant when they are not.

# # WARNING: p-values calculated in this manner behave correctly for un-penalized models or models with
# #          known smoothing parameters, but when smoothing parameters have been estimated, the p-values
# #          are typically lower than they should be, meaning that the tests reject the null too readily.
# # C:\Users\chuffard\AppData\Local\Temp\ipykernel_16868\3358381670.py:1: UserWarning: KNOWN BUG: p-values computed in this summary are likely much smaller than they should be. 

# #plot morning atn in relation to moon phase
# bins = [0, .25, .5, .75, 1]
# morningSESMARS['moonbins'] = pd.cut(morningSESMARS['mooninterp'], bins)

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



# #plot morning atn in relation to moon phase
# bins = [0, .25, .5, .75, 1]
# nightSESMARS['moonbins'] = pd.cut(nightSESMARS['mooninterp'], bins)

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



# mako4 = sns.color_palette("rocket", 4)
# sns.lineplot(data=morningSESMARS, 
#              x="mooninterp", 
#              y="atn", 
#              hue = 'moonbins',
#              palette = mako4)
# plt.legend(loc=(1.04, 0))
# plt.savefig('morning atn vs moon phaseI.pdf')


# sns.scatterplot(data=max_scores_night, x='mooninterp', y='atn', hue='NUMdayssincesolstice')
# plt.legend(loc=(1.04, 0))



# import numpy as np
# import seaborn as sns
# import matplotlib.pyplot as plt
# labels=np.array(['Midnight', '1', '2', '3', '4', '5', '6','7','8','9','10','11', '12', '13', '14', '15', '16','17','18','19','20','21', '22', '23'])

# #wide data
# WIDEuseatnmoonphasewideHour =  useatnmoonphaseCRAB.dropna().pivot(index='datetime', columns='hour', values='atn')

# # # number of variable
# # categories=list(df)[1:]
# # N = len(categories)
# # # What will be the angle of each axis in the plot? (we divide the plot / number of variable)
# # angles = [n / float(N) * 2 * pi for n in range(N)]
# # angles += angles[:1]
# # # Initialise the spider plot
# # ax = plt.subplot(111, polar=True)
# # # If you want the first axis to be on top:
# # ax.set_theta_offset(pi / 2)
# # ax.set_theta_direction(-1)
# # # Draw one axe per variable + add labels labels yet
# # plt.xticks(angles[:-1], categories)
# # # Draw ylabels
# # ax.set_rlabel_position(0)
# # plt.yticks([1,2,3,4], ["1","2","3","4"], color="grey", size=7)
# # plt.ylim(0,5)
# # # ------- PART 2: Add plots
# # # Name1
# # values=df.loc[0].drop('name').values.flatten().tolist()
# # values += values[:1]
# # ax.plot(angles, values, linewidth=1, linestyle='solid', label="John")
# # ax.fill(angles, values, 'b', alpha=0.1)
# # # Name2
# # values=df.loc[1].drop('name').values.flatten().tolist()
# # values += values[:1]
# # ax.plot(angles, values, linewidth=1, linestyle='solid', label="David")
# # ax.fill(angles, values, 'r', alpha=0.1)
# # # Add legend
# # plt.legend(loc='upper right', bbox_to_anchor=(0.1, 0.1))



#read csv from folder and combine into one csv file then save
#bring in fluoresence
# read in data; use first line as header

import pandas as pd
import os

# from statsmodels.tsa.tsatools import detrend


# Main_Hourly_SES_MARSComplete['detrended690Fluoro'] = detrend(Main_Hourly_SES_MARSComplete['Fluoro690'],  order=2)


# replace with your folder's path
folder_path = r'G:/Data/compiling'

all_files = os.listdir(folder_path)

# Filter out non-CSV files
csv_files = [f for f in all_files if f.endswith('.csv')]

# Create a list to hold the dataframes
df_list = []

for csv in csv_files:
    file_path = os.path.join(folder_path, csv)
    try:
        # Try reading the file using default UTF-8 encoding
        df = pd.read_csv(file_path)
        df_list.append(df)
    except UnicodeDecodeError:
        try:
            # If UTF-8 fails, try reading the file using UTF-16 encoding with tab separator
            df = pd.read_csv(file_path, sep='\t', encoding='utf-16')
            df_list.append(df)
        except Exception as e:
            print(f"Could not read file {csv} because of error: {e}")
    except Exception as e:
        print(f"Could not read file {csv} because of error: {e}")

# Concatenate all data into one DataFrame
big_df = pd.concat(df_list, ignore_index=True)

# Save the final result to a new CSV file
big_df.to_csv(os.path.join(folder_path, 'combined_fluoroMARS.csv'), index=False)



Main_SES_Fluoro = pd.read_csv('G:/Data/compiling/combined_fluoroMARS.csv', parse_dates=True, encoding='unicode_escape')

list(Main_SES_Fluoro.columns)
# ['SES Timestamp (ISO601)',
#  'Dark Ref Integration Time',
#  'Dark Ref 690 Counts',
#  'Dark Ref 590 Counts',
#  'Measure Ref Integration Time',
#  'Measure Ref 690 Counts',
#  'Measure Ref 590 Counts',
#  'Luminescense Integration Time',
#  'Luminescense 690 Counts',
#  'Luminescense 590 Counts',
#  'Fluorescense Integration Time',
#  'Fluorescense 690 Counts',
#  'Fluorescence690DATA',
#  'Fluorescense 590 Counts',
#  'Fluorescence590DATA',
#  'Phosphorescense Integration Time',
#  'Phosphorescense 690 Counts',
#  'Phosphorescense 590 Counts',
#  'Constant Target Fluorescense Integration Time',
#  'Constant Target Fluorescense 690 Counts',
#  'Constant Target Fluorescense 590 Counts',
#  'Unnamed: 21']

Main_SES_Fluoro['Fluoro690']= Main_SES_Fluoro['Fluorescense 690 Counts']-Main_SES_Fluoro['Measure Ref 690 Counts']
Main_SES_Fluoro['Fluoro590']= Main_SES_Fluoro['Fluorescense 590 Counts']-Main_SES_Fluoro['Measure Ref 590 Counts']
Main_SES_Fluoro =  Main_SES_Fluoro[['SES Timestamp (ISO601)', 'Fluoro690','Fluoro590']] 
Main_SES_Fluoro['datetime'] = pd.to_datetime(Main_SES_Fluoro['SES Timestamp (ISO601)'])
#Main_SES_Fluoro = Main_SES_Fluoro.loc[:'2023-08-14']

a = Main_SES_Fluoro[Main_SES_Fluoro['datetime'] <= datetime(2023, 8, 15)]

#get rid of two weird rows in the 1970's
Main_SES_Fluoro = a[a['datetime']>= datetime(2022, 1, 1)]
# Main_SES_Fluoro['hour'] = Main_SES_Fluoro['datetime'].dt.hour
# Main_SES_Fluoro['year'] = Main_SES_Fluoro['datetime'].dt.year
# Main_SES_Fluoro['month'] = Main_SES_Fluoro['datetime'].dt.month
# Main_SES_Fluoro['day'] = Main_SES_Fluoro['datetime'].dt.day
# Main_SES_Fluoro['Year-Week'] = Main_SES_Fluoro['datetime'].dt.strftime('%Y-%U')
# Main_SES_Fluoro['justdate'] = Main_SES_Fluoro['datetime'].dt.date

#find hourly average
Main_SES_FluoroHourly = Main_SES_Fluoro.resample('1H', on='datetime').mean()
  
MainSES_fluoro_atn_moon_sun = Main_Hourly_SES_MARS.merge(Main_SES_FluoroHourly, how='outer',on='datetime')

MainSES_fluoro_atn_moon_sun.to_csv('MainSES_fluoro_atn_moon_sun.csv')


#Import current speed
current_speed = pd.read_csv("G:/Deployments/MARS-230508-01/TIlt_meters/No_2_2212301/currents/2212301_MARS2023_tiltMeter_DATA.csv", parse_dates=True,index_col="Date_Time")
current_speed.dtypes
current_speed= current_speed.loc[current_speed.notna().all(axis='columns')]
current_speed['datetime'] = pd.to_datetime(current_speed['ISO 8601 Time'])
trimhcurrent_speed = current_speed.loc[:'2023-08-14']  

trimhcurrent_speed["Speed_cm_s"].plot()

#find hourly average of SES time series
hourly_current = trimhcurrent_speed.resample('1H').mean()
hourly_current = hourly_current.to_frame()

hourly_current['datetime'] = hourly_current.index

hourly_current['datetime'] = pd.to_datetime(hourly_current['datetime'])
#now merge with main dataset


MainSES_fluoro_atn_moon_sun_current = Main_Hourly_SES_MARS.merge(hourly_current, how='outer',on='datetime')

MainSES_fluoro_atn_moon_sun_current.to_csv("MainSES_fluoro_atn_moon_sun_current.csv")
#try spectray density welch
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html
fs = 1
f, Pxx_den = signal.welch(hourly_current['Speed'])
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 current speed acf hourly.pdf')