diff --git a/glidertest/tools.py b/glidertest/tools.py index fcc8ab6..83c5040 100644 --- a/glidertest/tools.py +++ b/glidertest/tools.py @@ -782,8 +782,8 @@ def plot_ts_histograms(ds: xr.Dataset, ax: plt.Axes = None, **kw: dict) -> tuple CT = gsw.CT_from_t(SA, temp, depth) # Reduce to middle 95% of values - temp_filtered = CT[(CT >= np.nanpercentile(temp, 2.5)) & (CT <= np.nanpercentile(CT, 97.5))] - sal_filtered = SA[(SA >= np.nanpercentile(sal, 2.5)) & (SA <= np.nanpercentile(sal, 97.5))] + temp_filtered = CT[(CT >= np.nanpercentile(CT, 2.5)) & (CT <= np.nanpercentile(CT, 97.5))] + sal_filtered = SA[(SA >= np.nanpercentile(SA, 2.5)) & (SA <= np.nanpercentile(SA, 97.5))] ax[0].hist(temp_filtered, bins=50, **kw) ax[0].set_xlabel('Conservative Temperature (°C)') @@ -817,4 +817,226 @@ def plot_ts_histograms(ds: xr.Dataset, ax: plt.Axes = None, **kw: dict) -> tuple ax[2].contour(xi, yi, zi, colors='black', alpha=0.5, linewidths=0.5) ax[2].clabel(ax[2].contour(xi, yi, zi, colors='black', alpha=0.5, linewidths=0.5), inline=True, fontsize=10) - return fig, ax \ No newline at end of file + return fig, ax + +def calc_DEPTH_Z(ds): + """ + Calculate the depth (Z position) of the glider using the gsw library to convert pressure to depth. + + Parameters + ---------- + ds (xarray.Dataset): The input dataset containing 'PRES', 'LATITUDE', and 'LONGITUDE' variables. + + Returns + ------- + xarray.Dataset: The dataset with an additional 'DEPTH_Z' variable. + """ + # Ensure the required variables are present + if 'PRES' not in ds.variables or 'LATITUDE' not in ds.variables or 'LONGITUDE' not in ds.variables: + raise ValueError("Dataset must contain 'PRES', 'LATITUDE', and 'LONGITUDE' variables.") + + # Initialize the new variable with the same dimensions as dive_num + ds['DEPTH_Z'] = (['N_MEASUREMENTS'], np.full(ds.dims['N_MEASUREMENTS'], np.nan)) + + # Calculate depth using gsw + depth = gsw.z_from_p(ds['PRES'], ds['LATITUDE']) + ds['DEPTH_Z'] = depth + + # Assign the calculated depth to a new variable in the dataset + ds['DEPTH_Z'].attrs = { + "units": "meters", + "positive": "up", + "standard_name": "depth", + "comment": "Depth calculated from pressure using gsw library, positive up.", + } + + return ds + +def calc_glider_w_from_depth(ds): + """ + Calculate the vertical velocity of a glider using changes in pressure with time. + + Parameters + ---------- + ds (xarray.Dataset): Dataset containing 'DEPTH' and 'TIME'. + - DEPTH (array-like): Array of depth measurements + - TIME (array-like): Array of time stamps + + Returns + ------- + ds (xarray.Dataset): Containing the new variable + - GLIDER_VERT_VELO_DZDT (array-like): with vertical velocities calculated from dz/dt + """ + # Ensure inputs are numpy arrays + time = ds.TIME.values + if 'DEPTH_Z' not in ds.variables and all(var in ds.variables for var in ['PRES', 'LATITUDE', 'LONGITUDE']): + ds = calc_DEPTH_Z(ds) + depth = ds.DEPTH_Z.values + + # Calculate the centered differences in pressure and time, i.e. instead of using neighboring points, + # use the points two steps away. This has a couple of advantages: one being a slight smoothing of the + # differences, and the other that the calculated speed will be the speed at the midpoint of the two + # points. + # For data which are evenly spaced in time, this will be equivalent to a centered difference. + # For data which are not evenly spaced in time, i.e. when a Seaglider sample rate changes from 5 + # seconds to 10 seconds, there may be some uneven weighting of the differences. + delta_z_meters = (depth[2:] - depth[:-2]) + delta_time_datetime64ns = (time[2:] - time[:-2]) + delta_time_sec = delta_time_datetime64ns / np.timedelta64(1, 's') # Convert to seconds + + # Calculate vertical velocity (rate of change of pressure with time) + vertical_velocity = delta_z_meters / delta_time_sec + + # Pad the result to match the original array length + vertical_velocity = np.pad(vertical_velocity, (1, 1), 'edge') + + # No - Convert vertical velocity from m/s to cm/s + vertical_velocity = vertical_velocity + + # Add vertical velocity to the dataset + ds = ds.assign(GLIDER_VERT_VELO_DZDT=(('N_MEASUREMENTS'), vertical_velocity, {'long_name': 'glider_vertical_speed_from_pressure', 'units': 'm s-1'})) + + return ds + +def calc_seawater_w(ds): + """ + Calculate the vertical seawater velocity and add it to the dataset. + + Parameters + ---------- + ds (xarray.Dataset): Dataset containing 'VERT_GLIDER_SPEED' and 'VERT_SPEED_DZDT'. + + Returns + ------- + ds (xarray.Dataset): Dataset with the new variable 'VERT_SW_SPEED', which is the inferred vertical seawater velocity. + + Eleanor's note: This could be bundled with calc_glider_w_from_depth, but keeping them separate allows for some extra testing/flexibility for the user. + """ + # Check if 'VERT_GLIDER_SPEED' is in the dataset + if 'GLIDER_VERT_VELO_MODEL' not in ds: + print("Error: 'GLIDER_VERT_VELO_MODEL' is not in the dataset.") + return ds + + # Calculate the vertical seawater velocity + vert_sw_speed = ds['GLIDER_VERT_VELO_DZDT'].values - ds['GLIDER_VERT_VELO_MODEL'].values + + # Add vertical seawater velocity to the dataset as a data variable + ds = ds.assign(VERT_CURR_MODEL=(('N_MEASUREMENTS'), vert_sw_speed, {'long_name': 'vertical_current_of_seawater_derived_from_glider_flight_model', 'units': 'm s-1'})) + return ds + + +def plot_vertical_speeds_with_histograms(ds, start_prof=None, end_prof=None): + """ + Plot vertical speeds with histograms for diagnostic purposes. + This function generates a diagnostic plot for the calculation of vertical seawater velocity. + It plots the modelled and computed (from dz/dt) vertical velocities as line plots, where these + should be similar. The difference between these velocities is the implied seawater velocity, + which should be closer to zero than the vehicle velocities. The histogram provides a visual + representation to identify any biases. The final calculation of the median should be close to + zero if a large enough sample of dives is input and if the glider flight model is well-tuned. + + Parameters + ---------- + ds (xarray.Dataset): The dataset containing the vertical speed data where + - VERT_GLIDER_SPEED is the modelled glider speed + - VERT_SPEED_DZDT is the computed glider speed from the pressure sensor + - VERT_SW_SPEED is the implied seawater velocity. + start_prof (int, optional): The starting profile number for subsetting the dataset. Defaults to first profile number. + end_prof (int, optional): The ending profile number for subsetting the dataset. Defaults to last profile number. + + Returns + ------- + fig, axs (tuple): The figure and axes objects for the plot. + """ + required_vars = ['GLIDER_VERT_VELO_MODEL', 'GLIDER_VERT_VELO_DZDT', 'VERT_CURR_MODEL'] + for var in required_vars: + if var not in ds: + print(f"Dataset must contain '{var}' to create this plot.") + return + + if start_prof is not None and end_prof is not None: + # Subset the dataset for the given profile range + ds = ds.where((ds['PROFILE_NUMBER'] >= start_prof) & (ds['PROFILE_NUMBER'] <= end_prof), drop=True) + + if start_prof is None: + start_prof = ds['PROFILE_NUMBER'].values.min() + + if end_prof is None: + end_prof = ds['PROFILE_NUMBER'].values.max() + + vert_curr = ds.VERT_CURR_MODEL.values * 100 # Convert to cm/s + vert_dzdt = ds.GLIDER_VERT_VELO_DZDT.values * 100 # Convert to cm/s + vert_model = ds.GLIDER_VERT_VELO_MODEL.values * 100 # Convert to cm/s + + + fig, axs = plt.subplots(2, 2, figsize=(14, 12), gridspec_kw={'width_ratios': [3, 1]}) + + # Upper left subplot for vertical velocity and glider speed + ax1 = axs[0, 0] + ax1.plot(ds['TIME'], vert_dzdt, label='Vertical Velocity (from dz/dt)') + ax1.set_xlabel('Time') + ax1.set_ylabel('Vertical Velocity (cm/s)') + ax1.legend(loc='lower left') + + ax1_twin = ax1.twinx() + ax1_twin.plot(ds['TIME'], vert_model, color='r', label='Vertical Glider Speed (model)') + ax1_twin.legend(loc='lower left') + + ax1_twin.plot(ds['TIME'], vert_curr, color='g', label='Vertical Water Speed (calculated)') + ax1_twin.legend(loc='lower right') + # Remove y-axis labels on the right + ax1_twin.set_yticklabels([]) + + # Upper right subplot for histogram of vertical velocity + ax1_hist = axs[0, 1] + ax1_hist.hist(vert_dzdt, bins=50, orientation='horizontal', alpha=0.5, color='blue', label='Vertical Velocity (from dz/dt)') + ax1_hist.hist(vert_model, bins=50, orientation='horizontal', alpha=0.5, color='red', label='Vertical Glider Speed (model)') + ax1_hist.hist(vert_curr, bins=50, orientation='horizontal', alpha=0.5, color='green', label='Vertical Water Speed (calculated)') + ax1_hist.set_xlabel('Frequency') + + # Determine the best location for the legend based on the y-axis limits and zero + y_upper_limit = ax1_hist.get_ylim()[1] + y_lower_limit = ax1_hist.get_ylim()[0] + + if abs(y_upper_limit) > abs(y_lower_limit): + legend_loc = 'upper right' + else: + legend_loc = 'lower right' + + ax1_hist.legend(loc=legend_loc) + + # Lower left subplot for vertical water speed + ax2 = axs[1, 0] + ax2.axhline(0, color='darkgray', linestyle='-', linewidth=0.5) # Add zero horizontal line + ax2.plot(ds['TIME'], vert_curr, 'g', label='Vertical Water Speed') + ax2.set_xlabel('Time') + ax2.set_ylabel('Vertical Water Speed (cm/s)') + ax2.legend(loc='upper left') + + # Lower right subplot for histogram of vertical water speed + ax2_hist = axs[1, 1] + ax2_hist.hist(vert_curr, bins=50, orientation='horizontal', alpha=0.5, color='green', label='Vertical Water Speed (calculated)') + ax2_hist.set_xlabel('Frequency') + + # Calculate and plot the median line + median_vert_sw_speed = np.nanmedian(vert_curr) + ax2_hist.axhline(median_vert_sw_speed, color='red', linestyle='dashed', linewidth=1, label=f'Median: {median_vert_sw_speed:.2f} cm/s') + + # Determine the best location for the legend based on the y-axis limits and median + y_upper_limit = ax2_hist.get_ylim()[1] + y_lower_limit = ax2_hist.get_ylim()[0] + median_vert_sw_speed = np.nanmedian(vert_curr) + + if abs(y_upper_limit - median_vert_sw_speed) > abs(y_lower_limit - median_vert_sw_speed): + legend_loc = 'upper right' + else: + legend_loc = 'lower right' + + ax2_hist.legend(loc=legend_loc) + + plt.tight_layout() + plt.show() + + return fig, axs + + diff --git a/notebooks/demo.ipynb b/notebooks/demo.ipynb index 427032a..68b074d 100644 --- a/notebooks/demo.ipynb +++ b/notebooks/demo.ipynb @@ -21,7 +21,7 @@ { "cell_type": "code", "execution_count": null, - "id": "fc10d805", + "id": "23dfc3fb", "metadata": {}, "outputs": [], "source": [ @@ -67,6 +67,32 @@ "cell_type": "markdown", "id": "6b93a2ef", "metadata": {}, + "source": [ + "### Vertical velocity" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02264a97", + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate vertical seawater velocity \n", + "# First, calculate the vertical speed of the glider from the depth data\n", + "ds = tools.calc_glider_w_from_depth(ds)\n", + "\n", + "# Next, calculate the vertical seawater speed by differencing the DZDT data and the modelled vertical glider speed\n", + "ds = tools.calc_seawater_w(ds)\n", + "\n", + "tools.plot_vertical_speeds_with_histograms(ds)\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "d22b9ca5", + "metadata": {}, "source": [ "# Basic statistics of dataset" ] @@ -101,7 +127,7 @@ "outputs": [], "source": [ "# Basic diagnostics of the watermass properties\n", - "tools.plot_ts_histograms(ds)" + "tools.plot_ts_histograms(ds)\n" ] }, {