figures for thesis

plotting/2014_2019_DO/Puget_Sound_volume.py

@@ -0,0 +1,88 @@
+
+# import things
+from subprocess import Popen as Po
+from subprocess import PIPE as Pi
+from matplotlib.markers import MarkerStyle
+import matplotlib.dates as mdates
+import numpy as np
+import xarray as xr
+import csv
+from datetime import datetime, timedelta
+from matplotlib.dates import DateFormatter
+from matplotlib.dates import MonthLocator
+import matplotlib.patches as patches
+from matplotlib.offsetbox import (OffsetImage, AnnotationBbox)
+import matplotlib.image as image
+import pandas as pd
+import cmocean
+import matplotlib.pylab as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import matplotlib.patheffects as PathEffects
+import pinfo
+
+from lo_tools import Lfun, zfun, zrfun
+from lo_tools import plotting_functions as pfun
+
+import sys
+from pathlib import Path
+pth = Path(__file__).absolute().parent.parent.parent.parent / 'LO' / 'pgrid'
+if str(pth) not in sys.path:
+    sys.path.append(str(pth))
+import gfun
+
+Gr = gfun.gstart()
+
+Ldir = Lfun.Lstart()
+
+##############################################################
+##                       USER INPUTS                        ##
+##############################################################
+
+remove_straits = True
+
+# which  model run to look at?
+gtagex = 'cas7_t0_x4b' # long hindcast (anthropogenic)
+
+##############################################################
+##              CALCULATE PUGET SOUND VOLUME                ##
+##############################################################
+
+# get percent hypoxic volume
+fp = Ldir['LOo'] / 'extract' / gtagex / 'box' / ('pugetsoundDO_2013.01.01_2013.12.31.nc')
+ds = xr.open_dataset(fp)
+if remove_straits:
+        print('    Removing Straits...')
+        lat_threshold = 48.14
+        lon_threshold = -122.76
+        # Create a mask for latitudes and longitudes in the Straits
+        mask = (ds['lat_rho'] > lat_threshold) & (ds['lon_rho'] < lon_threshold)
+        # Expand mask dimensions to match 'oxygen' dimensions
+        expanded_mask = mask.expand_dims(ocean_time=len(ds['ocean_time']), s_rho=len(ds['s_rho']))
+        # Apply the mask to the 'oxygen' variable
+        ds['z_w'] = xr.where(expanded_mask, np.nan, ds['z_w'])
+        ds['pm'] = xr.where(expanded_mask, np.nan, ds['pm'])
+        ds['pn'] = xr.where(expanded_mask, np.nan, ds['pn'])
+print('calculating vertical thickness')
+# get S for the whole grid
+Sfp = Ldir['data'] / 'grids' / 'cas7' / 'S_COORDINATE_INFO.csv'
+reader = csv.DictReader(open(Sfp))
+S_dict = {}
+for row in reader:
+    S_dict[row['ITEMS']] = row['VALUES']
+S = zrfun.get_S(S_dict)
+# get cell thickness
+h = ds['h'].values # height of water column
+z_rho, z_w = zrfun.get_z(h, 0*h, S) 
+dzr = np.diff(z_w, axis=0) # vertical thickness of all cells [m]
+print(dzr.shape)
+
+DZ = dzr/1000
+DX = (ds.pm.values)**-1
+DY = (ds.pn.values)**-1
+DA = DX*DY*(1/1000)*(1/1000) # get area, but convert from m^2 to km^2
+
+# calculate volume of Puget Sound
+print('Puget Sound volume with straits omitted [km^3]')
+vol_timeseries = np.sum(DZ * DA, axis=(1, 2))
+vol = np.nanmean(vol_timeseries)
+print(vol)

plotting/2014_2019_DO/budget_DO/make_terminalinlet_map.py

@@ -44,10 +44,15 @@
 plt.close('all')
 fig = plt.figure(figsize=(7.7,8.7))
 ax = fig.add_subplot(1,1,1)
-ax.tick_params(axis='both', labelsize=10)
+lon_low = -123.42
+lon_high =  -122
+lat_low = 46.93
+lat_high = 48.46
+ax.add_patch(Rectangle((lon_low, lat_low), lon_high-lon_low,lat_high-lat_low, facecolor='#EEEEEE'))
+ax.tick_params(axis='both', labelsize=12)
 plt.xticks(rotation=30)
 # plt.pcolormesh(plon, plat, zm, linewidth=0.5, vmin=-1.2, vmax=0, cmap=plt.get_cmap('Greys'))
