Transect rewrite

docs/api/transect.rst

@@ -1,24 +1,64 @@
-.. module:: emsarray.transect
 
 =================
 emsarray.transect
 =================
 
-.. currentmodule:: emsarray.transect
+.. module:: emsarray.transect
 
-Plot transects through your dataset.
-Transects are vertical slices along some path through your dataset.
+This module provides methods for extracting and plotting data
+along transects through your datasets.
+A transect path is represented as a :class:`shapely.LineString`.
+Data along the transect can be extracted in to a new :class:`xarray.Dataset`,
+or plotted using :meth:`Transect.make_artist`.
+
+Currently it is only possible to take transects through grids with polygonal geometry.
+Taking transects through other kinds of geometry is a planned future enhancement.
 
 Examples
---------
+========
+
+.. minigallery::
 
-.. minigallery:: ../examples/plot-kgari-transect.py
+   ../examples/plot-kgari-transect.py
+   ../examples/plot-animated-transect.py
 
-.. autofunction:: plot
+Transects
+=========
+
+These classes find the intersection of a :class:`shapely.LineString` with a dataset
+and provide methods to introspect this intersection, plot data along this path,
+and extract data along this path.
 
 .. autoclass:: Transect
    :members:
 
-.. autoclass:: TransectPoint
+.. autoclass:: TransectPoint()
+   :members:
+
+.. autoclass:: TransectSegment()
+   :members:
+
+Artists
+=======
+
+These classes plot data along a transect.
+Transect artists are normally created by calling :meth:`.Transect.make_artist`.
+
+.. module:: emsarray.transect.artists
+
+.. autoclass:: TransectArtist()
+   :members: set_data_array
+
+.. autoclass:: CrossSectionArtist()
+   :members: from_transect
+
+.. autoclass:: TransectStepArtist()
+   :members: from_transect
+
+Utilities
+=========
+
+.. currentmodule:: emsarray.transect
 
-.. autoclass:: TransectSegment
+.. autofunction:: setup_distance_axis
+.. autofunction:: setup_depth_axis

docs/conf.py

@@ -111,6 +111,7 @@ def setup(app):
     'examples_dirs': '../examples',
     'gallery_dirs': './examples',
     'filename_pattern': '/plot-',
-    'matplotlib_animations': True,
+    'matplotlib_animations': (True, 'jshtml'),
     'backreferences_dir': './examples/backreferences',
+    'remove_config_comments': True,
 }

examples/plot-animated-transect.py

