# -*- coding: utf-8 -*-
"""
Created on Thu Sep 19 16:07:51 2024

@author: chuffard
"""

#SES 2024 MARS data
# -*- coding: utf-8 -*-
"""
Created on Mon May 13 13:39:04 2024

@author: chuffard
"""
# libraries
import seaborn as sns
import matplotlib.pyplot as pl
import math
import pandas as pd
from shapely.geometry import Polygon
import datetime
import seaborn as sns
import statistics
import matplotlib.pyplot as plt
import os
import re
from scipy.cluster.hierarchy import dendrogram, linkage
import plotly.figure_factory as ff
import pandas as pd
import numpy as np
import scipy.cluster.hierarchy
from scipy.cluster.hierarchy import ward, dendrogram, leaves_list
from scipy.spatial.distance import pdist
import numpy.ma as ma
from scipy import stats
from scipy.stats import zscore
from numpy.lib.stride_tricks import as_strided
import scipy.stats
from numpy.lib import pad
import tanglegram as tg
import warnings
from matplotlib.dates import DateFormatter
import matplotlib.dates as mdates

def letter_annotation(ax, xoffset, yoffset, letter):
 ax.text(xoffset, yoffset, letter, transform=ax.transAxes,
         size=12, weight='bold')
 




#import attenuance

SES_RAW_sorted=pd.read_csv("G:/Analyses/2024 MARS deployment/SES_RAW.csv", parse_dates=True)
SES_RAW_sorted["date_time"] = pd.to_datetime(SES_RAW_sorted["datetime"])

# SES_RAW_sorted = pd.read_csv("G:/Analyses/2024 MARS deployment/SES_RAW_sorted.csv", parse_dates=True)

# SES_RAW_sorted["date_time"] = pd.to_datetime(SES_RAW_sorted["datetime"])
SES_RAW_sorted = SES_RAW_sorted.drop('datetime', axis=1)


#drop inf values
SES_RAW_sorted = SES_RAW_sorted.replace([np.inf, -np.inf], np.nan).dropna(subset=['atn'])
#filter out periods with 0 collect time
SES_RAW_sorted =SES_RAW_sorted[SES_RAW_sorted['collect_time'] == '5min']
SES_RAW_sorted['rm150'] = SES_RAW_sorted['atn'].rolling(150, min_periods=23).mean()

McLaneData = pd.read_excel('G:/Deployments/MARS-2024/Trigger trap/Carbon_data/Trigger_trap_filter_log.xlsx', sheet_name='Cups')
McLaneData['date_time'] = pd.to_datetime(McLaneData['datetime'])
McLaneData = McLaneData.drop('datetime', axis=1)

 
TriggerData = pd.read_excel('G:/Deployments/MARS-2024/Trigger trap/Carbon_data/Trigger_trap_filter_log.xlsx', sheet_name='TriggerDates')
TriggerData['date_time'] = pd.to_datetime(TriggerData['datetime'])
TriggerData = TriggerData.drop('datetime', axis=1)
TriggerData['hour'] = TriggerData['date_time'].dt.hour

TriggerData['date_time'] = TriggerData['date_time'].dt.floor('H')

SES_RAW_sorted= SES_RAW_sorted.merge(McLaneData, 'outer')

HourlySES_RAW_sorted=  SES_RAW_sorted.groupby([pd.Grouper(key='date_time', freq='H')]).mean() 
HourlySES_RAW_sorted = pd.merge(HourlySES_RAW_sorted, TriggerData, 'outer', on='date_time')
HourlySES_RAW_sorted.rename(columns={'hour_x': 'hour'}, inplace=True)
HourlySES_RAW_sorted.drop('hour_y', axis=1, inplace=True)
DailySES_RAW_sorted=  SES_RAW_sorted.groupby([pd.Grouper(key='date_time', freq='D')]).mean() 
WeeklySES_RAW_sorted=  SES_RAW_sorted.groupby([pd.Grouper(key='date_time', freq='W')]).mean() 
#lots of missing hours. Need to get to the bottom of that

#Carbon data






#Current meter path

Current_meter = pd.read_csv("G:/Deployments/MARS-2024/Current meters/SES current meter 2024/2212301_MARS2024SES_(0)_Current.csv", parse_dates=True)

