# -*- coding: utf-8 -*-
"""
Created on Mon Jan 13 09:17:33 2025

@author: chuffard
"""


import pandas as pd
import matplotlib.pyplot as plt
import scipy
import seaborn as sns


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

# # Create an observer object representing your location
# observer = ephem.Observer()
# observer.lat = '36.71171'  # Latitude
# observer.lon = '-122.1860967'  # Longitude

SESfpAggpresenceabsence = pd.read_excel('G:/Analyses/OTC conference/FecalAgg500.xlsx', sheet_name = "FPaggHourly")
SESfpAggpresenceabsence.dtypes
SESfpAggpresenceabsence = SESfpAggpresenceabsence.dropna()

Moon = pd.read_excel('G:/Analyses/OTC conference/MoonFraction2024.xlsx', sheet_name = "SESdeploy")


HourlySES = pd.read_csv('G:/Deployments/MARS-2024/Current meters/SES current meter 2024/HourlyMAIN.csv',parse_dates=True)
HourlySES.dtypes

SES_OTC= pd.merge(HourlySES, SESfpAggpresenceabsence, on=['YMD','hour'], how='outer')

SES_OTC.dtypes
SES_OTC=SES_OTC[['YMD','hour', 'presence agg', 'presence big fecal', 'atn', 'Speed (cm/s)', 'TempC']]
SES_OTC = SES_OTC.dropna()
SES_OTC = SES_OTC.apply(pd.to_numeric, errors='coerce')
SES_OTC['date'] =  pd.to_datetime(SES_OTC['YMD'], format='%y%m%d')
SES_OTC = SES_OTC[(SES_OTC['date'] > "2024-06-27") & (SES_OTC['date'] < "2024-07-29")]

SES_OTC= pd.merge(SES_OTC, Moon, on=['date'], how='outer')
SES_OTC = SES_OTC.dropna()

#CALC atn flux Estapa et al 2023: This uses the rate of change to be atn accumulated over five min (from 0 to image time)
SES_OTC['F_atn'] = SES_OTC['atn']/(5/60*24)

    

#calc POC flux from atn flux Estapa et al 2023
SES_OTC['POCmodel_mg']= (10**3.4 * SES_OTC['F_atn']**0.94) * 2 #because collection area in 0.5 m^2
SES_OTC['POCmodel_g']= SES_OTC['POCmodel_mg']/1000



summaryAGG = SES_OTC.groupby('hour')['presence agg'].agg(['mean', 'std'])
summaryAGG= summaryAGG.reset_index()
summaryFEC = SES_OTC.groupby('hour')['presence big fecal'].agg(['mean', 'std'])
summaryFEC= summaryFEC.reset_index()
summaryATN = SES_OTC.groupby('hour')['atn'].agg(['mean', 'std'])
summaryATN= summaryATN.reset_index()
summaryCUR = SES_OTC.groupby('hour')['Speed (cm/s)'].agg(['mean', 'std'])
summaryCUR= summaryCUR.reset_index()
summaryTemp = SES_OTC.groupby('hour')['TempC'].agg(['mean', 'std'])
summaryTemp= summaryTemp.reset_index()
summaryPOC_g = SES_OTC.groupby('hour')['POCmodel_g'].agg(['mean', 'std'])
summaryPOC_g= summaryPOC_g.reset_index()



summaryAGGdaily = SES_OTC.groupby('date')['presence agg'].agg(['mean', 'std'])
summaryAGGdaily= summaryAGGdaily.reset_index()
summaryFECdaily = SES_OTC.groupby('date')['presence big fecal'].agg(['mean', 'std'])
summaryFECdaily= summaryFECdaily.reset_index()
summaryATNdaily = SES_OTC.groupby('date')['atn'].agg(['mean', 'std'])
summaryATNdaily= summaryATNdaily.reset_index()
summaryCURdaily = SES_OTC.groupby('date')['Speed (cm/s)'].agg(['mean', 'std'])
summaryCURdaily= summaryCURdaily.reset_index()
summaryTempdaily = SES_OTC.groupby('date')['TempC'].agg(['mean', 'std'])
summaryTempdaily= summaryTempdaily.reset_index()
summaryMoondaily = SES_OTC.groupby('date')['MoonFraction'].agg(['mean', 'std'])
summaryMoondaily= summaryMoondaily.reset_index()
summaryPOC_g_daily = SES_OTC.groupby('date')['POCmodel_g'].agg(['mean', 'std'])
summaryPOC_g_daily= summaryPOC_g_daily.reset_index()