@@ -0,0 +1,110 @@
+"""
+=================
+Animated transect
+=================
+
+Transect and cross section plots can be animated using
+:meth:`.TransectArtist.set_data_array()` to update the data.
+"""
+
+import shapely
+import matplotlib.pyplot as plt
+import pandas
+import xarray
+from matplotlib.artist import Artist
+from matplotlib.animation import FuncAnimation
+from matplotlib.colors import LogNorm
+from matplotlib.ticker import ScalarFormatter
+
+from emsarray import transect, utils
+from emsarray.operations import depth
+
+# The eReefs GBR4 Biogeochemistry and Sediments v4.2 baseline catchment scenario dataset
+# contains one daily timestep per file.
+# For an animation we need to open multiple files using xarray.open_mfdataset().
+dataset_url_template = (
+    'https://thredds.nci.org.au/thredds/dodsC/fx3/gbr4_H4p0_ABARRAr2_OBRAN2020_FG2Gv3_B4p2_Cq5b_Dhnd/'
+    'gbr4_H4p0_ABARRAr2_OBRAN2020_FG2G_B4p2_Cq5b_Dhnd_simple_{date:%Y-%m-%d}.nc'
+)
+dataset_urls = [
+    dataset_url_template.format(date=date)
+    for date in pandas.date_range('2022-10-01', '2022-10-14')
+]
+dataset = xarray.open_mfdataset(
+    dataset_urls, combine='by_coords', coords=['time'],
+    compat='override', data_vars='minimal')
+
+# Select only the variables we want to plot.
+dataset = dataset.ems.select_variables(['botz', 'DIN'])
+
+# Fetch all the data.
+# This will only fetch the variable we selected above,
+# plus all the geometry and coordinate variables.
+dataset.load()
+
+# The depth coordinate has positive=up, while the bathymetry has positive=down.
+# This causes issues when drawing the ocean floor.
+# Lets fix the depth coordinate.
+dataset = depth.normalize_depth_variables(dataset, ['zc'], positive_down=True)
+
+# Cross section plots need bounds information, so lets invent some
+dataset = utils.estimate_bounds_1d(dataset, 'zc')
+
+# %%
+# Define a :class:`transect <emsarray.transect.Transect>`
+# following a path starting near Bowen, passing Mackay, and ending past Gladstone.
+reef_transect = transect.Transect(dataset, shapely.LineString([
+    [149.19358465107092, -19.745298160117585],
+    [150.16742824240822, -21.444631352206670],
+    [151.33964738012760, -22.175746836082112],
+    [152.21129750817641, -23.780650909462680]
+]))
+
+
+# %%
+# Now we set up a figure and add a cross section artist.
+
+# sphinx_gallery_defer_figures
+
+figure = plt.figure(figsize=(7.8, 3), layout='constrained', dpi=100)
+axes = figure.add_subplot()
+
+axes.set_title('Dissolved inorganic nitrogen')
+date_labels = [
+    utils.datetime_from_np_time(date).strftime('%Y-%m-%d')
+    for date in dataset['time'].values
+]
+
+din = dataset['DIN']
+din_artist = reef_transect.make_cross_section_artist(
+    axes, din.isel(time=0), cmap='plasma',
+    norm=LogNorm(0.01, .7), colorbar=False)
+figure.colorbar(
+    din_artist, ticks=[0.01, 0.025, 0.05, 0.1, 0.25, 0.5],
+    format=ScalarFormatter())
+
+reef_transect.make_ocean_floor_artist(
+    axes, dataset['botz'])
+
+date_annotation = axes.annotate(
+    date_labels[0],
+    xy=(5, 5), xycoords='axes points',
+    verticalalignment='bottom', horizontalalignment='left')
+
+transect.setup_distance_axis(
+    reef_transect, axes)
+transect.setup_depth_axis(
+    reef_transect, axes, data_array='DIN',
+    label='Depth', ylim=(80, -1.5))
+
+# %%
+# Finally we set up the animation.
+# The ``update()`` function is called every frame to update the plot with new data.
+# The :meth:`.TransectArtist.set_data_array()` function does all the hard work here.
+
+def update(frame: int) -> list[Artist]:
+    din_artist.set_data_array(din.isel(time=frame))
+    date_annotation.set_text(date_labels[frame])
+    return [din_artist, date_annotation]
+
+animation = FuncAnimation(figure, update, frames=din.sizes['time'])

examples/plot-kgari-transect.py

