diff --git a/terminal_inlet_DO/scripts/budget_barchart.py b/terminal_inlet_DO/scripts/budget_barchart.py index f6cf5a5..02bb94c 100644 --- a/terminal_inlet_DO/scripts/budget_barchart.py +++ b/terminal_inlet_DO/scripts/budget_barchart.py @@ -147,9 +147,7 @@ def budget_barchart(inlets,shallowlay_dict,deeplay_dict, # create a new dictionary of results oxy_dict = {} hyp_dict = {} - oxy_dict_std = {} - hyp_dict_std = {} - + for inlet in inlets: for attribute, measurement in deeplay_dict[inlet].items(): # skip variables we are not interested in @@ -185,8 +183,8 @@ def budget_barchart(inlets,shallowlay_dict,deeplay_dict, print(' alpha = 0.05') a = oxy_dict[attribute] b = hyp_dict[attribute] - ttest = ttest_ind(a, b, axis=0, equal_var=False) - print(' '+str(ttest)) + ttest,p_value = ttest_ind(a, b, axis=0, equal_var=False) + print(' p = {}'.format(round(p_value,3))) else: continue print('\n') @@ -229,14 +227,13 @@ def budget_barchart(inlets,shallowlay_dict,deeplay_dict, wiggle = -0.7 ha = 'center' # plot - offset = width * multiplier_deep1 if attribute == 'Storage': hatchcolor = 'white' else: hatchcolor = 'white' if i == 0: - rects = ax[2].bar(pos + shift, avg, width, zorder=5, align='center', edgecolor=hatchcolor,color=color, hatch='xx') - rects = ax[2].bar(pos + shift, avg, width, zorder=5, align='center', edgecolor=color,color='none') + ax[2].bar(pos + shift, avg, width, zorder=5, align='center', edgecolor=hatchcolor,color=color, hatch='xx') + ax[2].bar(pos + shift, avg, width, zorder=5, align='center', edgecolor=color,color='none') if attribute == 'd/dt(DO)': ax[2].text(pos + shift, 0+wiggle, str(round(avg,3)),horizontalalignment=ha,verticalalignment='center', color=color,fontsize=10, fontweight='bold', rotation=45) @@ -245,8 +242,7 @@ def budget_barchart(inlets,shallowlay_dict,deeplay_dict, color='gray',fontsize=10, rotation=45) multiplier_deep1 += 2 elif i == 1: - offset = width * multiplier_deep2 - rects = ax[2].bar(pos + shift, avg, width, zorder=5, align='center', edgecolor=color,color=color, label=label) + ax[2].bar(pos + shift, avg, width, zorder=5, align='center', edgecolor=color,color=color, label=label) if attribute == 'd/dt(DO)': ax[2].text(pos+shift, 0+wiggle, str(round(avg,3)),horizontalalignment=ha,verticalalignment='center', color=color,fontsize=10, fontweight='bold', rotation=45) @@ -306,8 +302,8 @@ def budget_barchart(inlets,shallowlay_dict,deeplay_dict, print(' alpha = 0.05') a = oxy_dict[attribute] b = hyp_dict[attribute] - ttest = ttest_ind(a, b, axis=0, equal_var=False) - print(' '+str(ttest)) + ttest,p_value = ttest_ind(a, b, axis=0, equal_var=False) + print(' p = {}'.format(round(p_value,3))) else: continue @@ -340,11 +336,10 @@ def budget_barchart(inlets,shallowlay_dict,deeplay_dict, if avg > 0: wiggle = -0.04 # plot - offset = width * multiplier_deep1 hatchcolor = 'white' if i == 0: - rects = ax[3].bar(pos + shift, avg, width, zorder=5, edgecolor=hatchcolor,color=color, hatch='xx') - rects = ax[3].bar(pos + shift, avg, width, zorder=5, edgecolor=color,color='none') + ax[3].bar(pos + shift, avg, width, zorder=5, edgecolor=hatchcolor,color=color, hatch='xx') + ax[3].bar(pos + shift, avg, width, zorder=5, edgecolor=color,color='none') if attribute == 'd/dt(DO)': ax[3].text(pos+shift, 0+wiggle, str(round(avg,3)),horizontalalignment='center',verticalalignment='center', color=color,fontsize=10, fontweight='bold') @@ -353,8 +348,7 @@ def budget_barchart(inlets,shallowlay_dict,deeplay_dict, color='gray',fontsize=10) multiplier_deep1 += 2 elif i == 1: - offset = width * multiplier_deep2 - rects = ax[3].bar(pos + shift, avg, width, zorder=5, edgecolor=color,color=color,label=label) + ax[3].bar(pos + shift, avg, width, zorder=5, edgecolor=color,color=color,label=label) if attribute == 'd/dt(DO)': ax[3].text(pos+shift, 0+wiggle, str(round(avg,3)),horizontalalignment='center',verticalalignment='center', color=color,fontsize=10, fontweight='bold') diff --git a/terminal_inlet_DO/scripts/helper_functions.py b/terminal_inlet_DO/scripts/helper_functions.py index bdf990a..24f0e25 100644 --- a/terminal_inlet_DO/scripts/helper_functions.py +++ b/terminal_inlet_DO/scripts/helper_functions.py @@ -5,6 +5,7 @@ """ import numpy as np +import pytz def lowpass(data, f='hanning', n=40, nanpad=True): """ @@ -64,4 +65,12 @@ def hanning_shape(n=40): ff = np.cos(np.linspace(-np.pi,np.pi,n+2))[1:-1] filt = (1 + ff)/2 filt = filt / filt.sum() - return filt \ No newline at end of file + return filt + +def get_dt_local(dt, tzl='US/Pacific'): + # take a model datetime (assumed to be UTC) and return local datetime + tz_utc = pytz.timezone('UTC') + tz_local = pytz.timezone(tzl) + dt_utc = dt.replace(tzinfo=tz_utc) + dt_local = dt_utc.astimezone(tz_local) + return dt_local \ No newline at end of file diff --git a/terminal_inlet_DO/scripts/main.py b/terminal_inlet_DO/scripts/main.py index 7a56a43..cb50793 100644 --- a/terminal_inlet_DO/scripts/main.py +++ b/terminal_inlet_DO/scripts/main.py @@ -6,31 +6,12 @@ January 2025 """ -from subprocess import Popen as Po -from subprocess import PIPE as Pi -from matplotlib.markers import MarkerStyle -import matplotlib.dates as mdates -import matplotlib.ticker as ticker -import numpy as np -import xarray as xr -from datetime import datetime, timedelta import pandas as pd -import matplotlib.patches as mpatches -from scipy.linalg import lstsq -import math -import matplotlib.patches as patches -from matplotlib.colors import ListedColormap -import csv -import cmocean -from scipy.stats import pearsonr -from scipy.stats import ttest_ind -import pingouin as pg import matplotlib.pylab as plt -import gsw import pickle -# import get_two_layer # import helper functions +import helper_functions import get_monthly_means import budget_error import budget_barchart @@ -38,11 +19,7 @@ import plot_monthly_means import dodeep_hypvol_timeseries import net_decrease_boxplots - -from lo_tools import Lfun, zfun -from lo_tools import plotting_functions as pfun - -Ldir = Lfun.Lstart() +import multiple_regression plt.close('all') @@ -107,10 +84,10 @@ # create time_vector dates_hrly = pd.date_range(start= startdate, end=enddate_hrly, freq= 'h') -dates_local_hrly = [pfun.get_dt_local(x) for x in dates_hrly] +dates_local_hrly = [helper_functions.get_dt_local(x) for x in dates_hrly] # crop time vector (because we only have jan 2 - dec 30) dates_daily = pd.date_range(start= startdate, end=enddate, freq= 'd')[2::] -dates_local_daily = [pfun.get_dt_local(x) for x in dates_daily] +dates_local_daily = [helper_functions.get_dt_local(x) for x in dates_daily] ########################################################## ## Get monthly means ## @@ -186,39 +163,10 @@ net_decrease_boxplots.net_decrease_boxplots(dimensions_dict,deeplay_dict, minday,maxday) -# ########################################################## -# ## Multiple regression (1) ## -# ########################################################## - -# # create array of predictors - -# input_array = np.array([mean_DOin, mean_Tflush, [1]*len(mean_DOin)]).T - - -# B,a,b,c = lstsq(input_array,deep_lay_DO) -# slope_DOin = B[0] -# slope_Tflush = B[1] -# intercept = B[2] - -# print('\nMean deep layer DO [mg/L] = {}*DOin + {}*Tflush + {}'.format( -# round(slope_DOin,2),round(slope_Tflush,2),round(intercept,2))) - - -# # calculate r^2 and p value -# r,p = pearsonr(mean_DOin,deep_lay_DO) -# print('\n===============================================') -# print('DO_deep dependence on DO_in') -# print(' r = {}'.format(r)) -# print(' R^2 = {}'.format(r**2)) -# print(' p = {}'.format(p)) -# print('===============================================\n') +########################################################## +## Multiple regression ## +########################################################## -# # calculate r^2 and p value -# predicted_DOdeep = slope_DOin * mean_DOin + slope_Tflush * mean_Tflush + intercept -# r,p = pearsonr(deep_lay_DO,predicted_DOdeep) -# print('\n===============================================') -# print('DO_deep dependence on DO_in and T_flush') -# print(' r = {}'.format(r)) -# print(' R^2 = {}'.format(r**2)) -# print(' p = {}'.format(p)) -# print('===============================================\n') +multiple_regression.multiple_regression(MONTHLYmean_DOdeep, + MONTHLYmean_DOin, + MONTHLYmean_Tflush) diff --git a/terminal_inlet_DO/scripts/multiple_regression.py b/terminal_inlet_DO/scripts/multiple_regression.py new file mode 100644 index 0000000..9545cff --- /dev/null +++ b/terminal_inlet_DO/scripts/multiple_regression.py @@ -0,0 +1,44 @@ +""" +Calculate multiple linear regression +of DOdeep dependece on DOin and Tflush +""" +import numpy as np +from scipy.linalg import lstsq +from scipy.stats import pearsonr + +def multiple_regression(MONTHLYmean_DOdeep, + MONTHLYmean_DOin, + MONTHLYmean_Tflush): + + + print('\n=====================Multiple Linear Regression=======================\n') + + # create array of predictors + input_array = np.array([MONTHLYmean_DOin, MONTHLYmean_Tflush, [1]*len(MONTHLYmean_DOin)]).T + + + B,a,b,c = lstsq(input_array,MONTHLYmean_DOdeep) + slope_DOin = B[0] + slope_Tflush = B[1] + intercept = B[2] + + print('\nMean deep layer DO [mg/L] = {}*DOin + {}*Tflush + {}\n'.format( + round(slope_DOin,2),round(slope_Tflush,2),round(intercept,2))) + + + # calculate r^2 and p value + r,p = pearsonr(MONTHLYmean_DOin,MONTHLYmean_DOdeep) + print('DO_deep dependence on DO_in') + print(' r = {}'.format(round(r,3))) + print(' R^2 = {}'.format(round((r**2),3))) + print(' p = {:.2e}'.format(p)) + + # calculate r^2 and p value + predicted_DOdeep = slope_DOin * MONTHLYmean_DOin + slope_Tflush * MONTHLYmean_Tflush + intercept + r,p = pearsonr(MONTHLYmean_DOdeep,predicted_DOdeep) + print('\nDO_deep dependence on DO_in and T_flush') + print(' r = {}'.format(round(r,3))) + print(' R^2 = {}'.format(round((r**2),3))) + print(' p = {:.2e}'.format(p)) + + return \ No newline at end of file diff --git a/terminal_inlet_DO/scripts/net_decrease_boxplots.py b/terminal_inlet_DO/scripts/net_decrease_boxplots.py index 1d126e8..57bfaea 100644 --- a/terminal_inlet_DO/scripts/net_decrease_boxplots.py +++ b/terminal_inlet_DO/scripts/net_decrease_boxplots.py @@ -91,14 +91,14 @@ def net_decrease_boxplots(dimensions_dict,deeplay_dict, 'inlet': np.repeat(stations_sorted, repeats=repeats)}) pingu = pg.welch_anova(data=df, dv='storage', between='inlet') print('Welch\'s ANOVA test of all thirteen inlets') - print(pingu) + print(' p = {:.2e}'.format((pingu['p-unc'].values[0]))) print('\n') # remove Dabob Bay from analysis and repeat ANOVA df_nodabob = df[df['inlet'] != 'dabob'] pingu_nodabob = pg.welch_anova(data=df_nodabob, dv='storage', between='inlet') print('OMITTING DABOB BAY: Welch\'s ANOVA test of all other twelve inlets') - print(pingu_nodabob) + print(' p = {}'.format(round(pingu_nodabob['p-unc'].values[0],3))) print('\n') plt.tight_layout()