SES_prelong = SES_OTC.set_index('date')

SES_prelongDaily =  SES_OTC.groupby([pd.Grouper(key='date', freq='D')]).mean() 
#subset
SES_prelong= SES_prelongDaily[["presence agg", 'presence big fecal', 'POCmodel_g']]
SES_OTCLong = SES_prelong.melt(ignore_index=False).reset_index()
SES_prelongplot = SES_prelong.reset_index()
SES_prelongplot['total'] = SES_prelongplot['presence agg']+ SES_prelongplot['presence big fecal']
SES_prelongplot['percAGG'] = SES_prelongplot['presence agg']/ SES_prelongplot['total']
SES_prelongplot['percFEC'] = SES_prelongplot['presence big fecal']/ SES_prelongplot['total']


######################ADD ERROR BARS
plt.rcParams.update({'font.size': 9})
fig,axs = plt.subplots(5, figsize=(6,11), dpi=300, sharex=True)

#axs[0].errorbar(summaryATN['hour'], summaryATN['mean'], yerr=summaryATN['std'], fmt='o' ,ecolor='0.8',c='black')
axs[0].plot(summaryPOC_g_daily['date'], summaryPOC_g_daily['mean'], color='black')
axs[0].set_ylabel('POC flux \n(g C m$^{-2}$ d$^{-1}$)')
axs[0].set_xlabel('')
axs[0].fill_between(summaryPOC_g_daily['date'], summaryPOC_g_daily['mean'] - summaryPOC_g_daily['std'], summaryPOC_g_daily['mean'] + summaryPOC_g_daily['std'], alpha=0.3, color='lightblue')
letter_annotation(axs[0], -.2, 1, 'a')

#axs[1].errorbar(summaryAGG['hour'], summaryAGG['mean'], yerr=summaryAGG['std'], fmt='o' ,ecolor='0.8',c='black')
axs[1].plot(summaryAGGdaily['date'], summaryAGGdaily['mean'], color='black')
axs[1].set_ylabel('aggregate\npresence-absence \n(average)')
axs[1].set_xlabel('')
axs[1].fill_between(summaryAGGdaily['date'], summaryAGGdaily['mean'] - summaryAGGdaily['std'], summaryAGGdaily['mean'] + summaryAGGdaily['std'], alpha=0.3, color='lightblue')
letter_annotation(axs[1], -.2, 1, 'b')

#axs[2].errorbar(summaryFEC['hour'], summaryFEC['mean'], yerr=summaryFEC['std'], fmt='o' ,ecolor='0.8',c='black')
axs[2].plot(summaryFECdaily['date'], summaryFECdaily['mean'], color='black')
axs[2].set_ylabel('large fecal pellet\n presence-absence \n(average)')
axs[2].set_xlabel('')
axs[2].fill_between(summaryFECdaily['date'], summaryFECdaily['mean'] - summaryFECdaily['std'], summaryFECdaily['mean'] + summaryFECdaily['std'], alpha=0.3, color='lightblue')
letter_annotation(axs[2], -.2, 1, 'c')


#axs[3].errorbar(summaryCUR['hour'], summaryCUR['mean'], yerr=summaryCUR['std'], fmt='o' ,ecolor='0.8',c='black')
axs[3].plot(summaryCURdaily['date'], summaryCURdaily['mean'], color='black')
axs[3].set_ylabel('current speed (cm/s)')
axs[3].set_xlabel('')
axs[3].fill_between(summaryCURdaily['date'], summaryCURdaily['mean'] - summaryCURdaily['std'], summaryCURdaily['mean'] + summaryCURdaily['std'], alpha=0.3, color='lightblue')
letter_annotation(axs[3], -.2, 1, 'd')