@@ -5,20 +5,31 @@
 """
 
 import shapely
-from matplotlib import pyplot
+import matplotlib.pyplot as plt
+from matplotlib import gridspec
+from matplotlib.colors import LogNorm, PowerNorm
 
 import emsarray
-from emsarray import plot, transect
+from emsarray import plot, transect, utils
+from emsarray.operations import depth
 
-dataset_url = 'https://thredds.nci.org.au/thredds/dodsC/fx3/gbr4_H4p0_ABARRAr2_OBRAN2020_FG2Gv3_Dhnd/gbr4_simple_2022-10-31.nc'
-dataset = emsarray.open_dataset(dataset_url).isel(time=-1)
-dataset = dataset.ems.select_variables(['botz', 'temp'])
+dataset_url = 'https://thredds.nci.org.au/thredds/dodsC/fx3/gbr4_H4p0_ABARRAr2_OBRAN2020_FG2Gv3_Dhnd/gbr4_simple_2022-10-01.nc'
+# dataset_url = '~/example-datasets/gbr4_simple_2022-10-31.nc'
+dataset = emsarray.open_dataset(dataset_url).isel(time=12)
+# Select only the variables we want to plot.
+dataset = dataset.ems.select_variables(['botz', 'temp', 'eta'])
+# The depth coordinate has positive=up, while the bathymetry has positive=down.
+# This causes issues when drawing the ocean floor.
+# Lets fix the depth coordinate.
+dataset = depth.normalize_depth_variables(dataset, ['zc'], positive_down=True)
+# Cross section plots need bounds information, so lets invent some
+dataset = utils.estimate_bounds_1d(dataset, 'zc')
 
 # %%
 # The following is a :mod:`transect <emsarray.transect>` path
 # starting in the Great Sandy Strait near K'gari,
 # heading roughly North out to deeper waters:
-line = shapely.LineString([
+north_transect = transect.Transect(dataset, shapely.LineString([
     [152.9768944, -25.4827962],
     [152.9701996, -25.4420345],
     [152.9727745, -25.3967620],
@@ -33,41 +44,94 @@
     [152.7607727, -24.3521012],
     [152.6392365, -24.1906056],
     [152.4792480, -24.0615124],
-])
+]))
 landmarks = [
     ('Round Island', shapely.Point(152.9262543, -25.2878719)),
     ('Lady Elliot Island', shapely.Point(152.7145958, -24.1129146)),
 ]
 
+
 # %%
-# Plot a transect showing temperature along this path.
+# Set up three axes: one showing the transect path,
+# one showing a temperature cross section along the transect,
+# and one showing the sea surface height along the transect.
+
+# sphinx_gallery_defer_figures
 
-figure = transect.plot(
-    dataset, line, dataset['temp'],
-    figsize=(7.9, 3),
-    bathymetry=dataset['botz'],
-    landmarks=landmarks,
-    title="Temperature",
-    cmap='Oranges_r')
-pyplot.show()
+figure = plt.figure(figsize=(7.8, 8), layout='constrained', dpi=100)
+gs_root = gridspec.GridSpec(3, 1, figure=figure, height_ratios=[3, 1, 1])
+path_axes = figure.add_subplot(gs_root[0], projection=dataset.ems.data_crs)
+temp_axes = figure.add_subplot(gs_root[1])
+eta_axes = figure.add_subplot(gs_root[2], sharex=temp_axes)
 
 # %%
-# The path of the transect can be plotted using matplotlib.
+# First make a plot showing the path of the transect overlayed on the bathymetry
+
+# sphinx_gallery_defer_figures
+# sphinx_gallery_capture_repr_block = ()
 
-# Plot the path of the transect
-figure = pyplot.figure(figsize=(5, 5), dpi=100)
-axes = figure.add_subplot(projection=dataset.ems.data_crs)
-axes.set_aspect(aspect='equal', adjustable='datalim')
-axes.set_title('Transect path')
+path_axes.set_aspect(aspect='equal', adjustable='datalim')
+path_axes.set_title('Transect path')
 dataset.ems.make_artist(
-    axes, 'botz', cmap='Blues', clim=(0, 2000), edgecolor='face',
+    path_axes, 'botz', cmap='Blues', clim=(0, 2000), edgecolor='face',
+    norm=PowerNorm(gamma=0.5),
     linewidth=0.5, zorder=0)
-axes = figure.axes[0]
-axes.set_extent(plot.bounds_to_extent(line.envelope.buffer(0.2).bounds))
-axes.plot(*line.coords.xy, zorder=2, c='orange', linewidth=4)
+path_axes.set_extent(plot.bounds_to_extent(north_transect.line.envelope.buffer(0.2).bounds))
+path_axes.plot(*north_transect.line.coords.xy, zorder=1, c='orange', linewidth=4)
+
+plot.add_coast(path_axes, zorder=1)
+plot.add_gridlines(path_axes)
+plot.add_landmarks(path_axes, landmarks)
+
+# %%
+# Now plot a cross section along the transect showing the ocean temperature.
+# As the temperature variable has a depth axis the cross section is two dimensional.
+
+# sphinx_gallery_defer_figures
+# sphinx_gallery_capture_repr_block = ()
+
+temp_axes.set_title('Temperature')
+dataset['temp'].attrs['units'] = '°C'
+dataset['zc'].attrs['long_name'] = 'Depth'
+
+north_transect.make_artist(
+    temp_axes, 'temp', cmap='plasma')
+north_transect.make_ocean_floor_artist(
+    temp_axes, dataset['botz'])
+# yaxis
+transect.setup_depth_axis(
+    north_transect, temp_axes, data_array='temp',
+    label='Depth', ylim=(50, -1.5))
+
+
+# %%
+# Now plot the sea surface height along the transect.
+# As the sea surface height does not have a depth axis
+# the transect is one dimensional.
+
+# sphinx_gallery_defer_figures
+# sphinx_gallery_capture_repr_block = ()
+
+eta_axes.set_title('Sea surface height')
+eta_artist = north_transect.make_artist(
+    eta_axes, data_array=dataset['eta'])
+# xaxis
+transect.setup_distance_axis(north_transect, eta_axes)
+# yaxis
+eta_axes.set_ylim(-0.5, 1.5)
+eta_axes.set_ylabel('Height above\nmean sea level')
+eta_axes.axhline(0, linestyle='--', color='lightgrey')
+eta_axes.yaxis.set_major_formatter("{x:.2g} m")
+
+# %%
+# The last step is to add some landmarks along the top border of the axes
+# to help viewers link the distance along transect path to geographic locations.
+
+top_axis = temp_axes.secondary_xaxis('top')
+top_axis.set_ticks(
+    [north_transect.distance_along_line(point) for label, point in landmarks],
+    [label for label, point in landmarks],
+)
 
-plot.add_coast(axes, zorder=1)
-plot.add_gridlines(axes)
-plot.add_landmarks(axes, landmarks)
 
-pyplot.show()
+plt.show()

src/emsarray/conventions/_base.py

@@ -251,9 +251,11 @@ def wind(
         axis : int, optional
             The axis number that should be wound.
             Optional, defaults to the last axis.
+            Mutually exclusive with the `linear_dimension` parameter.
         linear_dimension : Hashable, optional
-            The axis number that should be wound.
+            The name of the dimension in the data array that should be wound.
             Optional, defaults to the last dimension.
+            Mutually exclusive with the `axis` parameter.
 
         Returns
         -------
@@ -989,9 +991,11 @@ def wind(
         axis : int, optional
             The axis number that should be wound.
             Optional, defaults to the last axis.
+            Mutually exclusive with the `linear_dimension` parameter.
         linear_dimension : Hashable, optional
-            The axis number that should be wound.
+            The name of the dimension in the data array that should be wound.
             Optional, defaults to the last dimension.
+            Mutually exclusive with the `axis` parameter.
 
         Returns
         -------

src/emsarray/conventions/arakawa_c.py

@@ -372,7 +372,7 @@ def _make_geometry_centroid(self, grid_kind: ArakawaCGridKind) -> numpy.ndarray:
         return cast(numpy.ndarray, points)
 
     def get_all_geometry_names(self) -> list[Hashable]:
-        return [
+        return utils.geometry_plus_bounds(self.dataset, [
             self.face.longitude.name,
             self.face.latitude.name,
             self.node.longitude.name,
@@ -381,7 +381,7 @@ def get_all_geometry_names(self) -> list[Hashable]:
             self.left.latitude.name,
             self.back.longitude.name,
             self.back.latitude.name,
-        ]
+        ])
 
     def make_clip_mask(
         self,