Skip to content

Commit

Permalink
Merge pull request #85 from MOchiara/chiara-patch-24
Browse files Browse the repository at this point in the history 
Chiara patch 24
  • Loading branch information
@MOchiara
MOchiara authored Nov 11, 2024
2 parents db055c0 + 3fcff37 commit 079f404
Showing 1 changed file with 84 additions and 77 deletions.
161 changes: 84 additions & 77 deletions glidertest/tools.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
import gsw
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import warnings


def _necessary_variables_check(ds: xr.Dataset, vars: list):
Expand Down Expand Up @@ -89,9 +90,10 @@ def updown_bias(ds, var='PSAL', v_res=1):
varG, profG, depthG = grid2d(ds.PROFILE_NUMBER, ds.DEPTH, ds[var], p, z)

grad = np.diff(varG, axis=0) # Horizontal gradients

dc = np.nanmean(grad[0::2, :], axis=0) # Dive - CLimb
cd = np.nanmean(grad[1::2, :], axis=0) # Climb - Dive
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=RuntimeWarning)
dc = np.nanmean(grad[0::2, :], axis=0) # Dive - CLimb
cd = np.nanmean(grad[1::2, :], axis=0) # Climb - Dive

df = pd.DataFrame(data={'dc': dc, 'cd': cd, 'depth': depthG[0, :]})
else:
Expand Down Expand Up @@ -140,7 +142,9 @@ def find_cline(var, depth_array):
Find the depth of the maximum vertical difference for a specified variables
Input data has to be gridded
"""
clin = np.where(np.abs(np.diff(np.nanmean(var, axis=0))) == np.nanmax(np.abs(np.diff(np.nanmean(var, axis=0)))))
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=RuntimeWarning)
clin = np.where(np.abs(np.diff(np.nanmean(var, axis=0))) == np.nanmax(np.abs(np.diff(np.nanmean(var, axis=0)))))
return np.round(depth_array[0, clin[0]], 1)


Expand Down Expand Up @@ -187,55 +191,56 @@ def plot_basic_vars(ds, v_res=1, start_prof=0, end_prof=-1):
pycno = find_cline(denG, depthG)
print(
f'The thermocline, halocline and pycnocline are located at respectively {thermo}, {halo} and {pycno}m as shown in the plots as well')

fig, ax = plt.subplots(1, 2, figsize=(15, 5))
ax1 = ax[0].twiny()
ax2 = ax[0].twiny()
ax2.spines["top"].set_position(("axes", 1.2))
ax[0].plot(np.nanmean(tempG, axis=0), depthG[0, :], c='blue')
ax1.plot(np.nanmean(salG, axis=0), depthG[0, :], c='red')
ax2.plot(np.nanmean(denG, axis=0), depthG[0, :], c='black')
ax[0].axhline(thermo, linestyle='dashed', c='blue')
ax1.axhline(halo, linestyle='dashed', c='red')
ax2.axhline(pycno, linestyle='dashed', c='black')

ax[0].set(xlabel=f'Average Temperature [C] \nbetween profile {start_prof} and {end_prof}', ylabel='Depth (m)')
ax[0].tick_params(axis='x', colors='blue')
ax[0].xaxis.label.set_color('blue')
ax1.spines['bottom'].set_color('blue')
ax1.set(xlabel=f'Average Salinity [PSU] \nbetween profile {start_prof} and {end_prof}')
ax1.xaxis.label.set_color('red')
ax1.spines['top'].set_color('red')
ax1.tick_params(axis='x', colors='red')
ax2.spines['bottom'].set_color('black')
ax2.set(xlabel=f'Average Density [kg m-3] \nbetween profile {start_prof} and {end_prof}')
ax2.xaxis.label.set_color('black')
ax2.spines['top'].set_color('black')
ax2.tick_params(axis='x', colors='black')

if 'CHLA' in ds.variables:
chlaG, profG, depthG = grid2d(ds.PROFILE_NUMBER, ds.DEPTH, ds.CHLA, p, z)
chlaG = chlaG[start_prof:end_prof, :]
ax2_1 = ax[1].twiny()
ax2_1.plot(np.nanmean(chlaG, axis=0), depthG[0, :], c='green')
ax2_1.set(xlabel=f'Average Chlorophyll-a [mg m-3] \nbetween profile {start_prof} and {end_prof}')
ax2_1.xaxis.label.set_color('green')
ax2_1.spines['top'].set_color('green')
ax2_1.tick_params(axis='x', colors='green')
else:
ax[1].text(0.3, 0.7, 'Chlorophyll data unavailable', va='top', transform=ax[1].transAxes)

if 'DOXY' in ds.variables:
oxyG, profG, depthG = grid2d(ds.PROFILE_NUMBER, ds.DEPTH, ds.DOXY, p, z)
oxyG = oxyG[start_prof:end_prof, :]
ax[1].plot(np.nanmean(oxyG, axis=0), depthG[0, :], c='orange')
ax[1].set(xlabel=f'Average Oxygen [mmol m-3] \nbetween profile {start_prof} and {end_prof}')
ax[1].xaxis.label.set_color('orange')
ax[1].spines['top'].set_color('orange')
ax[1].tick_params(axis='x', colors='orange')
ax[1].spines['bottom'].set_color('orange')
else:
ax[1].text(0.3, 0.5, 'Oxygen data unavailable', va='top', transform=ax[1].transAxes)
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=RuntimeWarning)
fig, ax = plt.subplots(1, 2, figsize=(15, 5))
ax1 = ax[0].twiny()
ax2 = ax[0].twiny()
ax2.spines["top"].set_position(("axes", 1.2))
ax[0].plot(np.nanmean(tempG, axis=0), depthG[0, :], c='blue')
ax1.plot(np.nanmean(salG, axis=0), depthG[0, :], c='red')
ax2.plot(np.nanmean(denG, axis=0), depthG[0, :], c='black')
ax[0].axhline(thermo, linestyle='dashed', c='blue')
ax1.axhline(halo, linestyle='dashed', c='red')
ax2.axhline(pycno, linestyle='dashed', c='black')

ax[0].set(xlabel=f'Average Temperature [C] \nbetween profile {start_prof} and {end_prof}', ylabel='Depth (m)')
ax[0].tick_params(axis='x', colors='blue')
ax[0].xaxis.label.set_color('blue')
ax1.spines['bottom'].set_color('blue')
ax1.set(xlabel=f'Average Salinity [PSU] \nbetween profile {start_prof} and {end_prof}')
ax1.xaxis.label.set_color('red')
ax1.spines['top'].set_color('red')
ax1.tick_params(axis='x', colors='red')
ax2.spines['bottom'].set_color('black')
ax2.set(xlabel=f'Average Density [kg m-3] \nbetween profile {start_prof} and {end_prof}')
ax2.xaxis.label.set_color('black')
ax2.spines['top'].set_color('black')
ax2.tick_params(axis='x', colors='black')

if 'CHLA' in ds.variables:
chlaG, profG, depthG = grid2d(ds.PROFILE_NUMBER, ds.DEPTH, ds.CHLA, p, z)
chlaG = chlaG[start_prof:end_prof, :]
ax2_1 = ax[1].twiny()
ax2_1.plot(np.nanmean(chlaG, axis=0), depthG[0, :], c='green')
ax2_1.set(xlabel=f'Average Chlorophyll-a [mg m-3] \nbetween profile {start_prof} and {end_prof}')
ax2_1.xaxis.label.set_color('green')
ax2_1.spines['top'].set_color('green')
ax2_1.tick_params(axis='x', colors='green')
else:
ax[1].text(0.3, 0.7, 'Chlorophyll data unavailable', va='top', transform=ax[1].transAxes)

if 'DOXY' in ds.variables:
oxyG, profG, depthG = grid2d(ds.PROFILE_NUMBER, ds.DEPTH, ds.DOXY, p, z)
oxyG = oxyG[start_prof:end_prof, :]
ax[1].plot(np.nanmean(oxyG, axis=0), depthG[0, :], c='orange')
ax[1].set(xlabel=f'Average Oxygen [mmol m-3] \nbetween profile {start_prof} and {end_prof}')
ax[1].xaxis.label.set_color('orange')
ax[1].spines['top'].set_color('orange')
ax[1].tick_params(axis='x', colors='orange')
ax[1].spines['bottom'].set_color('orange')
else:
ax[1].text(0.3, 0.5, 'Oxygen data unavailable', va='top', transform=ax[1].transAxes)

[a.set_ylim(depthG.max() + 10, -5) for a in ax]
[a.grid() for a in ax]
Expand Down Expand Up @@ -340,7 +345,7 @@ def sunset_sunrise(time, lat, lon):

bluffton = []
for i in range(len(grp_avg.lat)):
bluffton.append(api.wgs84.latlon(grp_avg.lat[i], grp_avg.lon[i]))
bluffton.append(api.wgs84.latlon(grp_avg.lat.iloc[i], grp_avg.lon.iloc[i]))
bluffton = np.array(bluffton)

sunrise = []
Expand Down Expand Up @@ -662,9 +667,9 @@ def plot_profIncrease(ds: xr.DataArray, ax: plt.Axes = None, **kw: dict, ) -> tu
ax[1].scatter(ds.TIME[np.where((np.diff(ds.PROFILE_NUMBER) != 0) & (np.diff(ds.PROFILE_NUMBER) != 1))],
ds.DEPTH[np.where((np.diff(ds.PROFILE_NUMBER) != 0) & (np.diff(ds.PROFILE_NUMBER) != 1))],
s=10, c='red', label='Depth at which we have issues \n with the profile number assigned')
ax[1].legend()
ax[1].set(ylabel='Depth')
ax[1].invert_yaxis()
ax[1].legend()
ax[1].xaxis.set_major_locator(plt.MaxNLocator(8))
[a.grid() for a in ax]
return fig, ax
Expand Down Expand Up @@ -1145,29 +1150,31 @@ def findbetw(arr, bounds):
w_upper = np.zeros(len(zgrid))

# Cycle through the bins and calculate the mean vertical velocity
for zdo in range(len(zgrid)):
z1 = bin_edges[zdo]
z2 = bin_edges[zdo + 1]
ifind = findbetw(depth, [z1, z2])

CIlimits = .95 # Could be passed as a variable. 0.95 for 95% confidence intervals
if len(ifind):
meanw[zdo] = np.nanmean(ww[ifind])

# Confidence intervals
# Number of data points used in the mean at this depth (zgrid[zdo])
NNz[zdo] = len(ifind)
if NNz[zdo] > 1:
se = np.nanstd(ww[ifind]) / np.sqrt(NNz[zdo]) # Standard error
ci = se * stats.t.ppf((1 + CIlimits) / 2, NNz[zdo] - 1) # Confidence interval based on CIlimits
w_lower[zdo] = meanw[zdo] - ci
w_upper[zdo] = meanw[zdo] + ci
else:
w_lower[zdo] = np.nan
w_upper[zdo] = np.nan
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=RuntimeWarning)
for zdo in range(len(zgrid)):
z1 = bin_edges[zdo]
z2 = bin_edges[zdo + 1]
ifind = findbetw(depth, [z1, z2])

CIlimits = .95 # Could be passed as a variable. 0.95 for 95% confidence intervals
if len(ifind):
meanw[zdo] = np.nanmean(ww[ifind])

# Confidence intervals
# Number of data points used in the mean at this depth (zgrid[zdo])
NNz[zdo] = len(ifind)
if NNz[zdo] > 1:
se = np.nanstd(ww[ifind]) / np.sqrt(NNz[zdo]) # Standard error
ci = se * stats.t.ppf((1 + CIlimits) / 2, NNz[zdo] - 1) # Confidence interval based on CIlimits
w_lower[zdo] = meanw[zdo] - ci
w_upper[zdo] = meanw[zdo] + ci
else:
w_lower[zdo] = np.nan
w_upper[zdo] = np.nan

else:
meanw[zdo] = np.nan
else:
meanw[zdo] = np.nan


# Package the outputs into an xarray dataset
Expand Down

0 comments on commit 079f404

Please sign in to comment.