# -*- coding: utf-8 -*- """ Created on Thu Oct 8 11:37:10 2026 @author: chuffard """ import pandas as pd import matplotlib.pyplot as plt import seaborn as sns import numpy as np # Load the JSON file into a DataFrame SINKER = pd.read_json('W:/Analyses/SINKER comparison/flux_timeseries.json', lines=True) #Convert SINKER to local time SINKER['date_timeUTC']= pd.to_datetime(SINKER['date_time']) # Convert to Pacific Time SINKER['DateTime_Local'] = SINKER['date_timeUTC'].dt.tz_convert('US/Pacific') #TRIM these bad data periods # Start_time End_time Reason # 9/16/2026 18:00 9/18/2026 19:00 Recovery and redeploy period # 3/4/2026 0:00 3/25/2026 0:00 wiper misaligned # 7/15/2025 0:00 7/18/2025 3:00 not opening on purpose # 7/21/2025 15:00 7/22/2025 19:00 accuator jammed # 12/8/2025 12/9/2025 Bright image with poor detection at 23:00 cut_outs = [ {'start': '9/16/2026 18:00', 'end': '9/18/2026 19:00'}, {'start': '3/4/2026 0:00', 'end': '3/25/2026 0:00'}, {'start': '7/15/2025 0:00', 'end': '7/18/2025 3:00'}, {'start': '7/21/2025 15:00', 'end': '7/22/2025 19:00'}, {'start': '12/8/2025 00:00', 'end': '12/9/2025 0:00'} ] # 3. Build a combined condition dynamically # Start with all True (keep everything) keep_mask = pd.Series(True, index=SINKER.index) for period in cut_outs: # Identify rows inside this specific blackout window in_blackout = (SINKER['date_timeUTC'] >= pd.to_datetime(period['start'], utc=True)) & \ (SINKER['date_timeUTC'] <= pd.to_datetime(period['end'], utc=True)) # Update the mask: if it's in the blackout, we don't keep it keep_mask = keep_mask & ~in_blackout # 4. Filter the dataframe filtered_SINKER = SINKER[keep_mask] filtered_SINKER = filtered_SINKER.drop(columns=['date_timeUTC', 'date_time']) hourly_SINKER =filtered_SINKER.groupby( filtered_SINKER['DateTime_Local'].dt.tz_localize(None).dt.floor('h') ).mean(numeric_only=True) #Import SES 2025 MARS MARS2025SES= pd.read_csv("W:/Analyses/2025 MARS deployment/MARS2025_atn_POC_red_5only.csv",parse_dates=True) MARS2025SES['DateTime_Local']= pd.to_datetime(MARS2025SES['DateTime']) MARS2025SES['Location'] = "MARS" #Import SES 2026 Mooring Mooring2026_SES = pd.read_csv("W:/Analyses/2026 Mooring/Mooring2026atn_POC_red.csv", parse_dates=True) Mooring2026_SES['DateTime_Local']= pd.to_datetime(Mooring2026_SES['DateTime']) Mooring2026_SES['Location'] = "OffShore" SES_Merge = pd.concat([Mooring2026_SES, MARS2025SES], ignore_index=True) SES_Merge['date'] = SES_Merge['DateTime_Local'].dt.date hourly_SES =SES_Merge.groupby( SES_Merge['DateTime_Local'].dt.tz_localize(None).dt.floor('h') ).mean(numeric_only=True) hourly_SES['DateTime_Local_df'] = hourly_SES.index hourly_SES["date"] = hourly_SES.index.date hourly_SES = hourly_SES.merge(SES_Merge[['date', 'Location']], on='date', how='left') # Reset index of df1, merge on 'year', and set the original index back hourly_SES = hourly_SES.reset_index().merge(SES_Merge[['date']], on='date').set_index('DateTime_Local_df') hourly_SES.rename(columns={'Location': 'SES_Location'}, inplace=True) hourly_SES.index.name = "DateTime_Local" hourly_SES['DateTime_Local_df'] = hourly_SES.index #CHECK BEFORE DOING hourly_SES = hourly_SES.iloc[:, 3:] hourly_SES = hourly_SES.drop_duplicates() Hourly_all = pd.merge(hourly_SES, hourly_SINKER, left_index=True, right_index=True, how='outer') # 2. Extract the date (or floor to day) to group by day Hourly_all["date"] = Hourly_all.index.date # 3. Use transform('mean') to calculate daily averages broadcasted to every hourly row #daily_means = Hourly_all.groupby("date").transform("mean") Hourly_all.columns # Index(['atn', 'time', 'F_atn', 'POCmodel_mg', 'POCmodel_g', 'date', # 'SES_Location', 'DateTime_Local_df', 'PSD_slope', 'total_area_flux'], # dtype='object') target_cols = ['atn', 'F_atn', 'POCmodel_mg','total_area_flux','PSD_slope'] daily_means = Hourly_all.groupby("date")[target_cols].transform("mean") # 4. Add them as new columns (e.g., with a suffix like '_daily_avg') Hourly_all[[col + "_daily_avg" for col in daily_means.columns]] = daily_means #PAD date time index start_date = min(Hourly_all['date']) end_date = max(Hourly_all['date']) # Create a complete daily date range full_date_range = pd.date_range(start=start_date, end=end_date, freq='h') #Hourly_all= Hourly_all.set_index('date', inplace=False)#has columns like atn Hourly_all = Hourly_all.reindex(full_date_range) #Hourly_all.to_csv("Z:/STORE/Files/Manuscripts/Pyrosomes/2025/SCOC_daily_pad.csv") #Hourly_all = Hourly_all.reset_index() Hourly_all['DateTime_Local_df'] = Hourly_all.index Hourly_all.to_csv("W:/Analyses/SINKER comparison/Hourly_all.csv") Hourly_all.columns # Index(['atn', 'time', 'F_atn', 'POCmodel_mg', 'POCmodel_g', 'date', # 'SES_Location', 'DateTime_Local_df', 'PSD_slope', 'total_area_flux', # 'atn_daily_avg', 'F_atn_daily_avg', 'POCmodel_mg_daily_avg', # 'total_area_flux_daily_avg', 'PSD_slope_daily_avg'], # dtype='object') #%% ##Now add triggering import numpy as np from datetime import datetime import numpy as np import pandas as pd import sys import os.path import warnings import glob import seaborn as sns import matplotlib.pyplot as plt from matplotlib.backends.backend_pdf import PdfPages warnings.filterwarnings("ignore") atnSES = Hourly_all[['DateTime_Local_df','atn','SES_Location', 'atn_daily_avg']] areafluxSINKER= Hourly_all[['DateTime_Local_df','total_area_flux', 'total_area_flux_daily_avg']] #Trigger SES one channel atnSES["DateTime_Local_df"] = pd.to_datetime(atnSES['DateTime_Local_df']) atnSES = atnSES.sort_values(by='DateTime_Local_df') atnSES['atn'] = pd.to_numeric(atnSES['atn'], errors='coerce') atnSES['rm150atn'] = atnSES['atn'].rolling(150, min_periods=23).mean() atnSES['Slope_rm_150atn'] = atnSES['rm150atn'].diff() #triggering values lowtrigger_atn = 0.007 hightrigger_atn = 0.011 hightrigger_fluoro= 2.4 lowtrigger_fluoro= 1.8 use_min_period = 24 nrec=150 # Trigger code #atn trigger base atnSES["TriggerHigh_atn"] = np.where(atnSES['Slope_rm_150atn'] >= hightrigger_atn, 'Trigger', 'Wait') atnSES['Rate_HighTrigger_atn'] = (atnSES['TriggerHigh_atn'] == 'Trigger').rolling(144, min_periods=2).sum() atnSES['Max_HighTrigger_rate_atn'] = atnSES['Rate_HighTrigger_atn'].cummax() atnSES["TriggerLow_atn"] = np.where(atnSES['Slope_rm_150atn'] >= lowtrigger_atn, 'Trigger', 'Wait') atnSES['Rate_LowTrigger_atn'] = (atnSES['TriggerLow_atn'] == 'Trigger').rolling(144, min_periods=2).sum() atnSES['Max_LowTrigger_rate_atn'] = atnSES['Rate_LowTrigger_atn'].cummax() conditionsatn = [ (atnSES['Rate_LowTrigger_atn'] <= 1), #fire one low, then switch to high trigger (atnSES['Rate_HighTrigger_atn'] <= 1), #fire two high, then switch to higher trigger (atnSES['Rate_HighTrigger_atn'] ==3), #fire more 3 highs then even higher trigger (atnSES['Max_HighTrigger_rate_atn'] >= 4),#fire four or more highs, then use max trigger ((atnSES['Rate_HighTrigger_atn']).rolling(window=240, min_periods=1).sum()==0), #more than 1.5 week since last trigger, use the lower gear ] # Define corresponding choices (results) choicesatn = [np.where(atnSES['Slope_rm_150atn'] >= lowtrigger_atn, 'FireLow', 'Wait'), np.where(atnSES['Slope_rm_150atn'] >= hightrigger_atn, 'FireHigh', 'Wait'), np.where(atnSES['Slope_rm_150atn'] >= (hightrigger_atn*1.1), 'Fire1_1High', 'Wait'), np.where(atnSES['Slope_rm_150atn'] >= (hightrigger_atn*1.2), 'Fire1_2High', 'Wait'), np.where(atnSES['Slope_rm_150atn'] >= (hightrigger_atn), 'Fire_High', 'Wait')] atnSES["Fire_trigger_atn"] = np.select(conditionsatn, choicesatn, default='Wait') #Send trigger signal on "Get_down == "Fire!" atnSES["Get_Down"] = np.where( (atnSES['Fire_trigger_atn'] != 'Wait').rolling(window=48, min_periods=1).sum()==1, np.where((atnSES['Fire_trigger_atn'] != 'Wait'), 'Fire!', 'Wait'), np.where(((atnSES['Fire_trigger_atn'] == 'Wait').rolling(336, min_periods=2).sum()==0), 'Fire!', 'Wait'), ) #save results atnSES.to_csv('W:/Analyses/SINKER comparison/SES_triggerATN.csv') #%% SINKER triggering #Trigger SES one channel areafluxSINKER["DateTime_Local_df"] = pd.to_datetime(areafluxSINKER['DateTime_Local_df']) areafluxSINKER = areafluxSINKER.sort_values(by='DateTime_Local_df') areafluxSINKER['total_area_flux'] = pd.to_numeric(areafluxSINKER['total_area_flux'], errors='coerce') areafluxSINKER['rm150total_area_flux'] = areafluxSINKER['total_area_flux'].rolling(150, min_periods=23).mean() areafluxSINKER['Slope_rm_150total_area_flux'] = areafluxSINKER['rm150total_area_flux'].diff() #triggering values lowtrigger_total_area_flux = 0.007 hightrigger_total_area_flux = 0.011 hightrigger_fluoro= 2.4 lowtrigger_fluoro= 1.8 use_min_period = 24 nrec=150 # Trigger code #total_area_flux trigger base areafluxSINKER["TriggerHigh_total_area_flux"] = np.where(areafluxSINKER['Slope_rm_150total_area_flux'] >= hightrigger_total_area_flux, 'Trigger', 'Wait') areafluxSINKER['Rate_HighTrigger_total_area_flux'] = (areafluxSINKER['TriggerHigh_total_area_flux'] == 'Trigger').rolling(144, min_periods=2).sum() areafluxSINKER['Max_HighTrigger_rate_total_area_flux'] = areafluxSINKER['Rate_HighTrigger_total_area_flux'].cummax() areafluxSINKER["TriggerLow_total_area_flux"] = np.where(areafluxSINKER['Slope_rm_150total_area_flux'] >= lowtrigger_total_area_flux, 'Trigger', 'Wait') areafluxSINKER['Rate_LowTrigger_total_area_flux'] = (areafluxSINKER['TriggerLow_total_area_flux'] == 'Trigger').rolling(144, min_periods=2).sum() areafluxSINKER['Max_LowTrigger_rate_total_area_flux'] = areafluxSINKER['Rate_LowTrigger_total_area_flux'].cummax() conditionstotal_area_flux = [ (areafluxSINKER['Rate_LowTrigger_total_area_flux'] <= 1), #fire one low, then switch to high trigger (areafluxSINKER['Rate_HighTrigger_total_area_flux'] <= 1), #fire two high, then switch to higher trigger (areafluxSINKER['Rate_HighTrigger_total_area_flux'] ==3), #fire more 3 highs then even higher trigger (areafluxSINKER['Max_HighTrigger_rate_total_area_flux'] >= 4),#fire four or more highs, then use max trigger ((areafluxSINKER['Rate_HighTrigger_total_area_flux']).rolling(window=240, min_periods=1).sum()==0), #more than 1.5 week since last trigger, use the lower gear ] # Define corresponding choices (results) choicestotal_area_flux = [np.where(areafluxSINKER['Slope_rm_150total_area_flux'] >= lowtrigger_total_area_flux, 'FireLow', 'Wait'), np.where(areafluxSINKER['Slope_rm_150total_area_flux'] >= hightrigger_total_area_flux, 'FireHigh', 'Wait'), np.where(areafluxSINKER['Slope_rm_150total_area_flux'] >= (hightrigger_total_area_flux*1.1), 'Fire1_1High', 'Wait'), np.where(areafluxSINKER['Slope_rm_150total_area_flux'] >= (hightrigger_total_area_flux*1.2), 'Fire1_2High', 'Wait'), np.where(areafluxSINKER['Slope_rm_150total_area_flux'] >= (hightrigger_total_area_flux), 'Fire_High', 'Wait')] areafluxSINKER["Fire_trigger_total_area_flux"] = np.select(conditionstotal_area_flux, choicestotal_area_flux, default='Wait') #Send trigger signal on "Get_down == "Fire!" areafluxSINKER["Get_Down"] = np.where( (areafluxSINKER['Fire_trigger_total_area_flux'] != 'Wait').rolling(window=48, min_periods=1).sum()==1, np.where((areafluxSINKER['Fire_trigger_total_area_flux'] != 'Wait'), 'Fire!', 'Wait'), np.where(((areafluxSINKER['Fire_trigger_total_area_flux'] == 'Wait').rolling(336, min_periods=2).sum()==0), 'Fire!', 'Wait'), ) #save results areafluxSINKER.to_csv('W:/Analyses/SINKER comparison/SINKER_Trigger_areaflux.csv') #%% #Plot RM 150 triggers custom_markers = {"MARS": "o", "OffShore": "+"} palette =["#3776ab", "#FF8C00"] custom_sizes = {"Wait": 50, "Fire!": 400} # 3. Create the first scatterplot (Top) fig, axes = plt.subplots(2, 1, figsize=(10, 10)) # 2. Plot on the first axis (axes[0]) sns.scatterplot(data=atnSES, x='DateTime_Local_df', y='rm150atn', hue='Get_Down',size="Get_Down", style="SES_Location", sizes=custom_sizes, ax=axes[0], alpha = 0.3) # 3. Plot on the second axis (axes[1]) sns.scatterplot(data=areafluxSINKER, x='DateTime_Local_df', y='rm150total_area_flux', hue='Get_Down',size="Get_Down", sizes=custom_sizes, ax=axes[1], alpha = 0.3) # Adjust layout to prevent overlap plt.tight_layout() plt.show() fig.savefig('W:/Analyses/SINKER comparison/SES_Sinker_triggers.pdf', dpi=300, bbox_inches='tight') #%% Plot data with triggers #Plot custom_markers = {"MARS": "o", "OffShore": "+"} palette =["#3776ab", "#FF8C00"] custom_sizes = {"Wait": 50, "Fire!": 400} # 3. Create the first scatterplot (Top) fig, axes = plt.subplots(2, 1, figsize=(10, 10)) # 2. Plot on the first axis (axes[0]) sns.scatterplot(data=atnSES, x='DateTime_Local_df', y='atn', hue='Get_Down',size="Get_Down", style="SES_Location", sizes=custom_sizes, ax=axes[0], alpha = 0.2) sns.lineplot(data=atnSES, x='DateTime_Local_df', y='atn_daily_avg', ax=axes[0], color = 'k', style="SES_Location") # 3. Plot on the second axis (axes[1]) sns.scatterplot(data=areafluxSINKER, x='DateTime_Local_df', y='total_area_flux', hue='Get_Down',size="Get_Down", sizes=custom_sizes, ax=axes[1], alpha = 0.2) sns.lineplot(data=Hourly_all, x='DateTime_Local_df', y='total_area_flux_daily_avg', ax=axes[1], color = 'k') # Adjust layout to prevent overlap plt.tight_layout() plt.show() fig.savefig('W:/Analyses/SINKER comparison/SES_Sinker_raw.pdf', dpi=300, bbox_inches='tight') #%% #plot SES vs sinker when offshore and MARS Hourly_all['date'] = pd.to_datetime(Hourly_all['date']) # 3. Filter between two dates SES_MARS_Start = Hourly_all.index[Hourly_all['SES_Location'] == 'MARS'].min() SES_MARS_end = Hourly_all.index[Hourly_all['SES_Location'] == 'MARS'].max() SES_Offshore_start = Hourly_all.index[Hourly_all['SES_Location'] == 'OffShore'].min() SES_Offshore_end = Hourly_all.index[Hourly_all['SES_Location'] == 'OffShore'].max() SES_MARS_SINKER = Hourly_all.loc[SES_MARS_Start:SES_MARS_end].copy() SES_OffShore_SINKER= Hourly_all.loc[SES_Offshore_start:SES_Offshore_end].copy() SES_OffShore_SINKER.to_csv('W:/Analyses/SINKER comparison/SES_OffShore_SINKER.csv') SES_MARS_SINKER.to_csv('W:/Analyses/SINKER comparison/SES_MARS_SINKER.csv') SES_OffShore_SINKER.columns # Index(['atn', 'time', 'F_atn', 'POCmodel_mg', 'POCmodel_g', 'date', # 'SES_Location', 'DateTime_Local_df', 'PSD_slope', 'total_area_flux', # 'atn_daily_avg', 'F_atn_daily_avg', 'POCmodel_mg_daily_avg', # 'total_area_flux_daily_avg', 'PSD_slope_daily_avg'], # dtype='object') # 3. Create the first scatterplot (Top) fig, axes = plt.subplots(1, 2, figsize=(10, 5)) # 2. Plot on the first axis (axes[0]) sns.scatterplot(data=SES_OffShore_SINKER, x='atn_daily_avg', y='total_area_flux_daily_avg', ax=axes[0], alpha = 0.2) axes[0].set_title("Offshore daily comparison") # 3. Plot on the second axis (axes[1]) sns.scatterplot(data=SES_MARS_SINKER, x='atn_daily_avg', y='total_area_flux_daily_avg', ax=axes[1], alpha = 0.2) axes[1].set_title("MARS daily comparison") # Adjust layout to prevent overlap plt.tight_layout() plt.show() fig.savefig('W:/Analyses/SINKER comparison/SES_Sinker_daily.pdf', dpi=300, bbox_inches='tight') #%% Cross correlation sliding window lags # #Negative lag, positive correlations # import numpy as np # import pandas as pd # import matplotlib.pyplot as plt # from matplotlib.colors import Normalize # def sliding_ccf(ts_lead, ts_lag, win_width=60, # lag_max=30, plot_results=True, dates=None): # """ # Sliding-window cross-correlation. # Parameters # ---------- # ts_lead : array-like # Leading variable (e.g., POCmodel_mg). # ts_lag : array-like # Lagging variable (e.g., total_area_flux). # win_width : int # Number of records in each window. # lag_max : int # Maximum absolute lag to examine. # plot_results : bool # Whether to plot the lag and correlation profiles. # dates : array-like, optional # Dates for the x-axis. # Returns # ------- # DataFrame with columns: # MaxNegLag # MaxNegCorrelation # """ # # Convert inputs to numeric arrays # lead = np.asarray(ts_lead, dtype=float) # lag = np.asarray(ts_lag, dtype=float) # # Check inputs # if len(lead) != len(lag): # raise ValueError( # "Lead and lag time series must have equal lengths." # ) # n = len(lead) # if win_width < 2: # raise ValueError("win_width must be at least 2.") # if n < 2: # raise ValueError("At least two observations are required.") # if win_width > n: # print("Window exceeds data length; using full data length.") # win_width = n # if lag_max >= win_width: # print("lag_max exceeds window width; reducing it.") # lag_max = win_width - 1 # if lag_max < 0: # raise ValueError("lag_max cannot be negative.") # half = win_width // 2 # # Match the original R window alignment. # # For even windows, the center falls between two records. # start_pos = half - 1 if win_width % 2 == 0 else half # # Initialize outputs, retaining one result per original record. # max_corr = np.full(n, np.nan) # max_lag = np.full(n, np.nan) # # Calculate correlation at a particular lag. # # This convention matches R's ccf(x, y) lag definition: # # lag k compares lead[t+k] with lag[t]. # def lag_correlation(x, y, k): # if k < 0: # x_pair = x[:k] # y_pair = y[-k:] # elif k > 0: # x_pair = x[k:] # y_pair = y[:-k] # else: # x_pair = x # y_pair = y # # Pairwise removal of missing/non-finite observations # valid = np.isfinite(x_pair) & np.isfinite(y_pair) # x_pair = x_pair[valid] # y_pair = y_pair[valid] # if len(x_pair) < 2: # return np.nan # # A correlation is undefined for a constant series. # if np.std(x_pair) == 0 or np.std(y_pair) == 0: # return np.nan # return np.corrcoef(x_pair, y_pair)[0, 1] # # Only evaluate negative lags and zero, as in the R script. # lags_to_test = np.arange(-lag_max, 1) # # Sliding-window loop # for i in range(start_pos, n - half): # if win_width % 2 == 0: # window_start = i - half + 1 # else: # window_start = i - half # window_end = window_start + win_width # x = lead[window_start:window_end] # y = lag[window_start:window_end] # correlations = np.array([ # lag_correlation(x, y, k) # for k in lags_to_test # ]) # valid = np.isfinite(correlations) # if not valid.any(): # continue # # Select the largest ABSOLUTE correlation, # # but retain its original sign in the output. # valid_indices = np.flatnonzero(valid) # best_idx = valid_indices[ # np.argmax(np.abs(correlations[valid])) # ] # max_corr[i] = correlations[best_idx] # max_lag[i] = lags_to_test[best_idx] # # Return one row per original observation. # results = pd.DataFrame({ # "MaxNegLag": max_lag, # "MaxNegCorrelation": max_corr # }) # if plot_results: # x_axis = ( # np.asarray(dates) # if dates is not None # else np.arange(n) # ) # valid_results = ( # np.isfinite(max_lag) & np.isfinite(max_corr) # ) # # Only color points where the selected correlation is positive # positive_corr = valid_results & (max_corr > 0) # # Plot 1: profile colored by lag (unchanged) # fig, ax = plt.subplots(figsize=(13, 5)) # ax.plot( # x_axis, lag, # color="lightgray", # linewidth=1.5, # zorder=1 # ) # sc = ax.scatter( # np.asarray(x_axis)[positive_corr], # lag[positive_corr], # c=max_lag[positive_corr], # cmap="hot_r", # norm=Normalize(vmin=-lag_max, vmax=0), # s=18, # zorder=2 # ) # ax.set_title( # f"Lagging Variable Profile Colored by Lag, for positive correlations " # f"(window = {win_width} records)" # ) # ax.set_ylabel("SINKER area flux (Lagging variable)") # ax.set_xlabel("Date" if dates is not None else "Record") # fig.colorbar(sc, ax=ax, label="Selected lag (records)") # fig.tight_layout() # plt.show() # # Plot 2: profile colored by POSITIVE correlation only # fig, ax = plt.subplots(figsize=(13, 5)) # ax.plot( # x_axis, lag, # color="lightgray", # linewidth=1.5, # zorder=1 # ) # sc = ax.scatter( # np.asarray(x_axis)[positive_corr], # lag[positive_corr], # c=max_corr[positive_corr], # cmap="coolwarm", # norm=Normalize(vmin=0, vmax=1), # s=18, # zorder=2 # ) # ax.set_title( # f"Lagging Variable Profile Colored by Positive Correlation " # f"(window = {win_width} records)" # ) # ax.set_ylabel("SINKER area flux (Lagging variable)") # ax.set_xlabel("Date" if dates is not None else "Record") # if positive_corr.any(): # fig.colorbar( # sc, ax=ax, # label="Positive correlation" # ) # fig.tight_layout() # plt.show() # # Plot 3: Leading and lagging variables with two y-axes # fig, ax1 = plt.subplots(figsize=(13, 5)) # # Left y-axis: leading variable # line1 = ax1.plot( # x_axis, # lead, # color="tab:blue", # linewidth=1.5, # label="SES modeled POC flux (leading variable)" # ) # ax1.set_xlabel("Date" if dates is not None else "Record") # ax1.set_ylabel("SES modeled POC flux (leading variable)", color="tab:blue") # ax1.tick_params(axis="y", labelcolor="tab:blue") # # Right y-axis: lagging variable # ax2 = ax1.twinx() # line2 = ax2.plot( # x_axis, # lag, # color="tab:red", # linewidth=1.5, # label="SINKER area flux (Lagging variable)" # ) # ax2.set_ylabel("SINKER area flux (Lagging variable)", color="tab:red") # ax2.tick_params(axis="y", labelcolor="tab:red") # # Combined legend # lines = line1 + line2 # labels = [line.get_label() for line in lines] # ax1.legend(lines, labels, loc="best") # ax1.set_title( # "Leading and Lagging Variables Through Time" # ) # fig.tight_layout() # plt.show() # return results # #now run it # slidingdata = SES_OffShore_SINKER[['date', 'total_area_flux', 'POCmodel_mg_daily_avg']] # #slidingdatadaily = slidingdata.groupby("date").transform("mean") # slidingdatadaily = slidingdata.resample('D').mean() # myData = slidingdatadaily # # Use your actual date column name if one exists. # # Replace "date" below with its name. # myData["date"] = pd.to_datetime(myData["date"]) # mySlidingCCF = sliding_ccf( # ts_lead=myData["POCmodel_mg_daily_avg"], # ts_lag=myData["total_area_flux"], # win_width=60, # lag_max=30, # plot_results=True, # dates=myData["date"] # ) #%%All lags, positive correlations import numpy as np import pandas as pd import matplotlib.pyplot as plt from matplotlib.colors import Normalize def sliding_ccf(ts_lead, ts_lag, win_width=60, lag_max=30, plot_results=True, dates=None): """ Sliding-window cross-correlation for positive and negative lags. Returns ------- DataFrame with columns: MaxLag MaxCorrelation Lag convention: k > 0: compares lead[t+k] with lag[t] k < 0: compares lead[t] with lag[t+abs(k)] k = 0: compares lead[t] with lag[t] """ lead = np.asarray(ts_lead, dtype=float) lag = np.asarray(ts_lag, dtype=float) if len(lead) != len(lag): raise ValueError( "Lead and lag time series must have equal lengths." ) n = len(lead) if win_width < 2: raise ValueError("win_width must be at least 2.") if n < 2: raise ValueError("At least two observations are required.") if win_width > n: print("Window exceeds data length; using full data length.") win_width = n if lag_max < 0: raise ValueError("lag_max cannot be negative.") if lag_max >= win_width: print("lag_max exceeds window width; reducing it.") lag_max = win_width - 1 half = win_width // 2 # Preserve original window alignment. start_pos = half - 1 if win_width % 2 == 0 else half max_corr = np.full(n, np.nan) max_lag = np.full(n, np.nan) def lag_correlation(x, y, k): # Positive lag: lead occurs earlier than lag if k > 0: x_pair = x[:-k] y_pair = y[k:] # Negative lag: lead occurs later than lag elif k < 0: x_pair = x[-k:] y_pair = y[:k] else: x_pair = x y_pair = y # Pairwise removal of missing/non-finite values valid = np.isfinite(x_pair) & np.isfinite(y_pair) x_pair = x_pair[valid] y_pair = y_pair[valid] if len(x_pair) < 2: return np.nan if np.std(x_pair) == 0 or np.std(y_pair) == 0: return np.nan return np.corrcoef(x_pair, y_pair)[0, 1] # Test every lag from -lag_max through +lag_max lags_to_test = np.arange(-lag_max, lag_max + 1) # Sliding-window loop for i in range(start_pos, n - half): if win_width % 2 == 0: window_start = i - half + 1 else: window_start = i - half window_end = window_start + win_width x = lead[window_start:window_end] y = lag[window_start:window_end] correlations = np.array([ lag_correlation(x, y, k) for k in lags_to_test ]) valid = np.isfinite(correlations) if not valid.any(): continue # Select strongest absolute correlation, # retaining its original sign. valid_indices = np.flatnonzero(valid) best_idx = valid_indices[ np.argmax(np.abs(correlations[valid])) ] max_corr[i] = correlations[best_idx] max_lag[i] = lags_to_test[best_idx] results = pd.DataFrame({ "MaxLag": max_lag, "MaxCorrelation": max_corr }) if plot_results: x_axis = ( np.asarray(dates) if dates is not None else np.arange(n) ) valid_results = ( np.isfinite(max_lag) & np.isfinite(max_corr) ) # Only color points with positive selected correlations positive_corr = valid_results & (max_corr > 0) # Plot 1: lagging-variable profile colored by selected lag fig, ax = plt.subplots(figsize=(13, 5)) ax.plot( x_axis, lag, color="lightgray", linewidth=1.5, zorder=1 ) sc = ax.scatter( x_axis[positive_corr], lag[positive_corr], c=max_lag[positive_corr], cmap="coolwarm", norm=Normalize(vmin=-lag_max, vmax=lag_max), s=70, zorder=2 ) ax.set_title( f"Positive means SES leads SINKER. Negative means SINKER leads SES" f"(positive correlations only; window = {win_width})" ) ax.set_ylabel("SINKER area flux") ax.set_xlabel("Date" if dates is not None else "Record") fig.colorbar(sc, ax=ax, label="Selected lag (records)") fig.tight_layout() fig.savefig('W:/Analyses/SINKER comparison/SES_Sinker_lags.pdf', dpi=300, bbox_inches='tight') plt.show() # Plot 2: lagging-variable profile colored by correlation fig, ax = plt.subplots(figsize=(13, 5)) ax.plot( x_axis, lag, color="lightgray", linewidth=1.5, zorder=1 ) sc = ax.scatter( x_axis[positive_corr], lag[positive_corr], c=max_corr[positive_corr], cmap="coolwarm", norm=Normalize(vmin=0, vmax=1), s=70, zorder=2 ) ax.set_title( f"SINKER area flux colored by Positive Correlation " f"(window = {win_width})" ) ax.set_ylabel("SINKER area flux") ax.set_xlabel("Date" if dates is not None else "Record") if positive_corr.any(): fig.colorbar( sc, ax=ax, label="Positive correlation" ) fig.tight_layout() fig.savefig('W:/Analyses/SINKER comparison/SES_Sinker_correlations.pdf', dpi=300, bbox_inches='tight') plt.show() # Plot 3: leading and lagging variables with two y-axes fig, ax1 = plt.subplots(figsize=(13, 5)) line1 = ax1.plot( x_axis, lead, color="tab:blue", linewidth=1.5, label="SES modeled POC flux" ) ax1.set_xlabel("Date" if dates is not None else "Record") ax1.set_ylabel( "SES modeled POC flux", color="tab:blue" ) ax1.tick_params(axis="y", labelcolor="tab:blue") ax2 = ax1.twinx() line2 = ax2.plot( x_axis, lag, color="tab:red", linewidth=1.5, label="SINKER area flux" ) ax2.set_ylabel( "SINKER area flux", color="tab:red" ) ax2.tick_params(axis="y", labelcolor="tab:red") lines = line1 + line2 labels = [line.get_label() for line in lines] ax1.legend(lines, labels, loc="best") ax1.set_title("SES and SINKER fluxes") fig.tight_layout() fig.savefig('W:/Analyses/SINKER comparison/SES_Sinker_TS.pdf', dpi=300, bbox_inches='tight') plt.show() results.to_csv("W:/Analyses/SINKER comparison/Lags and correlations.csv") return results # Run the analysis slidingdata = SES_OffShore_SINKER[ ['date', 'total_area_flux', 'POCmodel_mg_daily_avg'] ].copy() slidingdata["date"] = pd.to_datetime(slidingdata["date"]) # Set date as the index before resampling slidingdatadaily = ( slidingdata .set_index("date") .resample("D") .mean() .reset_index() ) myData = slidingdatadaily mySlidingCCF = sliding_ccf( ts_lead=myData["POCmodel_mg_daily_avg"], ts_lag=myData["total_area_flux"], win_width=60, lag_max=30, plot_results=True, dates=myData["date"] )