#axs[4].errorbar(summaryTemp['hour'], summaryTemp['mean'], yerr=summaryTemp['std'], fmt='o' ,ecolor='0.8',c='black')
axs[4].plot(SES_OTC['date'], SES_OTC['MoonFraction'], color='black')
axs[4].set_ylabel('lunar illumination')
axs[4].set_xlabel('')
#axs[4].fill_between(summaryTempdaily['date'], summaryTempdaily['mean'] - summaryTempdaily['std'], summaryTempdaily['mean'] + summaryTempdaily['std'], alpha=0.3, color='lightblue')
letter_annotation(axs[4], -.2, 1, 'e')
plt.xticks(rotation=90)  
plt.tight_layout()
plt.savefig('G:/Analyses/OTC conference/daily.pdf')
plt.savefig('G:/Analyses/OTC conference/daily.png')
SES_OTC.to_csv('G:/Analyses/OTC conference/SES_OTC.csv') 




######################ADD ERROR BARS
plt.rcParams.update({'font.size': 9})
fig,axs = plt.subplots(4, figsize=(6,11), dpi=300)

#axs[0].errorbar(summaryATN['hour'], summaryATN['mean'], yerr=summaryATN['std'], fmt='o' ,ecolor='0.8',c='black')
axs[0].plot(summaryPOC_g['hour'], summaryPOC_g['mean'], color='black')
axs[0].set_ylabel('POC flux \n(g C m$^{-2}$ d$^{-1}$)')
axs[0].set_xlabel('hour')
axs[0].fill_between(summaryPOC_g['hour'], summaryPOC_g['mean'] - summaryPOC_g['std'], summaryPOC_g['mean'] + summaryPOC_g['std'], alpha=0.3, color='lightblue')
letter_annotation(axs[0], -.2, 1, 'a')

#axs[1].errorbar(summaryAGG['hour'], summaryAGG['mean'], yerr=summaryAGG['std'], fmt='o' ,ecolor='0.8',c='black')
axs[1].plot(summaryAGG['hour'], summaryAGG['mean'], color='black')
axs[1].set_ylabel('aggregate\npresence-absence \n(average)')
axs[1].set_xlabel('hour')
axs[1].fill_between(summaryAGG['hour'], summaryAGG['mean'] - summaryAGG['std'], summaryAGG['mean'] + summaryAGG['std'], alpha=0.3, color='lightblue')
letter_annotation(axs[1], -.2, 1, 'b')

#axs[2].errorbar(summaryFEC['hour'], summaryFEC['mean'], yerr=summaryFEC['std'], fmt='o' ,ecolor='0.8',c='black')
axs[2].plot(summaryFEC['hour'], summaryFEC['mean'], color='black')
axs[2].set_ylabel('large fecal pellet\n presence-absence \n(average)')
axs[2].set_xlabel('hour')
axs[2].fill_between(summaryFEC['hour'], summaryFEC['mean'] - summaryFEC['std'], summaryFEC['mean'] + summaryFEC['std'], alpha=0.3, color='lightblue')
letter_annotation(axs[2], -.2, 1, 'c')


#axs[3].errorbar(summaryCUR['hour'], summaryCUR['mean'], yerr=summaryCUR['std'], fmt='o' ,ecolor='0.8',c='black')
axs[3].plot(summaryCUR['hour'], summaryCUR['mean'], color='black')
axs[3].set_ylabel('current speed (cm/s)')
axs[3].set_xlabel('hour')
axs[3].fill_between(summaryCUR['hour'], summaryCUR['mean'] - summaryCUR['std'], summaryCUR['mean'] + summaryCUR['std'], alpha=0.3, color='lightblue')
letter_annotation(axs[3], -.2, 1, 'd')


# #axs[4].errorbar(summaryTemp['hour'], summaryTemp['mean'], yerr=summaryTemp['std'], fmt='o' ,ecolor='0.8',c='black')
# axs[4].plot(summaryTemp['hour'], summaryTemp['mean'], color='black')
# axs[4].set_ylabel('temperature (C)')
# axs[4].set_xlabel('hour')
# axs[4].fill_between(summaryTemp['hour'], summaryTemp['mean'] - summaryTemp['std'], summaryTemp['mean'] + summaryTemp['std'], alpha=0.3, color='lightblue')
# letter_annotation(axs[4], -.2, 1, 'E')
plt.tight_layout()
plt.savefig('G:/Analyses/OTC conference/by_hour.pdf')
plt.savefig('G:/Analyses/OTC conference/by_hour.png')





