diff --git a/glidertest/tools.py b/glidertest/tools.py index 1a854cb..96dd17f 100644 --- a/glidertest/tools.py +++ b/glidertest/tools.py @@ -10,7 +10,10 @@ from skyfield import almanac from skyfield import api from tqdm import tqdm - +import matplotlib.colors as mcolors +import gsw +import cartopy.crs as ccrs +import cartopy.feature as cfeature def grid2d(x, y, v, xi=1, yi=1): """ @@ -588,3 +591,211 @@ def plot_profIncrease(ds: xr.DataArray, ax: plt.Axes = None, **kw: dict, ) -> tu ax[1].xaxis.set_major_locator(plt.MaxNLocator(8)) [a.grid() for a in ax] return fig, ax + + +def plot_glider_track(ds: xr.Dataset, ax: plt.Axes = None, **kw: dict) -> tuple({plt.Figure, plt.Axes}): + """ + This function plots the glider track on a map, with latitude and longitude colored by time. + + Parameters + ---------- + ds: xarray in OG1 format with at least TIME, LATITUDE, and LONGITUDE. + ax: Optional; axis to plot the data. + kw: Optional; additional keyword arguments for the scatter plot. + + Returns + ------- + One plot with the map of the glider track. + fig: matplotlib.figure.Figure + ax: matplotlib.axes._subplots.AxesSubplot + """ + if ax is None: + fig, ax = plt.subplots(figsize=(10, 6), subplot_kw={'projection': ccrs.PlateCarree()}) + else: + fig = plt.gcf() + + latitudes = ds.LATITUDE.values + longitudes = ds.LONGITUDE.values + times = ds.TIME.values + + # Reduce the number of latitudes, longitudes, and times to no more than 2000 + if len(latitudes) > 2000: + indices = np.linspace(0, len(latitudes) - 1, 2000).astype(int) + latitudes = latitudes[indices] + longitudes = longitudes[indices] + times = times[indices] + + # Convert time to the desired format + time_labels = [pd.to_datetime(t).strftime('%Y-%b-%d') for t in times] + + # Plot latitude and longitude colored by time + sc = ax.scatter(longitudes, latitudes, c=times, cmap='viridis', **kw) + + # Add colorbar with formatted time labels + cbar = plt.colorbar(sc, ax=ax) #, label='Time') + cbar.ax.set_yticklabels([pd.to_datetime(t).strftime('%Y-%b-%d') for t in cbar.get_ticks()]) + + ax.set_extent([np.min(longitudes)-1, np.max(longitudes)+1, np.min(latitudes)-1, np.max(latitudes)+1], crs=ccrs.PlateCarree()) + + # Add features to the map + ax.add_feature(cfeature.LAND) + ax.add_feature(cfeature.OCEAN) + ax.add_feature(cfeature.COASTLINE) + + ax.set_xlabel('Longitude') + ax.set_ylabel('Latitude') + ax.set_title('Glider Track') + gl = ax.gridlines(draw_labels=True, color='black', alpha=0.5, linestyle='--') + gl.top_labels = False + gl.right_labels = False + plt.show() + + return fig, ax + +def plot_grid_spacing_histograms(ds: xr.Dataset, ax: plt.Axes = None, **kw: dict) -> tuple({plt.Figure, plt.Axes}): + """ + This function plots histograms of the grid spacing (diff(ds.DEPTH) and diff(ds.TIME)) where only the inner 99% of values are plotted. + + Parameters + ---------- + ds: xarray in OG1 format with at least TIME and DEPTH. + ax: Optional; axis to plot the data. + kw: Optional; additional keyword arguments for the histograms. + + Returns + ------- + Two histograms showing the distribution of grid spacing for depth and time. + fig: matplotlib.figure.Figure + ax: matplotlib.axes._subplots.AxesSubplot + """ + if ax is None: + fig, ax = plt.subplots(1, 2, figsize=(14, 6)) + else: + fig = plt.gcf() + + depth_diff = np.diff(ds.DEPTH) + orig_time_diff = np.diff(ds.TIME) / np.timedelta64(1, 's') # Convert to seconds + + depth_diff = depth_diff[np.isfinite(depth_diff)] + time_diff = orig_time_diff[np.isfinite(orig_time_diff)] + + depth_diff = depth_diff[(depth_diff >= np.nanpercentile(depth_diff, 0.5)) & (depth_diff <= np.nanpercentile(depth_diff, 99.5))] + time_diff = time_diff[(time_diff >= np.nanpercentile(time_diff, 0.5)) & (time_diff <= np.nanpercentile(time_diff, 99.5))] + + ax[0].hist(depth_diff, bins=50, **kw) + ax[0].set_xlabel('Depth Spacing (m)') + ax[0].set_ylabel('Frequency') + ax[0].set_title('Histogram of Depth Spacing') + ax[0].grid() + + median_neg_depth_diff = np.median(depth_diff[depth_diff < 0]) + median_pos_depth_diff = np.median(depth_diff[depth_diff > 0]) + max_depth_diff = np.max(depth_diff) + min_depth_diff = np.min(depth_diff) + + annotation_text_left = ( + f'Median Negative: {median_neg_depth_diff:.2f} m\n' + f'Median Positive: {median_pos_depth_diff:.2f} m\n' + f'Max: {max_depth_diff:.2f} m\n' + f'Min: {min_depth_diff:.2f} m' + ) + ax[0].annotate(annotation_text_left, xy=(0.96, 0.96), xycoords='axes fraction', + fontsize=12, ha='right', va='top', + bbox=dict(boxstyle='round,pad=0.3', edgecolor='black', facecolor='white', alpha=.5)) + + ax[1].hist(time_diff, bins=50, **kw) + ax[1].set_xlabel('Time Spacing (s)') + ax[1].set_ylabel('Frequency') + ax[1].set_title('Histogram of Time Spacing') + ax[1].grid() + + median_time_diff = np.median(orig_time_diff) + mean_time_diff = np.mean(orig_time_diff) + max_time_diff = np.max(orig_time_diff) + min_time_diff = np.min(orig_time_diff) + max_time_diff_hrs = max_time_diff/3600 + + annotation_text = ( + f'Median: {median_time_diff:.2f} s\n' + f'Mean: {mean_time_diff:.2f} s\n' + f'Max: {max_time_diff:.2f} s ({max_time_diff_hrs:.2f} hr)\n' + f'Min: {min_time_diff:.2f} s' + ) + ax[1].annotate(annotation_text, xy=(0.96, 0.96), xycoords='axes fraction', + fontsize=12, ha='right', va='top', + bbox=dict(boxstyle='round,pad=0.3', edgecolor='black', facecolor='white', alpha=.5)) + + return fig, ax + +def plot_ts_histograms(ds: xr.Dataset, ax: plt.Axes = None, **kw: dict) -> tuple({plt.Figure, plt.Axes}): + """ + This function plots histograms of temperature and salinity values (middle 95%), and a 2D histogram of salinity and temperature with density contours. + + Parameters + ---------- + ds: xarray in OG1 format with at least TEMP and PSAL. + ax: Optional; axis to plot the data. + kw: Optional; additional keyword arguments for the histograms. + + Returns + ------- + Three plots: histogram of temperature, histogram of salinity, and 2D histogram of salinity and temperature with density contours. + fig: matplotlib.figure.Figure + ax: matplotlib.axes._subplots.AxesSubplot + """ + + if ax is None: + fig, ax = plt.subplots(1, 3, figsize=(18, 6)) + else: + fig = plt.gcf() + + temp_orig = ds.TEMP.values + sal_orig = ds.PSAL.values + + # Reduce both to where both are finite + temp = temp_orig[np.isfinite(temp_orig)&np.isfinite(sal_orig)] + sal = sal_orig[np.isfinite(sal_orig)&np.isfinite(temp_orig)] + depth = ds.DEPTH[np.isfinite(sal_orig)&np.isfinite(temp_orig)] + long = ds.LONGITUDE[np.isfinite(sal_orig)&np.isfinite(temp_orig)] + lat = ds.LATITUDE[np.isfinite(sal_orig)&np.isfinite(temp_orig)] + + SA = gsw.SA_from_SP(sal, depth, long, lat) + 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))] + + ax[0].hist(temp_filtered, bins=50, **kw) + ax[0].set_xlabel('Conservative Temperature (°C)') + ax[0].set_ylabel('Frequency') + ax[0].set_title('Histogram of Temperature') + ax[0].grid() + + ax[1].hist(sal_filtered, bins=50, **kw) + ax[1].set_xlabel('Absolute Salinity ( )') + ax[1].set_ylabel('Frequency') + ax[1].set_title('Histogram of Salinity') + ax[1].grid() + + h = ax[2].hist2d(sal, temp, bins=50, cmap='viridis', norm=mcolors.LogNorm(), **kw) + cbar = plt.colorbar(h[3], ax=ax[2]) + cbar.set_label('Log Counts') + ax[2].set_xlabel('Absolute Salinity ( )') + ax[2].set_ylabel('Conservative Temperature (°C)') + ax[2].set_title('2D Histogram of Salinity and Temperature (Log Scale)') + + # Calculate density and add contours + SA = gsw.SA_from_SP(sal, depth, long, lat) + CT = gsw.CT_from_t(SA, temp, depth) + sigma0 = gsw.sigma0(SA, CT) + + xi = np.linspace(sal.min()-1, sal.max()+1, 100) + yi = np.linspace(temp.min()-1, temp.max()+1, 100) + xi, yi = np.meshgrid(xi, yi) + zi = gsw.sigma0(xi, yi) + + 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 diff --git a/notebooks/demo.ipynb b/notebooks/demo.ipynb index 1cf834b..b0493d8 100644 --- a/notebooks/demo.ipynb +++ b/notebooks/demo.ipynb @@ -60,7 +60,48 @@ "metadata": {}, "outputs": [], "source": [ - "ds" + "ds\n" + ] + }, + { + "cell_type": "markdown", + "id": "6b93a2ef", + "metadata": {}, + "source": [ + "# Basic statistics of dataset" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c6fd094f", + "metadata": {}, + "outputs": [], + "source": [ + "# Basic plot of the location of the dataset in space/time\n", + "tools.plot_glider_track(ds)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "620e5523", + "metadata": {}, + "outputs": [], + "source": [ + "# Basic diagnostics of the gridding in the dataset\n", + "tools.plot_grid_spacing_histograms(ds)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6bf7fb8d", + "metadata": {}, + "outputs": [], + "source": [ + "# Basic diagnostics of the watermass properties\n", + "tools.plot_ts_histograms(ds)\n" ] }, { diff --git a/requirements-dev.txt b/requirements-dev.txt index 78e2856..22ed076 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -16,3 +16,4 @@ nbsphinx nbconvert sphinx sphinx-rtd-theme +cartopy \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index 52ba604..dc66e0c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -11,3 +11,4 @@ skyfield tqdm gsw sgp4>=2.14 +cartopy \ No newline at end of file