Current_meter.dtypes
# ISO 8601 Time         object
# Speed (cm/s)         float64
# Heading (degrees)    float64
# Velocity-N (cm/s)    float64
# Velocity-E (cm/s)    float64
Current_meter["date_time"] = pd.to_datetime(Current_meter["ISO 8601 Time"])
Current_meter = Current_meter[(Current_meter['date_time'] > '2024-06-21') & (Current_meter['date_time'] < '2024-09-17')]

#decimate
Decimated_Current_meter = Current_meter.iloc[::10, :]
HourlyCurrents=  Current_meter.groupby([pd.Grouper(key='date_time', freq='H')]).mean() 
DailyCurrents=  Current_meter.groupby([pd.Grouper(key='date_time', freq='D')]).mean() 
WeeklyCurrents=  Current_meter.groupby([pd.Grouper(key='date_time', freq='W')]).mean() 



#pHsensor = pd.read_csv("G:/Deployments/MARS-2024/pH sensors/SES pH sensor/2024-06-19_114930_SES2024.txt" ,sep="\t", header=23, encoding='utf-8')

pHTemp = pd.read_csv("G:/Deployments/MARS-2024/pH sensors/SES pH sensor/SES2024csv.csv")
pHTemp.dtypes
pHTemp["date_time"] = pd.to_datetime(pHTemp["date_time"])
# pHTemp = pHsensor.iloc[:, [0,1,3,13]]
# pHTemp.rename(columns={ pHTemp.columns[0]: "date" }, inplace = True)
# pHTemp.rename(columns={ pHTemp.columns[1]: "time" }, inplace = True)
# pHTemp.rename(columns={ pHTemp.columns[2]: "pH" }, inplace = True)

#could not convert string to float: '<6.5'
#pHTemp =   pHTemp[pHTemp["pH"].str.contains("<") == False]
#pHTemp['pH'] = pHTemp['pH'].astype(float, errors = 'raise')


#pHTemp.rename(columns={ pHTemp.columns[3]: "TempC" }, inplace = True)
#pHTemp['date_time'] = pd.to_datetime(pHTemp['date'] + ' ' + pHTemp['time'])
pHTemp = pHTemp[(pHTemp['date_time'] > '2024-06-21') & (pHTemp['date_time'] < '2024-09-17')]

HourlypHTemp=  pHTemp.groupby([pd.Grouper(key='date_time', freq='H')]).mean() 
#HourlypHTemp['date_time'] = pd.to_datetime(HourlypHTemp['date_time'])
DailypHTemp=  pHTemp.groupby([pd.Grouper(key='date_time', freq='D')]).mean() 
#DailypHTemp['date_time'] = pd.to_datetime(DailypHTemp['date_time'])
WeeklypHTemp=  pHTemp.groupby([pd.Grouper(key='date_time', freq='W')]).mean() 
#WeeklypHTemp['date_time'] = pd.to_datetime(WeeklypHTemp['date_time'])


##fluorometer
MARSESfluoro2024 = pd.read_csv('G:/Deployments\MARS-2024/dep-2024.06.20-02-SED-AT-MARS-2024/fluorometer.csv')
MARSESfluoro2024["date_time"] = pd.to_datetime(MARSESfluoro2024["SES Timestamp (ISO601)"])
MARSESfluoro2024 = MARSESfluoro2024[(MARSESfluoro2024['date_time'] > '2024-06-21') & (MARSESfluoro2024['date_time'] < '2024-09-17')]
MARSESfluoro2024['Fluoro690']= MARSESfluoro2024['Fluorescense 690 Counts']-MARSESfluoro2024['Measure Ref 690 Counts']
MARSESfluoro2024['Fluoro590']= MARSESfluoro2024['Fluorescense 590 Counts']-MARSESfluoro2024['Measure Ref 590 Counts']
MARSESfluoro2024.dtypes
HourlyMARSESfluoro2024=  MARSESfluoro2024.groupby([pd.Grouper(key='date_time', freq='H')]).mean() 
DailyMARSESfluoro2024=  MARSESfluoro2024.groupby([pd.Grouper(key='date_time', freq='D')]).mean() 
WeeklyMARSESfluoro2024=  MARSESfluoro2024.groupby([pd.Grouper(key='date_time', freq='W')]).mean() 