-plt.pcolormesh(plon, plat, zm, vmin=-15, vmax=0, cmap=plt.get_cmap(cmocean.cm.ice))
+plt.pcolormesh(plon, plat, zm, vmin=-8, vmax=0, cmap=plt.get_cmap(cmocean.cm.ice))
 
 # get terminal inlet locations
 for stn,station in enumerate(sta_dict): # stations: 
@@ -64,19 +69,21 @@
     inlet_loc[jj,ii] = 20
     # add inlet locations
     # plt.pcolormesh(plon, plat, inlet_loc, linewidth=0.5, vmin=0, vmax=65, cmap=plt.get_cmap('coolwarm'))
-    plt.pcolormesh(plon, plat, inlet_loc, linewidth=0.5, vmin=0, vmax=28, cmap=plt.get_cmap(cmocean.cm.ice))
+    plt.pcolormesh(plon, plat, inlet_loc, linewidth=0.5, vmin=0, vmax=35, cmap=plt.get_cmap(cmocean.cm.ice))
 
 # format
 # ax.axes.xaxis.set_visible(False)
 # ax.axes.yaxis.set_visible(False)
-ax.set_xlim(-123.42, -122) # Puget Sound
-ax.set_ylim(46.93, 48.46) # Puget Sound
+# ax.set_xlim(-123.42, -122) # Puget Sound
+# ax.set_ylim(46.93, 48.46) # Puget Sound
+ax.set_xlim(lon_low, lon_high) # Puget Sound
+ax.set_ylim(lat_low, lat_high) # Puget Sound
 
 pfun.dar(ax)
 
-# # remove border
-# for border in ['top','right','bottom','left']:
-#     ax.spines[border].set_visible(False)
+# remove border
+for border in ['top','right','bottom','left']:
+    ax.spines[border].set_visible(False)
 
 ##########################################################
 ##                 Plot station locations               ##
@@ -152,8 +159,8 @@
 ax.text(-123.37,48.37,'Ecology monitoring\nstation (ID)',va='center',ha='left',fontweight='bold',color='mediumorchid',fontsize = 12)
 ax.text(-123.37,48.3,'No observations',va='center',ha='left',fontweight='bold',color='k',fontsize = 12)
 
-plt.title('Inlet Locations',fontsize = 16)
-fig.tight_layout
+plt.title('Inlet Locations',fontsize = 14)
+plt.tight_layout()
 
 
 # where to put output figures

plotting/2014_2019_DO/budget_DO/plot_budget_TEF_21INLET.py

@@ -2605,11 +2605,14 @@
     deep_lay_DO = zfun.lowpass(deep_lay_DO.values,n=30)
 
     # plot DO colored by flushing time
-    ax.plot(dates_local_daily,deep_lay_DO,linewidth=3.5,color='black', alpha=0.8)
-    ax.plot(dates_local_daily,deep_lay_DO,linewidth=1.7,color='silver')
+    # ax.plot(dates_local_daily,deep_lay_DO,linewidth=3.5,color='black', alpha=0.8)
+    # ax.plot(dates_local_daily,deep_lay_DO,linewidth=1.7,color='silver')
+
+    ax.plot(dates_local_daily,deep_lay_DO,linewidth=1,color='black',alpha=0.5)
 
 # format labels
 ax.set_xlim([dates_local[0],dates_local[-2]])
+ax.set_ylim([0,12])
 ax.set_ylabel(r'DO$_{deep}$ [mg/L]',fontsize=12)
 
 plt.show()

plotting/2014_2019_DO/hypoxic_days_and_volume.py

@@ -11,6 +11,7 @@
 import matplotlib.dates as mdates
 import numpy as np
 import xarray as xr
+import csv
 from datetime import datetime, timedelta
 from matplotlib.dates import DateFormatter
 from matplotlib.dates import MonthLocator
@@ -47,7 +48,7 @@
 DO_thresh = 2 # [mg/L]
 
 # Show WWTP locations?
-WWTP_loc = True
+WWTP_loc = False
 
 remove_straits = True
 
@@ -297,7 +298,7 @@ def LO2SSM_name(rname):
             plt.setp(legend.get_title(),fontsize=20)
 
         # add 10 km bar
-        lat0 = 46.94
+        lat0 = 47.0
         lon0 = -123.05
         lat1 = lat0
         lon1 = -122.91825
@@ -335,7 +336,7 @@ def LO2SSM_name(rname):
 
 # Add colormap title
 plt.suptitle('Days with bottom DO < {} [mg/L]'.format(str(DO_thresh)),
-            fontsize=44, fontweight='bold', y=0.95)
+            fontsize=36, y=0.95)
 
 # Generate plot
 plt.tight_layout
