# -*- 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 statistics
#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

import numpy as np
from scipy import stats

from matplotlib.backends.backend_pdf import PdfPages

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

# 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()

SEStrim['atnflux']=SEStrim["atn"] *0.0034722222222222*12



#fix weird time zone stuff
# for col in Main_Hourly_SES_MARS.select_dtypes(include=['datetime64[ns, UTC]']).columns:
#     Main_Hourly_SES_MARS[col] = Main_Hourly_SES_MARS[col].apply(lambda x: x.tz_localize(None))

#run a loop that uses running means from 275 to 1500 hours, and see how many triggers it comes up with

count=0
R_atncount=[]
R_atnmeancount=[]
n_count=[]
dfcount=[]
SErcount=[]
CIhighr=[]
trigger=[]
trigger_count = []
wait_count = []
triggerlist=[]

#4 days to 18 days
for x in range(288, 1512):
        data =SEStrim
        variablename = str(x)
        data[variablename] = data['atn'].rolling(x, min_periods=50).mean()
        R_atncount = data['atn'].rolling(x, min_periods=50).mean()
        R_atnmeancount = R_atncount.mean()
        n_count = len(R_atncount)
        dfcount = n_count-1
        SErcount = stats.sem(R_atncount.dropna())
        CIlevel = 0.9999
        CI99atnrcount = stats.t.interval(CIlevel, dfcount, R_atnmeancount, scale = SErcount)
        CIhighr = float(CI99atnrcount[1])
        triggercol = str(x)
        data[triggercol]=np.where(data[variablename] >= CIhighr, 'Trigger', 'Wait')
        triggercount = data[triggercol].value_counts()['Trigger']
        waitcount = data[triggercol].value_counts()['Wait']
        totalcount = triggercount+waitcount
        triggerlist.append(triggercount)
    
triggerlistdf = pd.DataFrame(triggerlist)

triggerlistdf= triggerlistdf.rename(columns={0: "number_of_triggers"})
triggerlistdf = triggerlistdf.reset_index()
triggerlistdf['rollingmean_count'] = triggerlistdf['index']+288
#triggerlistdf= triggerlistdf.rename(columns={"index": "rollingmean_count"})
triggerlistdf.to_csv("triggerlistdf.csv")

sns.lineplot(data=triggerlistdf, x='rollingmean_count', y='number_of_triggers')

#find the rolling mean with the fewest triggers
rowwithlottrigger = triggerlistdf[triggerlistdf.number_of_triggers==triggerlistdf.number_of_triggers.min()]
rollingmeanUSE = rowwithlottrigger.rollingmean_count.min()
print(rollingmeanUSE)
#406
#make a new time series of this rolling mean
SEStrim['rolling_atn_406'] = SEStrim['atn'].rolling(rollingmeanUSE, min_periods=200).mean()

#Now check the slope of this time series
#subset columns
SESMARS406atn = SEStrim[["atn", "rolling_atn_406", 'atnflux']]

#find proxy for slope (percent change from one sample to next)
SESMARS406atn['diff'] = SESMARS406atn['rolling_atn_406'].diff()
#SESMARS406atn['Percentchange'] = SESMARS406atn['rolling_atn_406'].astype(float).pct_change(1)

#subset to remove first 406 rows (to eliminate artifacts of using rolling mean)
SESslope406TRIM = SESMARS406atn.iloc[406:, ]

#find max slope after this
trigger406slope = max(SESslope406TRIM['diff'])
print(trigger406slope)
#0.005114922847290648

#make that new max the trigger
SESslope406TRIM["trigger406slope"] = np.where(SESslope406TRIM['diff'] >= trigger406slope, 'TriggerSlope', 'WaitSlope')

#this is a little high and picks up the peak. maybe we can get a little before the peak
#trigger 80% of max slope
trigger80pc406slope = 0.8 * trigger406slope
print(trigger80pc406slope)
#0.004091938277832519

#list which samples are trigger vs wait
SESslope406TRIM["trigger80pc406slope"] = np.where(SESslope406TRIM['diff'] >= trigger80pc406slope, 'TriggerSlope', 'WaitSlope')

#plot to verify.
sns.scatterplot(data=SESslope406TRIM, x="Date_time", y="atn", hue="trigger80pc406slope")
plt.xticks(rotation=45);
plt.show()

#how does this compare to hourly and daily means of attenuance?
#find hourly average of SES time series
hourly_SESatn = SESslope406TRIM['atn'].resample('1H').mean()

#turn into data frame
hourly_SESatndf = pd.DataFrame([hourly_SESatn]).T

#find daily average of SES time series
daily_SESatn = SESslope406TRIM['atn'].resample('1D').mean()

#turn into data frame
daily_SESatnatndf = pd.DataFrame([daily_SESatn]).T



#find hourly average of SES time series
hourly_SESatnFLUX = SESslope406TRIM['atnflux'].resample('1H').mean()

#turn into data frame
hourly_SESatndfFLUX = pd.DataFrame([hourly_SESatnFLUX]).T

#find daily average of SES time series
daily_SESatnFLUX = SESslope406TRIM['atnflux'].resample('1D').mean()

#turn into data frame
daily_SESatnatndfFLUX = pd.DataFrame([daily_SESatnFLUX]).T




testtrim5min = SESslope406TRIM.loc['2023-05-08':'2023-06-13']  
testtrim1H = hourly_SESatndfFLUX.loc['2023-05-08':'2023-06-13']  
testtrim1D = daily_SESatnatndf.loc['2023-05-08':'2023-06-13']  


fig, ax = plt.subplots()
sns.scatterplot(data=testtrim5min , x="Date_time", y="atnflux", ax=ax, color = 'gray', s = 10)
sns.lineplot(data=testtrim1H, x="Date_time", y="atnflux", ax=ax, color = "orange")
plt.xticks(rotation=45);





#How long does the elevated period last?



























#What's the POC flux from this?
#Using Equation 6 and Tabe 3's coefficients in Table 3 for FatnKd
#https://aslopubs.onlinelibrary.wiley.com/doi/full/10.1002/lom3.10592
a2 =3.4 
a1 = 0.94

SESslope406TRIM
SESslope406TRIM['FPOC_SES'] = (10**a2) * (SESslope406TRIM['atn']**a1)
testtrim = SESslope406TRIM.loc['2023-05-08':'2023-06-13']  


#plot to verify. See if it looks like the plots Colleen made earlier
sns.lineplot(data=testtrim, x="Date_time", y="FPOC_SES", hue="trigger80pc406slope")
plt.xticks(rotation=45);
plt.show()



#plot to verify. See if it looks like the plots Colleen made earlier
sns.lineplot(data=testtrim, x="datetime", y="FPOC_SES")
plt.xticks(rotation=45);
plt.show()