#Combine 
dflist = [HourlySES_RAW_sorted, DailySES_RAW_sorted, WeeklySES_RAW_sorted, HourlyCurrents, DailyCurrents,WeeklyCurrents,
          HourlypHTemp,DailypHTemp, WeeklypHTemp, HourlyMARSESfluoro2024,DailyMARSESfluoro2024, WeeklyMARSESfluoro2024]

for df in dflist:
    df.reset_index(drop=False, inplace=True)
    df['date_time'] = pd.to_datetime(df['date_time'])



HourlyMAIN= HourlySES_RAW_sorted.merge(HourlyCurrents, 'outer')
HourlyMAIN= HourlyMAIN.merge(HourlypHTemp, 'outer')
HourlyMAIN= HourlyMAIN.merge(MARSESfluoro2024, 'outer')

HourlyMAIN = HourlyMAIN[(HourlyMAIN['date_time'] > '2024-06-21') & (HourlyMAIN['date_time'] < '2024-09-17')]
HourlyMAIN = HourlyMAIN.groupby([pd.Grouper(key='date_time', freq='H')]).mean() 

#sort by datetime
HourlyMAIN = HourlyMAIN.sort_values(by='date_time')

#find the rolling mean= 150 samples, minimum sample size 23.
#HourlyMAIN['rm150'] = HourlyMAIN['atn'].rolling(150, min_periods=23).mean()

HourlyMAIN = HourlyMAIN.replace([np.inf, -np.inf], np.nan).dropna(subset=['atn'])

DailyMAIN= DailySES_RAW_sorted.merge(DailyCurrents, 'outer')
DailyMAIN= DailyMAIN.merge(DailypHTemp, 'outer')
DailyMAIN= DailyMAIN.merge(DailyMARSESfluoro2024, 'outer')

DailyMAIN = DailyMAIN[(DailyMAIN['date_time'] > '2024-06-21') & (DailyMAIN['date_time'] < '2024-09-17')]

WeeklyMAIN= WeeklySES_RAW_sorted.merge(WeeklyCurrents, 'outer')
WeeklyMAIN= WeeklyMAIN.merge(WeeklypHTemp, 'outer')

WeeklyMAIN = WeeklyMAIN[(WeeklyMAIN['date_time'] > '2024-06-21') & (WeeklyMAIN['date_time'] < '2024-09-17')]


HighTriggers = HourlySES_RAW_sorted[HourlySES_RAW_sorted['type'] == 'High']
HighTriggerDates = HighTriggers['date_time'].tolist()
HighTriggerHours = HighTriggers['hour'].tolist()

LowTriggers = HourlySES_RAW_sorted[HourlySES_RAW_sorted['type'] == 'Low']
LowTriggerDates = LowTriggers['date_time'].tolist()
LowTriggerHours = LowTriggers['hour'].tolist()

SchedTriggers = HourlySES_RAW_sorted[HourlySES_RAW_sorted['type'] == 'Scheduled']
SchedTriggerDates = SchedTriggers['date_time'].tolist()
SchedTriggeHourss = SchedTriggers['hour'].tolist()


#Melted hours
HourlyMAIN.dtypes
HourlyMAIN['value'] = HourlyMAIN['atn']


max_atn = HourlyMAIN['atn'].max()
print(max_atn)

max_mass_flux = HourlyMAIN['MassFluxmgm2d1'].max()
print(max_mass_flux)
correctionfactor = max_atn/max_mass_flux

HourlyMAIN['coplotMassFlux'] = HourlyMAIN['MassFluxmgm2d1']*correctionfactor

HoursMelt = pd.melt(HourlyMAIN, id_vars=['hour'], value_vars=['atn'])
HoursMeltCurrents = pd.melt(HourlyMAIN, id_vars=['hour'], value_vars=['Speed (cm/s)'])

#For plot by sediment trap cup number

HourlyMAIN.columns

HourlyMAIN['atnInterp'] = HourlyMAIN['atn'].interpolate(method='linear')

massdata = HourlyMAIN[['Sediment trap cup#', 'MassFluxmgm2d1', 'atnInterp']]

massdata = massdata.dropna()