@@ -439,10 +440,16 @@ def LO2SSM_name(rname):
     # plot hypoxic area timeseries
     # ax1.plot(dates_local,hyp_vol[year],color=colors[i],
     #          linewidth=3,alpha=0.5,label=year)
-    ax1.plot(dates_local,hyp_vol[year],color='white',
-             linewidth=3.5)
-    ax1.plot(dates_local,hyp_vol[year],color=colors[i],
-             linewidth=3,label=year)
+    if year == '2017':
+        ax1.plot(dates_local,hyp_vol[year],color='white',
+                linewidth=4,zorder=4)
+        ax1.plot(dates_local,hyp_vol[year],color='black',
+                linewidth=3.5,label=year,zorder=4)
+    else:
+        ax1.plot(dates_local,hyp_vol[year],color='white',
+                linewidth=2.5)
+        ax1.plot(dates_local,hyp_vol[year],color=colors[i],
+                linewidth=2,label=year)
 
 # format figure
 ax1.grid(visible=True, color='w')

plotting/2014_2019_DO/hypoxic_season_avgDOmap.py

@@ -43,7 +43,7 @@
 ##############################################################
 
 # Show WWTP locations?
-WWTP_loc = True
+WWTP_loc = False
 
 remove_straits = False
 
@@ -244,7 +244,7 @@ def LO2SSM_name(rname):
     plt.setp(legend.get_title(),fontsize=20)
 
 # add 10 km bar
-lat0 = 46.94
+lat0 = 47.0
 lon0 = -123.05
 lat1 = lat0
 lon1 = -122.91825
@@ -292,8 +292,8 @@ def LO2SSM_name(rname):
     ax.set_title(letters[i] + ' ' + year + ' - avg', fontsize=28, loc = 'left')
 
 # Add colormap title
-plt.suptitle(start+' to '+end+' average ' + stext + ' ' + vn + ' ' + units,
-            fontsize=44, fontweight='bold', y=0.95)
+plt.suptitle(start+' to '+end+' average ' + stext + ' ' + vn + ' [mg/L]',
+            fontsize=36, y=0.95)
 
 # Generate plot
 plt.tight_layout

plotting/model_bathymetry.py

@@ -6,6 +6,7 @@
 import matplotlib.pyplot as plt
 from matplotlib.patches import Rectangle
 import xarray as xr
+from matplotlib.ticker import MaxNLocator
 import cmocean
 from lo_tools import Lfun, zfun, zrfun
 from lo_tools import plotting_functions as pfun
@@ -43,7 +44,7 @@
 # Initialize figure
 plt.close('all')
 fs = 10
-pfun.start_plot(fs=fs, figsize=(10,7))
+pfun.start_plot(fs=fs, figsize=(11,7))
 fig = plt.figure()
 plt.subplots_adjust(wspace=0, hspace=0)
 # create colormap
@@ -60,12 +61,12 @@
 ax0 = fig.add_subplot(1,2,1)
 # cs = ax0.pcolormesh(plon, plat, zm, vmin=-5, vmax=0, cmap=newcmap)
 cs = ax0.pcolormesh(plon, plat, zm, vmin=-500, vmax=0, cmap=newcmap)
-cbar = plt.colorbar(cs,ax=ax0, location='left', pad=0.05)
-cbar.ax.tick_params(labelsize=11, color='#EEEEEE')#, rotation=30)
-cbar.ax.set_ylabel('Depth [m]', fontsize=11, color='#EEEEEE')
+cbar = plt.colorbar(cs,ax=ax0, location='left', pad=0.1)#pad=0.05)
+cbar.ax.tick_params(labelsize=11)#, color='#EEEEEE')#, rotation=30)
+cbar.ax.set_ylabel('Depth [m]', fontsize=11)#, color='#EEEEEE')
 cbar.outline.set_visible(False)
 cbar_yticks = plt.getp(cbar.ax.axes, 'yticklabels')
-plt.setp(cbar_yticks, color='#EEEEEE')
+plt.setp(cbar_yticks)#, color='#EEEEEE')
 # format figure
 pfun.dar(ax0)
 # pfun.add_coast(ax0, color='gray')
@@ -79,11 +80,15 @@
 print('Lat={},{}'.format(Y[j1],Y[j2]))
 
 # plt.xticks(rotation=30, color='gray')
+plt.xticks(rotation=30,horizontalalignment='right',color='gray',fontsize=12)
+plt.yticks(color='gray',fontsize=12)
+ax0.xaxis.set_major_locator(MaxNLocator(integer=True))
+ax0.yaxis.set_major_locator(MaxNLocator(integer=True))
 # plt.yticks(color='gray')
