# -*- coding: utf-8 -*- """ Created on Mon Oct 5 07:32:12 2026 @author: chuffard """ #%% Now convert to POC flux- Station M data remember collect time is different for each one import numpy as np import pandas as pd import matplotlib.pyplot as plt import os , sys import cv2 import glob from datetime import datetime import os.path import warnings from pathlib import Path import math import seaborn as sns import matplotlib.dates as mdates #%% atn60 = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/SES_out_atn60.csv') SEScolltime = 0.5 #in days # Extract the pattern into a new column atn60['file_name'] = atn60['file_name'].astype(str).str.strip() atn60[['date', 'time']] = atn60['file_name'].str.extract(r'(\d{6})-(\d{6})') # Convert to actual datetime objects for analysis atn60['timestamp'] = pd.to_datetime(atn60['date'] + atn60['time'], format='%y%m%d%H%M%S') atn60.columns atn60["DateTime"] = pd.to_datetime(atn60['timestamp']) #convert atn to POC flux #from Estapa et al. 2024 optical sediment trap calibration. also accounting for image and collect area since those aren't in our og atn calculation SEStraparea = 0.5; #imgarea = length(find(isfinite(bgr(:))))./(0.07194.^2)./(1000000.^2); % pix / (pix/um^2) / (um^2/m^2)#this was in Meg's example code but it's not for SES #This is adapted to pyton from matlab code Meg gave us SESimgarea = math.pi*(25/1000/2)**2 # diameter of lower end of funnel is 25 mm, compute area in m^2 atn60['F_atn']= atn60['atn']/SEScolltime*SESimgarea/SEStraparea #These are from SES at Station M in 2024 paper atn60['POCmodel_mg']= 10**3.4 * atn60['F_atn']**0.94 atn60['POCmodel_g']= atn60['POCmodel_mg']/1000 atn60['POCmodel_g'].mean()#(0.2753247346815628) atn60['POCmodel_g'].plot() plt.show() atn60.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/atn60_toPOCplux.csv') #### #%%Pulse 61 atn61 = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/SES_out_atn61.csv') SEScolltime = 3/24 #in days # Extract the pattern into a new column atn61['file_name'] = atn61['file_name'].astype(str).str.strip() atn61[['date', 'time']] = atn61['file_name'].str.extract(r'(\d{6})-(\d{6})') # Convert to actual datetime objects for analysis atn61['timestamp'] = pd.to_datetime(atn61['date'] + atn61['time'], format='%y%m%d%H%M%S') atn61.columns atn61["DateTime"] = pd.to_datetime(atn61['timestamp']) #convert atn to POC flux #from Estapa et al. 2024 optical sediment trap calibration. also accounting for image and collect area since those aren't in our og atn calculation SEStraparea = 0.5; #imgarea = length(find(isfinite(bgr(:))))./(0.07194.^2)./(1000000.^2); % pix / (pix/um^2) / (um^2/m^2)#this was in Meg's example code but it's not for SES #This is adapted to pyton from matlab code Meg gave us SESimgarea = math.pi*(25/1000/2)**2 # diameter of lower end of funnel is 25 mm, compute area in m^2 atn61['F_atn']= atn61['atn']/SEScolltime*SESimgarea/SEStraparea #These are from SES at Station M in 2024 paper atn61['POCmodel_mg']= 10**3.4 * atn61['F_atn']**0.94 atn61['POCmodel_g']= atn61['POCmodel_mg']/1000 atn61['POCmodel_g'].mean()#(0.2753247346815628) atn61['POCmodel_g'].plot() plt.show() atn61.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/atn61_toPOCplux.csv') #%% Pulse 62 atn62 = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/SES_out_atn62.csv') SEScolltime = 4/24 #in days # Extract the pattern into a new column atn62['file_name'] = atn62['file_name'].astype(str).str.strip() atn62[['date', 'time']] = atn62['file_name'].str.extract(r'(\d{6})-(\d{6})') # Convert to actual datetime objects for analysis atn62['timestamp'] = pd.to_datetime(atn62['date'] + atn62['time'], format='%y%m%d%H%M%S') atn62.columns atn62["DateTime"] = pd.to_datetime(atn62['timestamp']) #convert atn to POC flux #from Estapa et al. 2024 optical sediment trap calibration. also accounting for image and collect area since those aren't in our og atn calculation SEStraparea = 0.5; #imgarea = length(find(isfinite(bgr(:))))./(0.07194.^2)./(1000000.^2); % pix / (pix/um^2) / (um^2/m^2)#this was in Meg's example code but it's not for SES #This is adapted to pyton from matlab code Meg gave us SESimgarea = math.pi*(25/1000/2)**2 # diameter of lower end of funnel is 25 mm, compute area in m^2 atn62['F_atn']= atn62['atn']/SEScolltime*SESimgarea/SEStraparea #These are from SES at Station M in 2024 paper atn62['POCmodel_mg']= 10**3.4 * atn62['F_atn']**0.94 atn62['POCmodel_g']= atn62['POCmodel_mg']/1000 atn62['POCmodel_g'].mean()#(0.2753247346815628) atn62['POCmodel_g'].plot() plt.show() atn62.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/atn62_toPOCplux.csv') #%% Pulse 64 atn64 = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/SES_out_atn64.csv') SEScolltime = 4/24 #in days # Extract the pattern into a new column atn64['file_name'] = atn64['file_name'].astype(str).str.strip() atn64[['date', 'time']] = atn64['file_name'].str.extract(r'(\d{6})-(\d{6})') # Convert to actual datetime objects for analysis atn64['timestamp'] = pd.to_datetime(atn64['date'] + atn64['time'], format='%y%m%d%H%M%S') atn64.columns atn64["DateTime"] = pd.to_datetime(atn64['timestamp']) #convert atn to POC flux #from Estapa et al. 2024 optical sediment trap calibration. also accounting for image and collect area since those aren't in our og atn calculation SEStraparea = 0.5; #imgarea = length(find(isfinite(bgr(:))))./(0.07194.^2)./(1000000.^2); % pix / (pix/um^2) / (um^2/m^2)#this was in Meg's example code but it's not for SES #This is adapted to pyton from matlab code Meg gave us SESimgarea = math.pi*(25/1000/2)**2 # diameter of lower end of funnel is 25 mm, compute area in m^2 atn64['F_atn']= atn64['atn']/SEScolltime*SESimgarea/SEStraparea #These are from SES at Station M in 2024 paper atn64['POCmodel_mg']= 10**3.4 * atn64['F_atn']**0.94 atn64['POCmodel_g']= atn64['POCmodel_mg']/1000 atn64['POCmodel_g'].mean()#(0.2753247346815648) atn64['POCmodel_g'].plot() plt.show() atn64.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/atn64_toPOCplux.csv') #%% Pulse 66 atn66 = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/SES_out_atn66.csv') SEScolltime = 4/24 #in days # Extract the pattern into a new column atn66['file_name'] = atn66['file_name'].astype(str).str.strip() atn66[['date', 'time']] = atn66['file_name'].str.extract(r'(\d{6})-(\d{6})') # Convert to actual datetime objects for analysis atn66['timestamp'] = pd.to_datetime(atn66['date'] + atn66['time'], format='%y%m%d%H%M%S') atn66.columns atn66["DateTime"] = pd.to_datetime(atn66['timestamp']) #convert atn to POC flux #from Estapa et al. 2024 optical sediment trap calibration. also accounting for image and collect area since those aren't in our og atn calculation SEStraparea = 0.5; #imgarea = length(find(isfinite(bgr(:))))./(0.07194.^2)./(1000000.^2); % pix / (pix/um^2) / (um^2/m^2)#this was in Meg's example code but it's not for SES #This is adapted to pyton from matlab code Meg gave us SESimgarea = math.pi*(25/1000/2)**2 # diameter of lower end of funnel is 25 mm, compute area in m^2 atn66['F_atn']= atn66['atn']/SEScolltime*SESimgarea/SEStraparea #These are from SES at Station M in 2024 paper atn66['POCmodel_mg']= 10**3.4 * atn66['F_atn']**0.94 atn66['POCmodel_g']= atn66['POCmodel_mg']/1000 atn66['POCmodel_g'].mean()#(0.2753247346815668) atn66['POCmodel_g'].plot() plt.show() atn66.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/atn66_toPOCplux.csv') #%% Pulse 68 atn68 = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/SES_out_atn68.csv') SEScolltime = 4/24 #in days # Extract the pattern into a new column atn68['file_name'] = atn68['file_name'].astype(str).str.strip() atn68[['date', 'time']] = atn68['file_name'].str.extract(r'(\d{6})-(\d{6})') # Convert to actual datetime objects for analysis atn68['timestamp'] = pd.to_datetime(atn68['date'] + atn68['time'], format='%y%m%d%H%M%S') atn68.columns atn68["DateTime"] = pd.to_datetime(atn68['timestamp']) #convert atn to POC flux #from Estapa et al. 2024 optical sediment trap calibration. also accounting for image and collect area since those aren't in our og atn calculation SEStraparea = 0.5; #imgarea = length(find(isfinite(bgr(:))))./(0.07194.^2)./(1000000.^2); % pix / (pix/um^2) / (um^2/m^2)#this was in Meg's example code but it's not for SES #This is adapted to pyton from matlab code Meg gave us SESimgarea = math.pi*(25/1000/2)**2 # diameter of lower end of funnel is 25 mm, compute area in m^2 atn68['F_atn']= atn68['atn']/SEScolltime*SESimgarea/SEStraparea #These are from SES at Station M in 2024 paper atn68['POCmodel_mg']= 10**3.4 * atn68['F_atn']**0.94 atn68['POCmodel_g']= atn68['POCmodel_mg']/1000 atn68['POCmodel_g'].mean()#(0.2753247346815688) atn68['POCmodel_g'].plot() plt.show() atn68.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/atn68_toPOCplux.csv') #%% Pulse 70 atn70 = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/SES_out_atn70.csv') SEScolltime = 2/24 #in days # Extract the pattern into a new column atn70['file_name'] = atn70['file_name'].astype(str).str.strip() atn70[['date', 'time']] = atn70['file_name'].str.extract(r'(\d{6})-(\d{6})') # Convert to actual datetime objects for analysis atn70['timestamp'] = pd.to_datetime(atn70['date'] + atn70['time'], format='%y%m%d%H%M%S') atn70.columns atn70["DateTime"] = pd.to_datetime(atn70['timestamp']) #convert atn to POC flux #from Estapa et al. 2024 optical sediment trap calibration. also accounting for image and collect area since those aren't in our og atn calculation SEStraparea = 0.5; #imgarea = length(find(isfinite(bgr(:))))./(0.07194.^2)./(1000000.^2); % pix / (pix/um^2) / (um^2/m^2)#this was in Meg's example code but it's not for SES #This is adapted to pyton from matlab code Meg gave us SESimgarea = math.pi*(25/1000/2)**2 # diameter of lower end of funnel is 25 mm, compute area in m^2 atn70['F_atn']= atn70['atn']/SEScolltime*SESimgarea/SEStraparea #These are from SES at Station M in 2024 paper atn70['POCmodel_mg']= 10**3.4 * atn70['F_atn']**0.94 atn70['POCmodel_g']= atn70['POCmodel_mg']/1000 atn70['POCmodel_g'].mean()#(0.2753247347015708) atn70['POCmodel_g'].plot() plt.show() atn70.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/atn70_toPOCplux.csv') #%% Pulse 71 atn71 = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/SES_out_atn71.csv') SEScolltime = 2/24 #in days # Extract the pattern into a new column atn71['file_name'] = atn71['file_name'].astype(str).str.strip() atn71[['date', 'time']] = atn71['file_name'].str.extract(r'(\d{6})-(\d{6})') # Convert to actual datetime objects for analysis atn71['timestamp'] = pd.to_datetime(atn71['date'] + atn71['time'], format='%y%m%d%H%M%S') atn71.columns atn71["DateTime"] = pd.to_datetime(atn71['timestamp']) #convert atn to POC flux #from Estapa et al. 2024 optical sediment trap calibration. also accounting for image and collect area since those aren't in our og atn calculation SEStraparea = 0.5; #imgarea = length(find(isfinite(bgr(:))))./(0.07194.^2)./(1000000.^2); % pix / (pix/um^2) / (um^2/m^2)#this was in Meg's example code but it's not for SES #This is adapted to pyton from matlab code Meg gave us SESimgarea = math.pi*(25/1000/2)**2 # diameter of lower end of funnel is 25 mm, compute area in m^2 atn71['F_atn']= atn71['atn']/SEScolltime*SESimgarea/SEStraparea #These are from SES at Station M in 2024 paper atn71['POCmodel_mg']= 10**3.4 * atn71['F_atn']**0.94 atn71['POCmodel_g']= atn71['POCmodel_mg']/1000 atn71['POCmodel_g'].mean()#(0.2753247347115718) atn71['POCmodel_g'].plot() plt.show() atn71.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/atn71_toPOCplux.csv') #%% Pulse 73 atn73 = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/SES_out_atn73.csv') SEScolltime = 2/24 #in days # Extract the pattern into a new column atn73['file_name'] = atn73['file_name'].astype(str).str.strip() atn73[['date', 'time']] = atn73['file_name'].str.extract(r'(\d{6})-(\d{6})') # Convert to actual datetime objects for analysis atn73['timestamp'] = pd.to_datetime(atn73['date'] + atn73['time'], format='%y%m%d%H%M%S') atn73.columns atn73["DateTime"] = pd.to_datetime(atn73['timestamp']) #convert atn to POC flux #from Estapa et al. 2024 optical sediment trap calibration. also accounting for image and collect area since those aren't in our og atn calculation SEStraparea = 0.5; #imgarea = length(find(isfinite(bgr(:))))./(0.07394.^2)./(1000000.^2); % pix / (pix/um^2) / (um^2/m^2)#this was in Meg's example code but it's not for SES #This is adapted to pyton from matlab code Meg gave us SESimgarea = math.pi*(25/1000/2)**2 # diameter of lower end of funnel is 25 mm, compute area in m^2 atn73['F_atn']= atn73['atn']/SEScolltime*SESimgarea/SEStraparea #These are from SES at Station M in 2024 paper atn73['POCmodel_mg']= 10**3.4 * atn73['F_atn']**0.94 atn73['POCmodel_g']= atn73['POCmodel_mg']/1000 atn73['POCmodel_g'].mean()#(0.2753247347315738) atn73['POCmodel_g'].plot() plt.show() atn73.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/atn73_toPOCplux.csv') #%% Pulse 75--DID NOT HOME # atn75 = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/SES_out_atn75.csv') # SEScolltime = 2/24 #in days # # Extract the pattern into a new column # atn75['file_name'] = atn75['file_name'].astype(str).str.strip() # atn75[['date', 'time']] = atn75['file_name'].str.extract(r'(\d{6})-(\d{6})') # # Convert to actual datetime objects for analysis # atn75['timestamp'] = pd.to_datetime(atn75['date'] + atn75['time'], format='%y%m%d%H%M%S') # atn75.columns # atn75["DateTime"] = pd.to_datetime(atn75['timestamp']) # #convert atn to POC flux # #from Estapa et al. 2024 optical sediment trap calibration. also accounting for image and collect area since those aren't in our og atn calculation # SEStraparea = 0.5; # #imgarea = length(find(isfinite(bgr(:))))./(0.07594.^2)./(1000000.^2); % pix / (pix/um^2) / (um^2/m^2)#this was in Meg's example code but it's not for SES # #This is adapted to pyton from matlab code Meg gave us # SESimgarea = math.pi*(25/1000/2)**2 # diameter of lower end of funnel is 25 mm, compute area in m^2 # atn75['F_atn']= atn75['atn']/SEScolltime*SESimgarea/SEStraparea # #These are from SES at Station M in 2024 paper # atn75['POCmodel_mg']= 10**3.4 * atn75['F_atn']**0.94 # atn75['POCmodel_g']= atn75['POCmodel_mg']/1000 # atn75['POCmodel_g'].mean()#(0.2753247547515758) # atn75['POCmodel_g'].plot() # plt.show() # atn75.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/atn75_toPOCplux.csv') #%% Merge rbind # #Station_M_SES_atntoPOC = pd.concat([atn60, atn61, atn62, atn64, atn66, atn68, atn70,atn71,atn73], ignore_index=True) # #dates are a mess. easier to clean in excel # #Station_M_SES_atntoPOC.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/Station_M_SES_atntoPOC.csv') # Station_M_SES_atntoPOCdates = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/Station_M_SES_atntoPOC.csv',parse_dates=['date']) # #average by day # Station_M_SES_atntoPOCdates['pure_date'] = pd.to_datetime(Station_M_SES_atntoPOCdates['date']).dt.date # daily_mean_atn_POCSESStationM = ( # Station_M_SES_atntoPOCdates.groupby(Station_M_SES_atntoPOCdates['pure_date']) # .mean(numeric_only=True) # .rename(columns=lambda col: f"{col}_daily_mean") # .reset_index() # ) # column_list = daily_mean_atn_POCSESStationM.columns.tolist() # print(column_list) # # column_list = daily_mean_atn_POCSESStationM.columns.tolist() # # print(column_list) # #['pure_date', 'Unnamed: 0_daily_mean', 'atn_daily_mean', 'year_daily_mean', 'month_daily_mean', 'day_daily_mean', 'hour_daily_mean', 'F_atn_daily_mean', 'POCmodel_mg_daily_mean', 'POCmodel_g_daily_mean'] # daily_mean_atn_POCSESStationM['pure_date'] = pd.to_datetime(daily_mean_atn_POCSESStationM['pure_date']).dt.date # StationMDailyPlotTSd = daily_mean_atn_POCSESStationM.set_index('pure_date') # full_range = pd.date_range(start=StationMDailyPlotTSd.index.min(), end=StationMDailyPlotTSd.index.max(), freq='D') # # Reindex the DataFrame to include missing dates, leaving values blank (NaN) # StationMDailyPlotTS_padded = StationMDailyPlotTSd.reindex(full_range) # StationMDailyPlotTS_padded.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/StationMDailyPlotTS_padded.csv') # # Plot the 'value' column directly from the padded DataFrame # StationMDailyPlotTS_padded['POCmodel_mg_daily_mean'].plot(kind='line', title='POC flux modeled from non-red atn INCL CLOUDY') StationMDailyPlotTS_NO_CLOUDY = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/StationMDailyPlotTS_NO_CLOUDY.csv') #StationMDailyPlotTS_NO_CLOUDY['pure_date'] = pd.to_datetime(StationMDailyPlotTS_NO_CLOUDY['date']).dt.date StationMDailyPlotTS_NO_CLOUDY['pure_date'] = pd.to_datetime(StationMDailyPlotTS_NO_CLOUDY['date']).dt.strftime("%Y/%m/%d") StationMDailyPlotTS_NO_CLOUDY['pure_date'] = pd.to_datetime(StationMDailyPlotTS_NO_CLOUDY['pure_date']).dt.date StationMDailyPlotTS = StationMDailyPlotTS_NO_CLOUDY.set_index('pure_date') #Trim first 14 d of each deployment because of potential cloudy periods full_range = pd.date_range(start=StationMDailyPlotTS.index.min(), end=StationMDailyPlotTS.index.max(), freq='D') StationMDailyPlotTS['POCmodel_mg_daily_mean'].plot(kind='line', title='POC flux modeled from non-red atn') # Label axes and show grid plt.xlabel('Date') plt.ylabel('POC flux mg C m-2 d-1 from atn SES') plt.grid(True) # Display the plot plt.show() #%% import POC flux from mclane traps POCfluxTraps = pd.read_csv('Z:/STORE/Lab wrap up/final organization/Crissy working/Data papers/POC flux-DONE/FINAL SUBMITTED/Station_M_Particulate_fluxes_DATA_updated_units.csv') POCfluxTraps['pure_date'] = pd.to_datetime(POCfluxTraps['Collect_date']).dt.strftime("%Y/%m/%d") POCfluxTraps['pure_date'] = pd.to_datetime(POCfluxTraps['pure_date']).dt.date POCfluxTraps = POCfluxTraps.set_index('pure_date') POC_STation_M_SES_and_McLane = pd.merge(StationMDailyPlotTS, POCfluxTraps, left_index=True, right_index=True) column_list = POC_STation_M_SES_and_McLane.columns.tolist() print(column_list) # ['date', 'atn_daily_mean', 'year_daily_mean', 'month_daily_mean', 'day_daily_mean', 'hour_daily_mean', 'F_atn_daily_mean', # 'POCmodel_mg_daily_mean', 'POCmodel_g_daily_mean', 'Elim14dCloudy', 'PULSE', 'Station', 'Collect_date', 'Cup_number', # 'Duration open', 'Latitude', 'Longitude', 'Source of Lat Long', 'Bottom depth', 'Source of bottom depth', '50_mab_collect_depth', # '600_mab_collect_depth', '50_mab_mass_flux_mg_m-2_d-1', '50_mab_POC_flux_mg_C_m-2_d-1', '50_mab_Total_Nitrogen_flux_mg_m-2_d-1', # '50_mab_Total_Carbon_flux_mg_C_m-2_d-1', '600_mab_mass_flux_mg_m-2_d-1', '600_mab_POC_flux_mg_C_m-2_d-1', # '600_mab_Total_Nitrogen_flux_mg_m-2_d-1', '600_mab_Total_Carbon_flux_mg_C_m-2_d-1', 'GAPFILL600_Massflux_mg_m-2_d-1', # 'GAPFILL600_POCflux_mg_C_m-2_d-1', 'GAPFILL600_Total_Nitrogen_flux_mg_m-2_d-1', 'GAPFILL600_Total_Carbon_flux_mg_C_m-2_d-1', 'Notes'] POC_STation_M_SES_and_McLane_subset = POC_STation_M_SES_and_McLane[['POCmodel_mg_daily_mean', 'PULSE', 'Cup_number', 'GAPFILL600_Massflux_mg_m-2_d-1','GAPFILL600_POCflux_mg_C_m-2_d-1']] full_idx = pd.date_range(start=POC_STation_M_SES_and_McLane_subset.index.min(), end=POC_STation_M_SES_and_McLane_subset.index.max(), freq='D') # Reindex safely POC_STation_M_SES_and_McLane_subset_padded = POC_STation_M_SES_and_McLane_subset.reindex(full_idx) # Reindex the DataFrame to include missing dates, leaving values blank (NaN) Collect_times = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/SES collect times.csv') POC_STation_M_SES_and_McLane_subset_padded['Collect_HOURS'] = POC_STation_M_SES_and_McLane_subset_padded['PULSE'].map( Collect_times.set_index('PULSE')['Collect_HOURS'] ) POC_STation_M_SES_and_McLane_subset_padded["Pulse_Cup"] = POC_STation_M_SES_and_McLane_subset_padded["PULSE"].astype(str) + "_" + POC_STation_M_SES_and_McLane_subset_padded["Cup_number"].astype(str) POC_STation_M_SES_and_McLane_subset_padded.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/POC_STation_M_SES_and_McLane_subset_padded.csv') # Plot using df.index as the first argument plt.plot(POC_STation_M_SES_and_McLane_subset_padded.index, POC_STation_M_SES_and_McLane_subset_padded["POCmodel_mg_daily_mean"], label="SES", color="blue") plt.plot(POC_STation_M_SES_and_McLane_subset_padded.index, POC_STation_M_SES_and_McLane_subset_padded["GAPFILL600_POCflux_mg_C_m-2_d-1"], label="Trap", color="red", linestyle="--") plt.ylabel('POC flux mg C m-2 d-1') plt.title('Station M') plt.legend() plt.savefig("W:/Analyses/Station M data/Station M attenuance_Not red channel/SES_TrapComparisonStationM.pdf", bbox_inches="tight") plt.show() #%% #START HERE POC_STation_M_SES_and_McLane_subset_padded = pd.read_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/POC_STation_M_SES_and_McLane_subset_padded.csv') #%% ##Compare cups and SES directly POC_STation_M_SES_and_McLaneComplete = POC_STation_M_SES_and_McLane_subset_padded.dropna() #average by cup ID POC_STation_M_SES_and_McLaneCupAvg = POC_STation_M_SES_and_McLaneComplete.groupby('Pulse_Cup').mean(numeric_only=True) column_list = POC_STation_M_SES_and_McLaneCupAvg.columns.tolist() print(column_list) #['POCmodel_mg_daily_mean', 'PULSE', 'Cup_number', 'GAPFILL600_Massflux_mg_m-2_d-1', 'GAPFILL600_POCflux_mg_C_m-2_d-1', 'Collect_HOURS'] #Find diff between Trap and SES POC as 5 of trap POC POC_STation_M_SES_and_McLaneCupAvg['perc'] = ((POC_STation_M_SES_and_McLaneCupAvg['POCmodel_mg_daily_mean']/POC_STation_M_SES_and_McLaneCupAvg['GAPFILL600_POCflux_mg_C_m-2_d-1']))*100 POC_STation_M_SES_and_McLaneCupAvg['diff_abs'] = (POC_STation_M_SES_and_McLaneCupAvg['POCmodel_mg_daily_mean']-POC_STation_M_SES_and_McLaneCupAvg['GAPFILL600_POCflux_mg_C_m-2_d-1']) #box plot and jitter plot: x = collect hour, boxes are diff_perc and diff_abs colored by pulse, each dot is a cup POC_STation_M_SES_and_McLaneCupAvg.to_csv('W:/Analyses/Station M data/Station M attenuance_Not red channel/POC_STation_M_SES_and_McLaneCupAvg.csv') # Plot using 'hue' to automatically group by color POC_STation_M_SES_and_McLaneCupAvg.plot.scatter( x='Cup_number', y='diff_abs', s='Collect_HOURS', c='PULSE') plt.title("Scatter Plot Colored by Group ID") plt.show() sns.boxplot(data=POC_STation_M_SES_and_McLaneCupAvg, x="Cup_number", y="perc") plt.show() #Would need to do more complicated stats but initial explorations here and excel don't show strong cup number effect, though for some deployments it could be an issue GZ?) #Trap vs SES, line each with equation/r^2 for each collect time, colored and sized by cup number POC_STation_M_SES_and_McLaneCupAvg.plot.scatter( x='GAPFILL600_POCflux_mg_C_m-2_d-1', y='POCmodel_mg_daily_mean', s='Collect_HOURS', c='PULSE') x=POC_STation_M_SES_and_McLaneCupAvg['GAPFILL600_POCflux_mg_C_m-2_d-1'] y1=POC_STation_M_SES_and_McLaneCupAvg['POCmodel_mg_daily_mean'] y2=POC_STation_M_SES_and_McLaneCupAvg['POCmodel_mg_daily_mean'] y3=POC_STation_M_SES_and_McLaneCupAvg['POCmodel_mg_daily_mean'] # 2. Calculate Equations and R² for each series without loops # Series 1 (y1) m1, b1 = np.polyfit(POC_STation_M_SES_and_McLaneCupAvg['x'], POC_STation_M_SES_and_McLaneCupAvg['y1'], 1) r1_sq = np.corrcoef(POC_STation_M_SES_and_McLaneCupAvg['x'], POC_STation_M_SES_and_McLaneCupAvg['y1'])[0, 1] ** 2 # Series 2 (y2) m2, b2 = np.polyfit(POC_STation_M_SES_and_McLaneCupAvg['x'], POC_STation_M_SES_and_McLaneCupAvg['y2'], 1) r2_sq = np.corrcoef(POC_STation_M_SES_and_McLaneCupAvg['x'], POC_STation_M_SES_and_McLaneCupAvg['y2'])[0, 1] ** 2 # Series 3 (y3) m3, b3 = np.polyfit(POC_STation_M_SES_and_McLaneCupAvg['x'], POC_STation_M_SES_and_McLaneCupAvg['y3'], 1) r3_sq = np.corrcoef(POC_STation_M_SES_and_McLaneCupAvg['x'], POC_STation_M_SES_and_McLaneCupAvg['y3'])[0, 1] ** 2 # 3. Initialize Plot plt.figure(figsize=(11, 7)) # 4. Explicitly plot each scatter and trendline # Plot Series 1 plt.scatter(POC_STation_M_SES_and_McLaneCupAvg['x'], POC_STation_M_SES_and_McLaneCupAvg['y1'], s=POC_STation_M_SES_and_McLaneCupAvg['size_val'], c=POC_STation_M_SES_and_McLaneCupAvg['color_val'], cmap='viridis', marker='o', alpha=0.7, label=f"y1: y = {m1:.2f}x + {b1:.2f} ($R^2$={r1_sq:.3f})") plt.plot(POC_STation_M_SES_and_McLaneCupAvg['x'], m1 * POC_STation_M_SES_and_McLaneCupAvg['x'] + b1, linestyle='--', color='teal', alpha=0.8) # Plot Series 2 plt.scatter(POC_STation_M_SES_and_McLaneCupAvg['x'], POC_STation_M_SES_and_McLaneCupAvg['y2'], s=POC_STation_M_SES_and_McLaneCupAvg['size_val'], c=POC_STation_M_SES_and_McLaneCupAvg['color_val'], cmap='plasma', marker='s', alpha=0.7, label=f"y2: y = {m2:.2f}x + {b2:.2f} ($R^2$={r2_sq:.3f})") plt.plot(POC_STation_M_SES_and_McLaneCupAvg['x'], m2 * POC_STation_M_SES_and_McLaneCupAvg['x'] + b2, linestyle='--', color='magenta', alpha=0.8) # Plot Series 3 c_scat = plt.scatter(POC_STation_M_SES_and_McLaneCupAvg['x'], POC_STation_M_SES_and_McLaneCupAvg['y3'], s=POC_STation_M_SES_and_McLaneCupAvg['size_val'], c=POC_STation_M_SES_and_McLaneCupAvg['color_val'], cmap='cool', marker='^', alpha=0.7, label=f"y3: y = {m3:.2f}x + {b3:.2f} ($R^2$={r3_sq:.3f})") plt.plot(POC_STation_M_SES_and_McLaneCupAvg['x'], m3 * df['x'] + b3, linestyle='--', color='blue', alpha=0.8) # 5. Visual Formatting plt.title("Multi-Series Scatter Plot (Loop-Free implementation)", fontsize=14, pad=15) plt.xlabel("X Axis Variable", fontsize=12) plt.ylabel("Y Axis Variables", fontsize=12) plt.grid(True, linestyle=':', alpha=0.6) # Map the colorbar based on the final layer mapping cbar = plt.colorbar(c_scat) cbar.set_label('Color Gradient Scale (color_val)', fontsize=10) # Render legend outside the main plot area plt.legend(loc='upper left', bbox_to_anchor=(1.2, 1), title="Regression Equations") plt.tight_layout() plt.show()