plt.rcParams.update({'font.size': 9})
fig,axs = plt.subplots(3,2, figsize=(7,7), dpi=300)

slopeAGG, interceptAGG, r_valueAGG, p_valueAGG, std_errAGG = scipy.stats.linregress(summaryAGGdaily['mean'], summaryPOC_g_daily['mean'])
print("slopePOCAGG:", slopeAGG)
print("r_valuePOCAGG:", r_valueAGG)
print("p_valuePOCAGG:", p_valueAGG)
lineAGG = slopeAGG*summaryAGGdaily['mean']+interceptAGG
#plt.plot(summaryPOC_g_daily['mean'], line, 'r', label='y={:.2f}x+{:.2f}'.format(slopeAGG,interceptAGG))
#end

axs[0,0].errorbar(summaryAGGdaily['mean'], summaryPOC_g_daily['mean'],yerr=summaryPOC_g_daily['std'], fmt='o' ,ecolor='0.8',c='black', alpha = 0.5)
axs[0,0].set_xlabel('aggregate presence-absence \n(average)')
axs[0,0].set_ylabel('POC flux \n(g C m$^{-2}$ d$^{-1}$)')
axs[0,0].plot(summaryAGGdaily['mean'], lineAGG, 'r', label='y={:.2f}x+{:.2f}'.format(slopeAGG,interceptAGG))
letter_annotation(axs[0,0], -.4, 1, 'a')

slopeFEC, interceptFEC, r_valueFEC, p_valueFEC, std_errFEC = scipy.stats.linregress(summaryFECdaily['mean'], summaryPOC_g_daily['mean'])
print("slopePOCFEC:", slopeFEC)
print("r_valuePOCFEC:", r_valueFEC)
print("p_valuePOCFEC:", p_valueFEC)
lineFEC = slopeFEC*summaryFECdaily['mean']+interceptFEC
#plt.plot(summaryPOC_g_daily['mean'], line, 'r', label='y={:.2f}x+{:.2f}'.format(slopeFEC,interceptFEC))
#end

axs[0,1].errorbar(summaryFECdaily['mean'], summaryPOC_g_daily['mean'],yerr=summaryPOC_g_daily['std'], fmt='o' ,ecolor='0.8',c='black', alpha = 0.5)
axs[0,1].set_xlabel('large fecal pellet \n presence-absence \n(average)')
axs[0,1].set_ylabel('POC flux \n(g C m$^{-2}$ d$^{-1}$)')
axs[0,1].plot(summaryFECdaily['mean'], lineFEC, 'r', label='y={:.2f}x+{:.2f}'.format(slopeFEC,interceptFEC))
letter_annotation(axs[0,1], -.4, 1, 'b')


# Plot the mean line
#axs[0,0].scatter(summaryCUR['mean'], summaryAGG['mean'], color='black')
axs[1,0].errorbar(summaryCURdaily['mean'], summaryAGGdaily['mean'],yerr=summaryAGGdaily['std'], fmt='o' ,ecolor='0.8',c='black', alpha = 0.5)
axs[1,0].set_ylabel('aggregate presence-absence \n(average)')
axs[1,0].set_xlabel('current speed (cm/s)')
letter_annotation(axs[1,0], -.4, 1, 'c')

#axs[0,1].scatter(summaryCUR['mean'], summaryFEC['mean'], color='black')
axs[1, 1].errorbar(summaryCURdaily['mean'], summaryFECdaily['mean'],yerr=summaryFECdaily['std'], fmt='o' ,ecolor='0.8',c='black', alpha = 0.5)
axs[1,1].set_ylabel('large fecal pellet \n presence-absence \n(average)')
axs[1,1].set_xlabel('current speed (cm/s)')
letter_annotation(axs[1,1], -.4, 1, 'd')