-ax0.set_yticklabels([])
-ax0.set_xticklabels([])
+# ax0.set_yticklabels([])
+# ax0.set_xticklabels([])
 # add title
-ax0.set_title('Salish Sea',fontsize=16,color='#EEEEEE')
+ax0.set_title('(a) Salish Sea',fontsize=14, loc='left')#,color='#EEEEEE')
 # add 10 km bar
 lat0 = 47
 lon0 = -124.63175
@@ -96,24 +101,24 @@
          horizontalalignment='center')
 # draw box around Puget Sound
 bordercolor = '#EEEEEE'
-ax0.add_patch(Rectangle((-123.2, 46.93), 1.1, 1.52,
+ax0.add_patch(Rectangle((-123.3, 46.93), 1.2, 1.52,
              edgecolor = bordercolor, facecolor='none', lw=2))
 
 # add major cities
 # Seattle
-ax0.scatter([-122.3328],[47.6061],s=[250],color='lightpink',
+ax0.scatter([-122.3328],[47.6061],s=[250],color='pink',
             marker='*',edgecolors='darkred')
-ax0.text(-122.3328 + 0.1,47.6061,'Seattle',color='lightpink', rotation=90,
+ax0.text(-122.3328 + 0.1,47.6061,'Seattle',color='pink', rotation=90,
          horizontalalignment='left',verticalalignment='center', size=12)
 # Vancouver
-ax0.scatter([-123.1207],[49.2827],s=[250],color='lightpink',
+ax0.scatter([-123.1207],[49.2827],s=[250],color='pink',
             marker='*',edgecolors='darkred')
-ax0.text(-123.1207 + 0.1,49.2827,'Vancouver',color='lightpink', rotation=0,
+ax0.text(-123.1207 + 0.1,49.2827,'Vancouver',color='pink', rotation=0,
          horizontalalignment='left',verticalalignment='center', size=12)
 
 # add major water bodies
-ax0.text(-124.937095,47.827238,'Pacific Ocean',color='black', rotation=-75,
-         horizontalalignment='left',verticalalignment='center', size=12)
+ax0.text(-124.937095,47.782238,'Pacific Ocean',color='black', rotation=-75,
+         horizontalalignment='left',verticalalignment='center', size=12,fontweight='bold')
 
 # Puget Sound ----------------------------------------------------------
 ax1 = fig.add_subplot(1,2,2)
@@ -125,11 +130,11 @@
 # cs = ax1.pcolormesh(plon, plat, zm, vmin=-5, vmax=0, cmap=newcmap)
 cs = ax1.pcolormesh(plon, plat, zm, vmin=-250, vmax=0, cmap=newcmap)
 cbar = plt.colorbar(cs,ax=ax1, location='right', pad=0.05)
-cbar.ax.tick_params(labelsize=11, color='#EEEEEE')#, rotation=30)
-cbar.ax.set_ylabel('Depth [m]', fontsize=11, color='#EEEEEE')
+cbar.ax.tick_params(labelsize=11)#, color='#EEEEEE')#, rotation=30)
+cbar.ax.set_ylabel('Depth [m]', fontsize=11)#, color='#EEEEEE')
 cbar.outline.set_visible(False)
 cbar_yticks = plt.getp(cbar.ax.axes, 'yticklabels')
-plt.setp(cbar_yticks, color='#EEEEEE')
+plt.setp(cbar_yticks)#, color='#EEEEEE')
 # format figure
 pfun.dar(ax1)
 # pfun.add_coast(ax1, color='gray')
@@ -142,10 +147,13 @@
 ax1.set_ylim([46.93,48.45])
 ax1.set_yticklabels([])
 ax1.set_xticklabels([])
+# plt.xticks(rotation=30,horizontalalignment='right',color='gray')
+# ax1.xaxis.set_major_locator(MaxNLocator(integer=True))
+# ax1.yaxis.set_major_locator(MaxNLocator(integer=True))
 # plt.xticks(rotation=30, color='gray')
 # plt.yticks(color='gray')
 # add title
-ax1.set_title('Puget Sound',fontsize=16, color='#EEEEEE')
+ax1.set_title('(b) Puget Sound',fontsize=14,loc='left')# color='#EEEEEE')
 # add 10 km bar
 lat0 = 47
 lon0 = -123.2
@@ -158,4 +166,4 @@
          horizontalalignment='center')
 
 plt.subplots_adjust(hspace = 0.01)
-plt.savefig(out_dir / ('model_bathy.png'),transparent='True')
+plt.savefig(out_dir / ('model_bathy.png'))#,transparent='True')