cupatn_means = massdata.groupby('Sediment trap cup#')['atnInterp'].mean()

cupatn_means= cupatn_means.reset_index()
#rename column
cupatn_means.columns.values[1] = 'atninterpMEAN'
massdata= massdata.merge(cupatn_means, 'outer')

massdata['step'] = massdata['atninterpMEAN']**0.87
massdata['modeledmass'] = massdata['step']*10300

#_____________________________________PLOTS__________________________________________

#____________________________________TIME SERIES PANEL


# instanciate figure and axes instances
fig,axes = plt.subplots(5,1, figsize=(8,11), dpi=300)
params = {'mathtext.default': 'regular' } 

# first plot

sns.scatterplot(data=SES_RAW_sorted, x="date_time", y = 'atn',color = 'black',ax=axes[0], alpha = 0.05)
sns.lineplot(data=SES_RAW_sorted, x="date_time", y = 'rm150',color = 'red',ax=axes[0]) 
sns.lineplot(data=DailyMAIN, x="date_time", y = 'atn', color = 'blue',ax=axes[0]) 
#sns.scatterplot(data=Triggers, x="date_time", y = 'atn', color = 'red',ax=axes[0],s=50)  #first position
sns.scatterplot(data=HighTriggers, x="date_time", y = 'atn', color = 'red',ax=axes[0],s=70,marker='s') 
sns.scatterplot(data=LowTriggers, x="date_time", y = 'atn', color = 'blue',ax=axes[0],s=70,marker='s') 
sns.scatterplot(data=SchedTriggers, x="date_time", y = 'atn', color = 'black',ax=axes[0],s=70)  #first position
sns.scatterplot(data=HourlyMAIN, x="date_time", y = 'coplotMassFlux', color = 'cyan',ax=axes[0],s=40, marker="^") 
axes[0].set_xlabel("")
axes[0].set_ylabel("Attenuance \nmass flux (g m$^{-2}$ d$^{-1}$) \nmax = 9.6 ")
axes[0].xaxis.set_major_locator(mdates.MonthLocator())
#axes[0].xaxis.set_major_formatter(xfmt)
letter_annotation(axes[0], -.1, 1, 'A')

#second plot

sns.scatterplot(data=HourlyMAIN, x="date_time", y = 'Fluoro690',color = 'black',ax=axes[1], alpha = 0.05)
sns.lineplot(data=DailyMAIN, x="date_time", y = 'Fluoro690',color = 'blue',ax=axes[1], alpha = 0.05)
sns.lineplot(data=DailyMAIN, x="date_time", y = 'Fluoro690', color = 'blue',ax=axes[1]) 
axes[1].set_xlabel("")
axes[1].set_ylabel('Fluorescence 690')
axes[1].xaxis.set_major_locator(mdates.MonthLocator())
letter_annotation(axes[1], -.1, 1, 'B')


sns.scatterplot(data=HourlyMAIN, x="date_time", y = 'TempC',color = 'black',ax=axes[2], alpha = 0.05) 
sns.lineplot(data=DailyMAIN, x="date_time", y = 'TempC', color = 'blue',ax=axes[2]) 
axes[2].set_xlabel("")
axes[2].set_ylabel("TempC")
axes[2].xaxis.set_major_locator(mdates.MonthLocator())
letter_annotation(axes[2], -.1, 1, 'C')

#second plot

sns.scatterplot(data=HourlyMAIN, x="date_time", y = 'Speed (cm/s)',color = 'black',ax=axes[3], alpha = 0.05)
sns.lineplot(data=DailyMAIN, x="date_time", y = 'Speed (cm/s)', color = 'blue',ax=axes[3]) 
axes[3].set_xlabel("")
axes[3].set_ylabel("Current speed \n(cm/s)")
axes[3].xaxis.set_major_locator(mdates.MonthLocator())
letter_annotation(axes[3], -.1, 1, 'D')



#fifth plot

sns.scatterplot(data=HourlyMAIN, x="date_time", y = 'pH',color = 'black',ax=axes[4], alpha = 0.05)
sns.lineplot(data=DailyMAIN, x="date_time", y = 'pH', color = 'blue',ax=axes[4]) 
axes[4].set_xlabel("")
axes[4].set_ylabel("pH")
axes[4].xaxis.set_major_locator(mdates.MonthLocator())
letter_annotation(axes[4], -.1, 1, 'E')