#axs[1,0].scatter(summaryTemp['mean'], summaryAGG['mean'], color='black')
axs[2,0].errorbar(summaryMoondaily['mean'], summaryAGGdaily['mean'],yerr=summaryAGGdaily['std'], fmt='o' ,ecolor='0.8',c='black', alpha = 0.5)
axs[2,0].set_ylabel('aggregate presence-absence \n(average)')
axs[2,0].set_xlabel('lunar illumination')
letter_annotation(axs[2,0], -.4, 1, 'e')


#axs[1,1].scatter(summaryTemp['mean'], summaryFEC['mean'], color='black')
axs[2,1].errorbar(summaryMoondaily['mean'], summaryFECdaily['mean'],yerr=summaryFECdaily['std'], fmt='o' ,ecolor='0.8',c='black', alpha = 0.5)
axs[2,1].set_ylabel('large fecal pellet \n presence-absence \n(average)')
axs[2,1].set_xlabel('lunar illumination')
letter_annotation(axs[2,1], -.4, 1, 'f')


plt.tight_layout(pad=1.5)
plt.savefig('G:/Analyses/OTC conference/scatterplots.pdf')
plt.savefig('G:/Analyses/OTC conference/scatterplots.png')

SES_prelongDaily.dtypes
# YMD                   float64
# hour                  float64
# presence agg          float64
# presence big fecal    float64
# atn                   float64
# Speed (cm/s)          float64
# TempC                 float64
# Day                   float64
# MoonFraction          float64
# F_atn                 float64
# POCmodel_mg           float64
# POCmodel_g            float64
# dtype: object


# Create subplots
fig, axes = plt.subplots(3, 2, figsize=(7, 7), dpi=300)

# Plot on the first subplot
sns.regplot(data=SES_prelongDaily, x='presence agg', y='POCmodel_g', ax=axes[0,0],scatter_kws={'alpha': 0.5})
axes[0,0].set_xlabel('aggregate \n presence-absence \n(average)')
axes[0,0].set_ylabel('POC flux \n(g C m$^{-2}$ d$^{-1}$)')
letter_annotation(axes[0,0], -.4, 1, 'a')


sns.regplot(data=SES_prelongDaily, x='presence big fecal', y='POCmodel_g', ax=axes[0,1],scatter_kws={'alpha': 0.5})
axes[0,1].set_xlabel('large fecal pellet \n presence-absence \n(average)')
axes[0,1].set_ylabel('POC flux \n(g C m$^{-2}$ d$^{-1}$)')
letter_annotation(axes[0,1], -.4, 1, 'b')

sns.regplot(data=SES_prelongDaily, y='presence agg', x='MoonFraction', order=2,ax=axes[1,0],scatter_kws={'alpha': 0.5})
axes[1,0].set_ylabel('aggregate \n presence-absence \n(average)')
axes[1,0].set_xlabel('lunar illumination')
letter_annotation(axes[1,0], -.4, 1, 'c')

sns.regplot(data=SES_prelongDaily, y='presence big fecal', x='MoonFraction', order=2,ax=axes[1,1],scatter_kws={'alpha': 0.5})
axes[1,1].set_ylabel('large fecal pellet \n presence-absence \n(average)')
axes[1,1].set_xlabel('lunar illumination')
letter_annotation(axes[1,1], -.4, 1, 'd')

sns.regplot(data=SES_prelongDaily, y='presence agg', x='Speed (cm/s)', ax=axes[2,0],scatter_kws={'alpha': 0.5})
axes[2,0].set_ylabel('aggregate \n presence-absence \n(average)')
axes[2,0].set_xlabel('current speed (cm/s)')
letter_annotation(axes[2,0], -.4, 1, 'e')


sns.regplot(data=SES_prelongDaily, y='presence big fecal', x='Speed (cm/s)', ax=axes[2,1],scatter_kws={'alpha': 0.5})
axes[2,1].set_ylabel('large fecal pellet \n presence-absence \n(average)')
axes[2,1].set_xlabel('current speed (cm/s)')
letter_annotation(axes[2,1], -.4, 1, 'f')



plt.tight_layout(pad=2)
plt.savefig('G:/Analyses/OTC conference/regplots.pdf')
plt.savefig('G:/Analyses/OTC conference/regplots.png')