budget figures refactored

ajleeson · Jan 14, 2025 · 99d30ba · 99d30ba
1 parent 218ff46
commit 99d30ba
Show file tree

Hide file tree

Showing 5 changed files with 77 additions and 82 deletions.
diff --git 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
@@ -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
+    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
diff --git a/terminal_inlet_DO/scripts/main.py b/terminal_inlet_DO/scripts/main.py
@@ -6,43 +6,20 @@
 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
 import QinDOin_correl_consumption
 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
@@ -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
diff --git 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()