plt.tight_layout()

plt.savefig('G:/Analyses/2024 MARS deployment/SES_TS.pdf')





fig,axes = plt.subplots(2,1, figsize=(8,4), dpi=300)
date_form = DateFormatter("%b")
    

# first plot 
sns.lineplot(data=HoursMelt, x="hour", y = 'value',ax=axes[0])
sns.lineplot(data=HourlyMAIN, x="hour", y="atn", hue="month",ax=axes[0], ci=None)
axes[0].set_xlabel("hour")
axes[0].set_ylabel("Attenuance")
letter_annotation(axes[0], -.25, 1, 'A')
axes[0].legend(bbox_to_anchor=(1.05, 1), loc='upper left',title='Month')



#second plot
sns.lineplot(data=HoursMeltCurrents, x="hour", y = 'value',ax=axes[1])
sns.lineplot(data=HourlyMAIN, x="hour", y="Speed (cm/s)", hue="month",ax=axes[1], ci=None)
axes[1].set_xlabel("hour")
axes[1].set_ylabel("Current speed \n(cm/s)")
letter_annotation(axes[1], -.25, 1, 'B')
axes[1].legend(bbox_to_anchor=(1.05, 1), loc='upper left',title='Month')

for ax in axes:
    ax.yaxis.set_label_coords(-.15, 0.5)

plt.tight_layout()


plt.savefig('G:/Analyses/2024 MARS deployment/hourly.pdf')


#Direct mass flux and atn values per cup

# fig,axes = plt.subplots(2,1, figsize=(8,6), dpi=300)


# # first plot 
# sns.lineplot(data = massdata, x = "modeledmass", y = "atninterpMEAN",ax=axes[0] ,color='black',alpha = 0.3,  linestyle='--')
# sns.scatterplot(data=massdata, x='MassFluxmgm2d1', y='atnInterp', hue='Sediment trap cup#', palette='Set2',ax=axes[0])
# sns.scatterplot(data=massdata, x='MassFluxmgm2d1', y='atninterpMEAN', hue='Sediment trap cup#',  s=50, edgecolor='black',palette='Set2',ax=axes[0])
# sns.scatterplot(data = massdata, x = "modeledmass", y = "atninterpMEAN",ax=axes[0] ,marker='D', alpha = 0.3,  edgecolor='black',hue='Sediment trap cup#', palette='Set2')
# sns.regplot(data=massdata, x="MassFluxmgm2d1", y = 'atninterpMEAN',ax=axes[0], ci=None, marker='')

# handles, labels = axes[0].get_legend_handles_labels()

# unique_labels = {}
# unique_handles = []
# for handle, label in zip(handles, labels):
#     if label not in unique_labels:
#         unique_labels[label] = True
#         unique_handles.append(handle)

# # Update the legend
# axes[0].legend(unique_handles, unique_labels.keys(),bbox_to_anchor=(1.05, 1), loc='upper left',title='Sediment trap cup #')
# axes[0].set_xlabel("Mass flux (mg m$^{-2}$ d$^{-1}$)")
# axes[0].set_ylabel("Attenuance")
# #axes[0].set_xlim(0, 12000)

# letter_annotation(axes[0], -.25, 1, 'A')
# #axes[0].legend(bbox_to_anchor=(1.05, 1), loc='upper left',title='Sediment trap cup #')
# #axes[0].set_xscale('log')


# # #second plot
# # sns.regplot(data=massdata, x="MassFluxmgm2d1", y = 'atninterpMEAN',ax=axes[1])
# # #sns.lineplot(data = massdata, x = "modeledmass", y = "atninterpMEAN",ax=axes[1] ,color='black', linestyle='--')
# # # sns.scatterplot(data=massdata, x='MassFluxmgm2d1', y='atninterpMEAN', hue='Sediment trap cup#', marker='D',  edgecolor='black',palette='Set2',ax=axes[1])
# # # sns.lineplot(data = massdata, x = "modeledmass", y = "atninterpMEAN",ax=axes[1] ,color='black', linestyle='--')
# # # axes[1].set_xlabel("Mass flux (mg m$^{-2}$ d$^{-1}$)")
# # axes[1].set_ylabel("Attenuance")
# # letter_annotation(axes[1], -.25, 1, 'B')
# # #axes[1].set_xlim(0, 12000)
# # #axes[1].set_xscale('log')

# #Direct mass flux and atn values per cup

#fig,axes = plt.subplots(2,1, figsize=(8,6), dpi=300)
plt.figure(figsize=(8, 3), dpi=300) 

sns.lineplot(data = massdata, x = "modeledmass", y = "atninterpMEAN",color='black',alpha = 0.5,  linestyle='--')
sns.scatterplot(data=massdata, x='MassFluxmgm2d1', y='atnInterp', hue='Sediment trap cup#', palette='Set2', label='_no_legend')
sns.scatterplot(data=massdata, x='MassFluxmgm2d1', y='atninterpMEAN', hue='Sediment trap cup#',  s=50, edgecolor='black',palette='Set2', label='_no_legend')
sns.scatterplot(data = massdata, x = "modeledmass", y = "atninterpMEAN",marker='D', alpha = 0.5,  edgecolor='black',hue='Sediment trap cup#', palette='Set2', label='_no_legend')
sns.regplot(data=massdata, x="MassFluxmgm2d1", y = 'atninterpMEAN',ci=None, marker='')
plt.xlabel("Mass flux (mg m$^{-2}$ d$^{-1}$)")
plt.ylabel("Attenuance")

# Get handles and labels
handles, labels = plt.gca().get_legend_handles_labels()

# Remove duplicate labels
new_labels, new_handles = [], []
for handle, label in zip(handles, labels):
    if label not in new_labels:
        new_labels.append(label)
        new_handles.append(handle)
# Set new legend
plt.legend(new_handles, new_labels,bbox_to_anchor=(1.05, 1), loc='upper left',title='Sediment trap cup #')

plt.tight_layout()

plt.savefig('G:/Analyses/2024 MARS deployment/Mass vs attenuance.pdf')






fig,axes = plt.subplots(2,1, figsize=(8,6), dpi=300)


# first plot 
sns.lineplot(data = massdata, x = "modeledmass", y = "atninterpMEAN",ax=axes[0] ,color='black',alpha = 0.3,  linestyle='--')
sns.scatterplot(data=massdata, x='MassFluxmgm2d1', y='atnInterp', hue='Sediment trap cup#', palette='Set2',ax=axes[0])
sns.scatterplot(data=massdata, x='MassFluxmgm2d1', y='atninterpMEAN', hue='Sediment trap cup#',  s=50, edgecolor='black',palette='Set2',ax=axes[0])
sns.scatterplot(data = massdata, x = "modeledmass", y = "atninterpMEAN",ax=axes[0] ,marker='D', alpha = 0.3,  edgecolor='black',hue='Sediment trap cup#', palette='Set2')
sns.regplot(data=massdata, x="MassFluxmgm2d1", y = 'atninterpMEAN',ax=axes[0], ci=None, marker='')

handles, labels = axes[0].get_legend_handles_labels()

unique_labels = {}
unique_handles = []
for handle, label in zip(handles, labels):
    if label not in unique_labels:
        unique_labels[label] = True
        unique_handles.append(handle)

# Update the legend
axes[0].legend(unique_handles, unique_labels.keys(),bbox_to_anchor=(1.05, 1), loc='upper left',title='Sediment trap cup #')
axes[0].set_xlabel("Mass flux (mg m$^{-2}$ d$^{-1}$)")
axes[0].set_ylabel("Attenuance")
#axes[0].set_xlim(0, 12000)

letter_annotation(axes[0], -.25, 1, 'A')

plt.tight_layout()

plt.savefig('G:/Analyses/2024 MARS deployment/Mass vs attenuance.pdf')

#STATS
res = stats.spearmanr(massdata['MassFluxmgm2d1'], massdata['atninterpMEAN'])
res.statistic #0.5245283018867926
res.pvalue # 0.03064707215397448
#_____________

#Attenuance and currents
#Welch doesn't show anything interesting
#filter out periods with 0 collect time

import numpy as np

import pandas as pd


import matplotlib.pyplot as plt


from datetime import datetime

import scipy as SP

import scipy.signal as signal

import seaborn as sns

import statsmodels.tsa.stattools as sts


Hourlyatn = HourlyMAIN[['date_time', 'atn']]

Hourlyatn =Hourlyatn[Hourlyatn['atn'].notna()]

Hourlyatn = Hourlyatn.set_index('date_time')

Hourlyatn = Hourlyatn.resample('H').interpolate()



#Try autocorrelation
#from https://machinelearningmastery.com/gentle-introduction-autocorrelation-partial-autocorrelation/
from pandas import read_csv
from matplotlib import pyplot
from statsmodels.graphics.tsaplots import plot_acf


#Now for current speed

Hourlycurrentspeed = HourlyMAIN[['date_time', 'Speed (cm/s)']]

Hourlycurrentspeed =Hourlycurrentspeed[Hourlycurrentspeed['Speed (cm/s)'].notna()]

Hourlycurrentspeed = Hourlycurrentspeed.set_index('date_time')

Hourlycurrentspeed = Hourlycurrentspeed.resample('H').interpolate()



#try spectray density welch
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html
fig,axs = plt.subplots(4,1, figsize=(8,8), dpi=300)
fs = 1
f, Pxx_den = signal.welch(Hourlyatn['atn'])
axs[0].semilogy(f, Pxx_den)
axs[0].set_xlabel('frequency [per hour]')
axs[0].set_ylabel('PSD [V**2/per hour]')
axs[0].set_title('5 min collect, attenuance MARS 2024')
letter_annotation(axs[0], -.1, 1, 'A')

#try spectray density welch
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html
fs = 1
f, Pxx_den = signal.welch(Hourlycurrentspeed['Speed (cm/s)'])
axs[1].semilogy(f, Pxx_den)
axs[1].set_xlabel('frequency [per hour]')
axs[1].set_ylabel('PSD [V**2/per hour]')
axs[1].set_title('5 min collect, current speed MARS 2024')
letter_annotation(axs[1], -.1, 1, 'B')

plot_acf(Hourlyatn['atn'], lags=72, ax=axs[2])
axs[2].set_ylim(-0.2, 1.01) 
axs[2].set_title('5 minute collect, attenuance MARS 2024')
letter_annotation(axs[2], -.1, 1, 'C')

plot_acf(Hourlycurrentspeed['Speed (cm/s)'], lags=72, ax=axs[3])
axs[3].set_ylim(-0.2, 1.01) 
axs[3].set_title('5 minute collect, current speed MARS 2024')
letter_annotation(axs[3], -.1, 1, 'D')

plt.tight_layout()

plt.savefig('G:/Analyses/2024 MARS deployment/periodicity.pdf')


#cross correlations

#DO THESE IN R I don't have a lot of confidence in this' to find lag

#AND DO SLIDING WINDOW ROUTINE I MADE FOR THE PYROSOME PAPER
# fields = ['date_time','Speed (cm/s)','atn'] # mdct is datetime 
# x = HourlyMAIN[fields]
# df_new = x.dropna()

# def ccf_values(series1, series2):
#     p = series1
#     q = series2
#     p = (p - np.mean(p)) / (np.std(p) * len(p))
#     q = (q - np.mean(q)) / (np.std(q))  
#     c = np.correlate(p, q, 'full')
#     return c
    
# ccf_ielts = ccf_values(df_new['Speed (cm/s)'], df_new['atn'])
# ccf_ielts


# lags = signal.correlation_lags(len(df_new['Speed (cm/s)']), len(df_new['atn']))

# def ccf_plot(lags, ccf):
#     fig, ax =plt.subplots(figsize=(9, 6))
#     ax.plot(lags, ccf)
#     ax.axvline(x = 0, color = 'black', lw = 1)
#     ax.axhline(y = 0, color = 'black', lw = 1)
#     ax.set(ylim = [-.5, .5])
#     ax.set_ylabel('Correlation Coefficients', weight='bold', 
#     fontsize = 12)
#     ax.set(xlim = [-30, 30])
#     ax.set_xlabel('Time Lags', weight='bold', fontsize = 12)
#     plt.legend()
    
# ccf_plot(lags, ccf_ielts)
# plt.savefig('G:/Analyses/2024 MARS deployment/No lag between current speed and attenuance.pdf')