From 5fed1e47bb83369cdfc7ad6810ed72aab721ac95 Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Tue, 13 Oct 2026 17:59:12 -0500 Subject: [PATCH 1/5] Project contour lines once, instead of letting GeoViews project each one With a projection, almost all the time of a contour plot went into GeoViews projecting the lines one at a time through Cartopy while rendering: ~13 of the ~17 s for a noisy field on a 786k-face mesh. GeoViews does this for any element that has a crs, even one already in the plot's projection, and a plain HoloViews element would lose GeoViews' map extents and options such as global_extent. _project_lines replaces _drop_lines_outside_projection. It projects all points and segment midpoints at once and applies Cartopy's own test, from cartopy/trace.pyx, to every segment: the projected midpoint within the projection's threshold of the segment, every point on the map, and the segment covered by the map when the map is not convex. Lines that pass are exactly what Cartopy gives; only the others, such as those that cross the edge of the map or a seam of the projection, are projected by Cartopy. The result matched Cartopy line by line, with no difference in any coordinate, for 10 projections on four datasets and both methods. The lines are a subclass of gv.Contours with crs=projection, plotted by GeoViews' own GeoContourPlot with its projection step turned off. It is named Contours to share the options of contour lines, created on first use since GeoViews is optional, and found by pickle through a module __getattr__. Noisy field, 786k faces, contour() + render: Robinson, Bokeh 17.1 s -> 2.3 s Robinson, Matplotlib 17.1 s -> 2.3 s Orthographic, Bokeh 16.1 s -> 5.7 s GeoViews now frames a plot of the lines alone on the projected lines, rather than on the projection of their longitude-latitude bounding box, so the extent of some maps (Lambert conformal, orthographic) changes slightly; the lines themselves are the same. Co-Authored-By: Claude Opus 5.5 --- test/test_plot.py | 36 +++++++--- uxarray/plot/accessor.py | 34 ++++++--- uxarray/plot/contour.py | 151 +++++++++++++++++++++++++++++++++------ 3 files changed, 178 insertions(+), 43 deletions(-) diff --git a/test/test_plot.py b/test/test_plot.py index 89c804a41..4614d696a 100644 --- a/test/test_plot.py +++ b/test/test_plot.py @@ -638,6 +638,16 @@ def test_contour_projection(gridpath): assert isinstance(lines, gv.Contours) and isinstance(labels, gv.Labels) hv.renderer("matplotlib").get_plot(uxda.plot(projection=projection, backend="matplotlib") * contours) + # the lines are already in the coordinates of the projection, so GeoViews plots them + # as a map without projecting them again, on both backends + x_limit = projection.x_limits[1] + assert lines.crs == projection + assert x_limit / 2 < max(np.abs(x).max() for _, x, _ in _contour_lines(lines)) <= x_limit + for backend in ["matplotlib", "bokeh"]: + lines = uxda.plot.contour(levels=levels, projection=projection, backend=backend, labels=False) + plot = hv.renderer(backend).get_plot(lines) + assert plot.geographic and plot._project_operation is None + # seen from above the north pole, the southern line is not on the map. # GeoViews raises an IndexError for such a line, so it is dropped for method in ["interpolated", "edges"]: @@ -650,21 +660,27 @@ def test_contour_projection(gridpath): @pytest.mark.parametrize("projection", [ ccrs.Orthographic(0, 90), ccrs.PlateCarree(central_longitude=180), ccrs.LambertConformal(), ccrs.Robinson(-150) ]) -def test_contour_projection_drops_same_lines_as_projecting_each(gridpath, projection): - """Lines found on the map from their points alone are the lines whose projection is not empty.""" +def test_contour_projected_lines_match_cartopy(gridpath, projection): + """Lines projected all at once are the lines Cartopy gives when it projects them one at a time.""" import shapely - from uxarray.plot.contour import _compute_contours, _drop_lines_outside_projection + from uxarray.plot.contour import _compute_contours, _project_lines uxda = _face_latitude(ux.open_grid(gridpath("ugrid", "outCSne30", "outCSne30.ug"))) for method in ["interpolated", "edges"]: contours = _compute_contours(uxda, levels=[-45.5, 0.5, 10.5, 45.5], method=method) - expected = [ - line for _, line in contours - if not projection.project_geometry(shapely.LineString(line), ccrs.PlateCarree()).is_empty - ] - kept = [line for _, line in _drop_lines_outside_projection(contours, projection)] - assert len(kept) == len(expected) - assert all(a is b for a, b in zip(kept, expected)) + expected, expected_on_map = [], [] + for level, line in contours: + geometry = projection.project_geometry(shapely.LineString(line), ccrs.PlateCarree()) + parts = [(level, np.asarray(part.coords)) for part in shapely.get_parts(geometry) if len(part.coords) > 1] + expected += parts + expected_on_map.append(len(parts) > 0) + + projected, on_map = _project_lines(contours, projection) + assert on_map.tolist() == expected_on_map + assert len(projected) == len(expected) + for (level, line), (expected_level, expected_line) in zip(projected, expected): + assert level == expected_level + np.testing.assert_allclose(line, expected_line, rtol=1e-12) def test_contour_frame(gridpath): diff --git a/uxarray/plot/accessor.py b/uxarray/plot/accessor.py index 9a87107bd..1a1be86e2 100644 --- a/uxarray/plot/accessor.py +++ b/uxarray/plot/accessor.py @@ -596,9 +596,10 @@ def contour( ------- contours : hv.Overlay, hv.Contours or gv.Contours The contour lines, one path per level with the level as its value - and its lines separated by NaN, overlaid with their labels. With ``labels=False`` the lines - alone: an ``hv.Contours``, or a ``gv.Contours`` if a ``projection`` is - given. + and its lines separated by NaN, overlaid with their labels. With + ``labels=False`` the lines alone: an ``hv.Contours``, or a + ``gv.Contours`` if a ``projection`` is given. Those lines are already + in the coordinates of the projection, which is their ``crs``. Raises ------ @@ -646,20 +647,23 @@ def contour( from uxarray.plot.contour import ( _apply_options, _compute_contours, - _drop_lines_outside_projection, _grid_extents, _label_fraction, _label_points, _paths_by_level, + _project_lines, + _projected_contours, ) label_fraction = _label_fraction(labels) plotting_backend.assign(backend) contours = _compute_contours(self._uxda, levels=levels, method=method) + labeled = contours if projection is not None: _raise_hint_if_optional_deps_missing("cartopy", "geoviews") - contours = _drop_lines_outside_projection(contours, projection) + contours, on_map = _project_lines(contours, projection) + labeled = [line for line, shown in zip(labeled, on_map) if shown] # The level has a fixed name, so data named "x" or "y" cannot replace a # coordinate; the name of the data is what hover and the colorbar show. @@ -672,19 +676,27 @@ def contour( hv.Dimension("y", label="Latitude"), ] paths = _paths_by_level(contours, level_dim.name) - label_x, label_y, label_text = _label_points(contours, label_fraction) + label_x, label_y, label_text = _label_points(labeled, label_fraction) label_data = {"x": label_x, "y": label_y, "text": label_text} extents = _grid_extents(self._uxda.uxgrid) if projection is not None: import cartopy.crs as ccrs import geoviews as gv - - crs = ccrs.PlateCarree() - lines = gv.Contours( - paths, kdims=kdims, vdims=[level_dim], crs=crs, extents=extents + from geoviews.util import project_extents + + # The lines are already projected, so the area of the grid is projected + # as GeoViews projects the extents of an element in longitude and latitude + lines = _projected_contours()( + paths, + kdims=kdims, + vdims=[level_dim], + crs=projection, + extents=project_extents(extents, ccrs.PlateCarree(), projection), + ) + text = gv.Labels( + label_data, kdims=kdims, vdims=["text"], crs=ccrs.PlateCarree() ) - text = gv.Labels(label_data, kdims=kdims, vdims=["text"], crs=crs) else: lines = hv.Contours(paths, kdims=kdims, vdims=[level_dim], extents=extents) text = hv.Labels(label_data, kdims=kdims, vdims=["text"]) diff --git a/uxarray/plot/contour.py b/uxarray/plot/contour.py index 4e5d40de4..c7e600a75 100644 --- a/uxarray/plot/contour.py +++ b/uxarray/plot/contour.py @@ -526,38 +526,145 @@ def _grid_extents(uxgrid: Grid) -> tuple[float, float, float, float]: return tuple(float(np.float32(value)) for value in (west, south, east, north)) -def _drop_lines_outside_projection( +def _project_lines( contours: list[tuple[float, np.ndarray]], projection -) -> list[tuple[float, np.ndarray]]: - """Removes the lines that have no extent in a map projection, such as a - line that lies along the edge of the map. GeoViews cannot project them.""" +) -> tuple[list[tuple[float, np.ndarray]], np.ndarray]: + """Projects contour lines from longitude and latitude to a map projection. + + Cartopy projects a line segment by segment, keeping a segment as it is if it + is inside the map and its projection is close to straight, and otherwise + splitting it or cutting it where it leaves the map. That is slow for many + lines, so the lines whose segments would all be kept are found by applying + the same test to all their points at once, and only the other lines, such as + those that cross the edge of the map or a seam of the projection, are + projected by Cartopy. + + Returns + ------- + projected : list of (level, line) + The parts of the lines on the map, in the coordinates of the projection. + on_map : np.ndarray + Whether each of the lines has a part on the map. + """ _raise_hint_if_optional_deps_missing("cartopy") import cartopy.crs as ccrs import shapely if not contours: - return contours + return [], np.zeros(0, dtype=bool) source = ccrs.PlateCarree() - # A line with a point well inside the map has extent on it. Projecting all the - # points at once is much faster than projecting each line, so only the other - # lines, such as those along the edge of the map, are projected as lines. - points = np.concatenate([line for _, line in contours]) - x, y = projection.transform_points(source, points[:, 0], points[:, 1])[:, :2].T - finite = np.isfinite(x) & np.isfinite(y) + lines = [line for _, line in contours] + counts = np.array([len(line) for line in lines]) + ends = np.cumsum(counts) + starts = ends - counts + points = np.concatenate(lines) + xy = projection.transform_points(source, points[:, 0], points[:, 1])[:, :2] + + # Every point on the map, which for Cartopy includes its edge domain = projection.domain - margin = 1e-6 * np.ptp(np.reshape(domain.bounds, (2, 2)), axis=0).max() - inside = np.zeros(len(points), dtype=bool) - inside[finite] = shapely.contains_xy(domain.buffer(-margin), x[finite], y[finite]) - starts = np.cumsum([0] + [len(line) for _, line in contours[:-1]]) - on_map = np.logical_or.reduceat(inside, starts) + finite = np.isfinite(xy).all(axis=1) + keep = np.zeros(len(points), dtype=bool) + keep[finite] = shapely.intersects_xy(domain, xy[finite, 0], xy[finite, 1]) + + # Every segment close to straight once projected, by Cartopy's test: the + # projection of its middle (Cartopy interpolates linearly in longitude and + # latitude) lies within the projection's threshold of the projected segment, + # less so towards its ends + middle = (points[:-1] + points[1:]) / 2 + xy_middle = projection.transform_points(source, middle[:, 0], middle[:, 1])[:, :2] + segment = xy[1:] - xy[:-1] + to_middle = xy_middle - xy[:-1] + with np.errstate(divide="ignore", invalid="ignore"): + length_sq = (segment**2).sum(axis=1) + along = (segment * to_middle).sum(axis=1) / length_sq + separation = np.abs( + to_middle[:, 0] * segment[:, 1] - to_middle[:, 1] * segment[:, 0] + ) / np.sqrt(length_sq) + straight = np.isnan(along) | ( + (along > 0.0) + & (along < 1.0) + & (separation <= projection.threshold * 2.0 * (0.5 - np.abs(0.5 - along))) + ) + # The segments between the end of a line and the start of the next do not count + straight[ends[:-1] - 1] = True + keep[:-1] &= straight + + # On a map that is not convex, a segment between two points on it can leave it + if not domain.equals(domain.convex_hull): + covered = np.ones(len(points) - 1, dtype=bool) + within_line = keep[:-1] & keep[1:] + within_line[ends[:-1] - 1] = False + pairs = np.stack([xy[:-1][within_line], xy[1:][within_line]], axis=1) + covered[within_line] = shapely.covers(domain, shapely.linestrings(pairs)) + keep[:-1] &= covered + + as_is = np.logical_and.reduceat(keep, starts) + + projected = [] + on_map = as_is.copy() + for i, (level, line) in enumerate(contours): + if as_is[i]: + projected.append((level, xy[starts[i] : ends[i]])) + continue + geometry = projection.project_geometry(shapely.LineString(line), source) + for part in shapely.get_parts(geometry): + if len(part.coords) > 1: + projected.append((level, np.asarray(part.coords)[:, :2])) + on_map[i] = True + return projected, on_map - return [ - (level, line) - for (level, line), certain in zip(contours, on_map) - if certain - or not projection.project_geometry(shapely.LineString(line), source).is_empty - ] + +_projected_contours_class = None + + +def _projected_contours(): + """The GeoViews element for contour lines that are already projected. + + GeoViews projects the lines of a ``Contours`` element one at a time, through + Cartopy, every time it plots them, even when they are already in the + coordinates of the plot. For the thousands of lines of a noisy field that is + most of the time of a plot. This subclass is plotted by GeoViews' own plot + classes with that step turned off, so map extents, options such as + ``global_extent`` and overlays behave as for any GeoViews element. It is + named ``Contours`` to share the options of other contour lines, and created + on first use since GeoViews is optional. + """ + global _projected_contours_class + if _projected_contours_class is None: + _raise_hint_if_optional_deps_missing("geoviews", "holoviews") + import geoviews as gv + import geoviews.plotting.bokeh as gv_bokeh + import geoviews.plotting.mpl as gv_mpl + import holoviews as hv + + contours = type( + "Contours", + (gv.Contours,), + { + "__module__": __name__, + "__qualname__": "ProjectedContours", + "__doc__": "GeoViews contour lines already in the coordinates of their crs.", + }, + ) + # Importing the GeoViews plotting modules registers both backends + # without changing the active one + for backend, plot in [ + ("bokeh", gv_bokeh.GeoContourPlot), + ("matplotlib", gv_mpl.GeoContourPlot), + ]: + hv.Store.registry[backend][contours] = type( + f"Projected{plot.__name__}", (plot,), {"_project_operation": None} + ) + _projected_contours_class = contours + return _projected_contours_class + + +def __getattr__(name): + # Lets pickle find the class of projected contour lines, which is created on first use + if name == "ProjectedContours": + return _projected_contours() + raise AttributeError(f"module {__name__!r} has no attribute {name!r}") def _paths_by_level( From 252cfabc53e5f465046af7d3ea2af077ec2c6a79 Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Tue, 13 Oct 2026 18:04:57 -0500 Subject: [PATCH 2/5] Credit Cartopy for the segment test in _project_lines The test that decides whether a projected segment can be kept as it is was adapted from straightAndDomain() in lib/cartopy/trace.pyx of Cartopy 0.25.0: the position along the segment and the separation of the projected midpoint from it, the tolerance that narrows towards the ends of the segment, and the requirement that the map cover the segment. Cartopy is released under the BSD 3-Clause License, which requires that redistributed source keep its copyright notice, conditions and disclaimer. The notice is reproduced verbatim above _project_lines, where it travels with the source into every distribution, and the docstring and comments name where the test comes from. Co-Authored-By: Claude Opus 5.5 --- uxarray/plot/contour.py | 51 +++++++++++++++++++++++++++++++++++++---- 1 file changed, 46 insertions(+), 5 deletions(-) diff --git a/uxarray/plot/contour.py b/uxarray/plot/contour.py index c7e600a75..a83db1895 100644 --- a/uxarray/plot/contour.py +++ b/uxarray/plot/contour.py @@ -526,6 +526,38 @@ def _grid_extents(uxgrid: Grid) -> tuple[float, float, float, float]: return tuple(float(np.float32(value)) for value in (west, south, east, north)) +# The segment test in _project_lines() is adapted from straightAndDomain() in +# lib/cartopy/trace.pyx of Cartopy 0.25.0, which is released under this license: +# +# BSD 3-Clause License +# +# Copyright (c) Crown and Cartopy contributors. +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, this +# list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. def _project_lines( contours: list[tuple[float, np.ndarray]], projection ) -> tuple[list[tuple[float, np.ndarray]], np.ndarray]: @@ -545,6 +577,14 @@ def _project_lines( The parts of the lines on the map, in the coordinates of the projection. on_map : np.ndarray Whether each of the lines has a part on the map. + + Notes + ----- + The test for a segment is adapted from ``straightAndDomain()`` in + ``lib/cartopy/trace.pyx`` of Cartopy 0.25.0, for a segment that starts on + the map, and applied here to all segments at once. Cartopy is Copyright (c) + Crown and Cartopy contributors, and is released under the BSD 3-Clause + License, whose notice is above this function. """ _raise_hint_if_optional_deps_missing("cartopy") import cartopy.crs as ccrs @@ -567,10 +607,10 @@ def _project_lines( keep = np.zeros(len(points), dtype=bool) keep[finite] = shapely.intersects_xy(domain, xy[finite, 0], xy[finite, 1]) - # Every segment close to straight once projected, by Cartopy's test: the - # projection of its middle (Cartopy interpolates linearly in longitude and - # latitude) lies within the projection's threshold of the projected segment, - # less so towards its ends + # Every segment close to straight once projected, by the test adapted from + # Cartopy's straightAndDomain() (see Notes): the projection of its middle + # (Cartopy interpolates linearly in longitude and latitude) lies within the + # projection's threshold of the projected segment, less so towards its ends middle = (points[:-1] + points[1:]) / 2 xy_middle = projection.transform_points(source, middle[:, 0], middle[:, 1])[:, :2] segment = xy[1:] - xy[:-1] @@ -590,7 +630,8 @@ def _project_lines( straight[ends[:-1] - 1] = True keep[:-1] &= straight - # On a map that is not convex, a segment between two points on it can leave it + # On a map that is not convex, a segment between two points on it can leave + # it, so Cartopy also requires the map to cover the segment if not domain.equals(domain.convex_hull): covered = np.ones(len(points) - 1, dtype=bool) within_line = keep[:-1] & keep[1:] From 3d4d19eaa36b31109d916dfc2d370e4ee82e84d2 Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Tue, 13 Oct 2026 18:10:17 -0500 Subject: [PATCH 3/5] Drop contour lines that cannot be projected without asking Cartopy On an orthographic map about half the lines of a field are on the far side of the globe, and none of their points can be projected. They were still given to Cartopy one at a time, only to come back empty. For a segment whose ends both cannot be projected, Cartopy's bisection in lib/cartopy/trace.pyx accepts the end at once and _project_segment() adds nothing, without looking at the points in between. A line with no point that can be projected is therefore never on the map, and is now dropped from the points already projected. The result still matches Cartopy line by line for 10 projections on four datasets and both methods, and projecting the lines of a 197k-face noisy field on an orthographic map takes 0.85 s instead of 1.9 s. The rule is credited to Cartopy next to the segment test, under the same BSD 3-Clause notice. Co-Authored-By: Claude Opus 5.5 --- test/test_plot.py | 20 ++++++++++++++++++++ uxarray/plot/contour.py | 22 +++++++++++++++++----- 2 files changed, 37 insertions(+), 5 deletions(-) diff --git a/test/test_plot.py b/test/test_plot.py index 4614d696a..963396e6c 100644 --- a/test/test_plot.py +++ b/test/test_plot.py @@ -3,6 +3,7 @@ import holoviews as hv import pytest import numpy as np +import shapely import matplotlib matplotlib.use("Agg") @@ -683,6 +684,25 @@ def test_contour_projected_lines_match_cartopy(gridpath, projection): np.testing.assert_allclose(line, expected_line, rtol=1e-12) +def test_contour_lines_on_far_side_are_dropped_without_cartopy(): + """A line none of whose points can be projected is not on the map, as Cartopy decides, without asking Cartopy.""" + from uxarray.plot.contour import _project_lines + + projection = ccrs.Orthographic(0, 90) + far_side = np.column_stack([np.linspace(-170, 170, 50), np.full(50, -45.5)]) + near_side = np.column_stack([np.linspace(-170, 170, 50), np.full(50, 45.5)]) + assert not np.isfinite(projection.transform_points(ccrs.PlateCarree(), *far_side.T)).all(axis=1).any() + assert projection.project_geometry(shapely.LineString(far_side), ccrs.PlateCarree()).is_empty + + def no_cartopy(*args, **kwargs): + raise AssertionError("the line was given to Cartopy") + + projection.project_geometry = no_cartopy + projected, on_map = _project_lines([(-45.5, far_side), (45.5, near_side)], projection) + assert on_map.tolist() == [False, True] + assert [level for level, _ in projected] == [45.5] + + def test_contour_frame(gridpath): """Lines alone get the frame of .plot(): equal axis scaling on a map, 700 x 300 without a projection. Overlaid on a shaded plot, the shaded plot's frame and area are kept.""" diff --git a/uxarray/plot/contour.py b/uxarray/plot/contour.py index a83db1895..ee57f6020 100644 --- a/uxarray/plot/contour.py +++ b/uxarray/plot/contour.py @@ -526,8 +526,10 @@ def _grid_extents(uxgrid: Grid) -> tuple[float, float, float, float]: return tuple(float(np.float32(value)) for value in (west, south, east, north)) -# The segment test in _project_lines() is adapted from straightAndDomain() in -# lib/cartopy/trace.pyx of Cartopy 0.25.0, which is released under this license: +# The segment test in _project_lines() is adapted from straightAndDomain(), and +# its rule for lines with no point that can be projected from bisect() and +# _project_segment(), in lib/cartopy/trace.pyx of Cartopy 0.25.0, which is +# released under this license: # # BSD 3-Clause License # @@ -582,9 +584,11 @@ def _project_lines( ----- The test for a segment is adapted from ``straightAndDomain()`` in ``lib/cartopy/trace.pyx`` of Cartopy 0.25.0, for a segment that starts on - the map, and applied here to all segments at once. Cartopy is Copyright (c) - Crown and Cartopy contributors, and is released under the BSD 3-Clause - License, whose notice is above this function. + the map, and applied here to all segments at once. The rule that a line + with no point that can be projected is not on the map follows ``bisect()`` + and ``_project_segment()`` in the same file. Cartopy is Copyright (c) Crown + and Cartopy contributors, and is released under the BSD 3-Clause License, + whose notice is above this function. """ _raise_hint_if_optional_deps_missing("cartopy") import cartopy.crs as ccrs @@ -642,12 +646,20 @@ def _project_lines( as_is = np.logical_and.reduceat(keep, starts) + # A line none of whose points can be projected, such as one on the far side + # of an orthographic map, is not on the map. For a segment whose ends both + # cannot be projected, Cartopy's bisection accepts the end at once and adds + # nothing, without looking at the points in between (see Notes) + off_map = ~np.logical_or.reduceat(finite, starts) + projected = [] on_map = as_is.copy() for i, (level, line) in enumerate(contours): if as_is[i]: projected.append((level, xy[starts[i] : ends[i]])) continue + if off_map[i]: + continue geometry = projection.project_geometry(shapely.LineString(line), source) for part in shapely.get_parts(geometry): if len(part.coords) > 1: From 580f42a779684e5853d24274c1baf2b449f8940a Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Tue, 13 Oct 2026 19:03:32 -0500 Subject: [PATCH 4/5] Contour the interpolated method with a numba kernel instead of tricontour The "interpolated" method handed its triangulation to Matplotlib's tricontour, which spent ~2.2 s of its ~2.6 s on a 786k-face mesh setting up the neighbours of every triangle before contouring. _crossed_edges, compiled with numba (cache=True, nogil=True), finds the two edges where a level crosses each triangle that it crosses. Each crossing point is known by its edge, which the triangles on either side share, so the pieces are joined into lines by the existing _join_segments kernel, and the points are interpolated along their edges all at once. The triangulation is unchanged. The lines are the same segments as tricontour's, to within 1e-9, on eight datasets including node-centered data and faces with repeated nodes; tricontour also left a few slivers ~1e-14 long on the antimeridian, which are gone. A closed loop may start at a different point, and a loop that crosses the antimeridian is then cut there instead, which draws the same. Contouring a noisy field on the 786k-face mesh takes 0.68 s instead of 3.4 s, of which 0.3 s is the triangulation. The splitting of lines at the antimeridian, shared with the edges method, moves to _split_lines_at_antimeridian. Co-Authored-By: Claude Opus 5.5 --- test/test_plot.py | 34 +++++++++++++ uxarray/plot/contour.py | 108 ++++++++++++++++++++++++++++------------ 2 files changed, 110 insertions(+), 32 deletions(-) diff --git a/test/test_plot.py b/test/test_plot.py index 963396e6c..d67f764b3 100644 --- a/test/test_plot.py +++ b/test/test_plot.py @@ -801,6 +801,40 @@ def test_contour_data_not_on_grid(gridpath): plt.close("all") +def test_contour_interpolated_matches_matplotlib(gridpath, datasetpath): + """The interpolated lines are the segments Matplotlib's tricontour gives on the same triangles.""" + from matplotlib.tri import Triangulation + from uxarray.plot.contour import _interpolated_contours, _split_at_antimeridian, _spatial_values, _triangulation + + geoflow = ux.open_dataset(gridpath("ugrid", "geoflow-small", "grid.nc"), datasetpath("ugrid", "geoflow-small", "v1.nc")) + cases = [ + (_face_latitude(ux.open_grid(gridpath("ugrid", "outCSne30", "outCSne30.ug"))), [-30.5, 10.5, 45.5]), + (geoflow["v1"].isel(time=0, meshLayers=0), [-12.0, -6.0, 6.0, 12.0]), + ] + + def segments(lines): + # the ends of each segment of non-zero length, in either order + out = [] + for level, line in lines: + for start, end in zip(line[:-1], line[1:]): + if not np.allclose(start, end, rtol=0, atol=1e-9): + out.append(sorted([(level, *start), (level, *end)])) + return np.array(out).reshape(-1, 6) + + for uxda, levels in cases: + values, dim = _spatial_values(uxda) + lon, lat, triangles, index = _triangulation(uxda.uxgrid, dim) + contour_set = plt.figure().add_subplot().tricontour(Triangulation(lon, lat, triangles), values[index], levels=levels) + expected = segments([(level, part) for level, lines in zip(levels, contour_set.allsegs) + for line in lines if len(line) > 1 for part in _split_at_antimeridian(line)]) + found = segments(_interpolated_contours(values, uxda.uxgrid, dim, np.array(levels))) + + assert len(found) == len(expected) + distance = np.abs(found[:, None, :] - expected[None, :, :]).max(axis=2) + assert distance.min(axis=1).max() < 1e-9 and distance.min(axis=0).max() < 1e-9 + plt.close("all") + + def test_contour_faces_with_repeated_nodes(gridpath): """Faces padded by repeating a node do not give triangles that Matplotlib cannot contour.""" from uxarray.plot.contour import _triangulation diff --git a/uxarray/plot/contour.py b/uxarray/plot/contour.py index ee57f6020..4be78a430 100644 --- a/uxarray/plot/contour.py +++ b/uxarray/plot/contour.py @@ -154,7 +154,8 @@ def _fan_triangles( ) -> np.ndarray: """Splits each face into triangles that share one of its corners. - Used by ``_triangulation``: Matplotlib contours data on triangles only.""" + Used by ``_triangulation``: the method ``"interpolated"`` contours data on + triangles only.""" # A fan from one corner covers a face exactly only if that corner can see every other # corner. Every convex face is like that, but dented faces occur, mostly among the # cells of the dual mesh of irregular meshes (FESOM-style ocean meshes, overlap meshes). @@ -303,36 +304,86 @@ def _interpolated_contours( ) -> list[tuple[float, np.ndarray]]: """Contour lines from linear interpolation on a triangulation of the grid. - The method ``"interpolated"``. The data values sit at the corners of the triangles - from ``_triangulation``, and Matplotlib's ``tricontour`` draws each line straight - across every triangle it passes, where the values change linearly. The lines are - smooth but do not follow the edges of the grid.""" - _raise_hint_if_optional_deps_missing("matplotlib") - from matplotlib.figure import Figure - from matplotlib.tri import Triangulation + The method ``"interpolated"``. The data values sit at the corners of the + triangles from ``_triangulation``, where they change linearly. + Inside each triangle that a level crosses, the line runs between the points + where the level crosses two of its edges. Each point is known by its edge, + which the triangles on either side share, so the pieces are joined into + lines by their edges. The lines are smooth but do not follow the edges of + the grid. + """ lon, lat, triangles, index = _triangulation(uxgrid, dim) values = values[index] # Triangles with a missing value at a corner cannot be contoured finite = np.isfinite(values) - triangles = triangles[finite[triangles].all(axis=1)] + triangles = np.ascontiguousarray( + triangles[finite[triangles].all(axis=1)], dtype=INT_DTYPE + ) if len(triangles) == 0 or len(levels) == 0: return [] - # The contours are computed on an axes that is never displayed - ax = Figure().add_subplot() - contour_set = ax.tricontour( - Triangulation(lon, lat, triangles), np.where(finite, values, 0.0), levels=levels - ) + n_corner = len(values) + contours = [] + for level in levels: + # Each crossed edge as a single number, from its two corners + edges = _crossed_edges(triangles, values, float(level), n_corner) + if len(edges) == 0: + continue + crossed, edge_at = np.unique(edges, return_inverse=True) + edge_nodes, offsets = _join_segments(edge_at[0::2], edge_at[1::2]) + + # The point where the level crosses each edge + first, second = np.divmod(crossed, n_corner) + fraction = (level - values[first]) / (values[second] - values[first]) + x = lon[first] + fraction * (lon[second] - lon[first]) + y = lat[first] + fraction * (lat[second] - lat[first]) + points = np.column_stack([x[edge_nodes], y[edge_nodes]]) + contours.extend(_split_lines_at_antimeridian(float(level), points, offsets)) + return contours - return [ - (float(level), part) - for level, lines in zip(contour_set.levels, contour_set.allsegs) - for line in lines - if len(line) > 1 - for part in _split_at_antimeridian(line) - ] + +@njit(cache=True, nogil=True) +def _crossed_edges(triangles, values, level, n_corner): + """The edges where a level crosses each triangle that it crosses, two per + triangle, each given by its corners ``first < second`` as + ``first * n_corner + second``. A corner is above the level if its value is + greater than it.""" + n_crossed = 0 + for t in range(triangles.shape[0]): + a, b, c = triangles[t, 0], triangles[t, 1], triangles[t, 2] + above = values[a] > level + if (values[b] > level) != above or (values[c] > level) != above: + n_crossed += 1 + + edges = np.empty(2 * n_crossed, dtype=np.int64) + k = 0 + for t in range(triangles.shape[0]): + a, b, c = triangles[t, 0], triangles[t, 1], triangles[t, 2] + for first, second in ((a, b), (b, c), (c, a)): + if (values[first] > level) != (values[second] > level): + if first > second: + first, second = second, first + edges[k] = first * n_corner + second + k += 1 + return edges + + +def _split_lines_at_antimeridian( + level: float, points: np.ndarray, offsets: np.ndarray +) -> list[tuple[float, np.ndarray]]: + """Splits lines stored one after the other, which start at ``offsets[:-1]``, + where they run past -180 or 180 degrees of longitude.""" + starts, ends = offsets[:-1], offsets[1:] + crosses = (np.minimum.reduceat(points[:, 0], starts) < -180.0) | ( + np.maximum.reduceat(points[:, 0], starts) > 180.0 + ) + contours = [] + for line, crossing in zip(np.split(points, ends[:-1]), crosses): + parts = _split_at_antimeridian(line) if crossing else [line] + contours.extend((level, part) for part in parts) + return contours def _join_segments(start: np.ndarray, end: np.ndarray) -> tuple[np.ndarray, np.ndarray]: @@ -461,19 +512,12 @@ def _edge_contours( if len(line_nodes) == 0: continue - # The lines are processed together, as one array of points - starts, ends = offsets[:-1], offsets[1:] + # The lines are processed together, as one array of points. Unwrapped, + # a line that crosses the antimeridian runs past 180 degrees points = np.column_stack( - [_unwrap_lines(lon[line_nodes], starts), lat[line_nodes]] - ) - - # Unwrapped, a line that crosses the antimeridian runs past 180 degrees - crosses = (np.minimum.reduceat(points[:, 0], starts) < -180.0) | ( - np.maximum.reduceat(points[:, 0], starts) > 180.0 + [_unwrap_lines(lon[line_nodes], offsets[:-1]), lat[line_nodes]] ) - for line, crossing in zip(np.split(points, ends[:-1]), crosses): - parts = _split_at_antimeridian(line) if crossing else [line] - contours.extend((float(level), part) for part in parts) + contours.extend(_split_lines_at_antimeridian(float(level), points, offsets)) return contours From f4c2d1f038f7c468ba370d799c3dce25682d9939 Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Mon, 19 Oct 2026 12:20:55 -0500 Subject: [PATCH 5/5] Keep projected contour lines inside the frame of the shaded plot contour-maps frames the lines by the area of the grid, set as their extents, so that they never widen a shaded plot they are overlaid on. Since the lines are projected before GeoViews sees them, HoloViews framed the lines and their labels by the projected line data instead, which reach closer to the antimeridian than the shaded faces: a Robinson map of a field whose lines reach the edge of the map was 34 km wider on each side than the shaded plot. The projected area of the grid is now also the default xlim and ylim of the lines, which HoloViews keeps to; options passed to .plot.contour() still win. The frame then equals the shaded plot's own on Robinson, Robinson centred on 180, orthographic, Lambert conformal and PlateCarree maps, alone and overlaid, with the frame sizes and label positions of contour-maps. Co-Authored-By: Claude Opus 5.5 --- test/test_plot.py | 20 ++++++++++++++++++++ uxarray/plot/accessor.py | 10 +++++----- 2 files changed, 25 insertions(+), 5 deletions(-) diff --git a/test/test_plot.py b/test/test_plot.py index d67f764b3..c6f8e5205 100644 --- a/test/test_plot.py +++ b/test/test_plot.py @@ -735,6 +735,26 @@ def ranges(element): assert ranges(shading * uxda.plot.contour(levels=[-30.5, 30.5], backend="bokeh")) == ranges(shading) +def test_contour_frame_with_lines_at_the_edge_of_the_map(gridpath): + """Projected lines that reach closer to the edge of the map than the shaded faces do not + widen the frame of the shaded plot, alone or overlaid.""" + uxgrid = ux.open_grid(gridpath("mpas", "QU", "oQU480.231010.nc")) + lon, lat = np.radians(uxgrid.face_lon.values), np.radians(uxgrid.face_lat.values) + uxda = ux.UxDataArray(xr.DataArray(np.sin(2 * lat) * np.cos(3 * lon), dims="n_face", name="wave"), uxgrid=uxgrid) + projection = ccrs.Robinson() + + def ranges(element): + figure = hv.renderer("bokeh").get_plot(element).state + return figure.x_range.start, figure.x_range.end, figure.y_range.start, figure.y_range.end + + shading = uxda.plot(projection=projection, backend="bokeh") + contours = uxda.plot.contour(levels=[-0.6, -0.2, 0.2, 0.6], projection=projection, backend="bokeh") + # the lines do reach past the shaded faces + assert np.nanmax(np.abs(contours.get(0).dimension_values("x"))) > max(np.abs(ranges(shading)[:2])) + np.testing.assert_allclose(ranges(contours), ranges(shading), rtol=1e-12) + np.testing.assert_allclose(ranges(shading * contours), ranges(shading), rtol=1e-12) + + def test_contour_missing_values(gridpath): """Faces without data are left out instead of raising.""" uxgrid = ux.open_grid(gridpath("ugrid", "outCSne30", "outCSne30.ug")) diff --git a/uxarray/plot/accessor.py b/uxarray/plot/accessor.py index 1a1be86e2..201d50f26 100644 --- a/uxarray/plot/accessor.py +++ b/uxarray/plot/accessor.py @@ -687,12 +687,9 @@ def contour( # The lines are already projected, so the area of the grid is projected # as GeoViews projects the extents of an element in longitude and latitude + extents = project_extents(extents, ccrs.PlateCarree(), projection) lines = _projected_contours()( - paths, - kdims=kdims, - vdims=[level_dim], - crs=projection, - extents=project_extents(extents, ccrs.PlateCarree(), projection), + paths, kdims=kdims, vdims=[level_dim], crs=projection, extents=extents ) text = gv.Labels( label_data, kdims=kdims, vdims=["text"], crs=ccrs.PlateCarree() @@ -709,6 +706,9 @@ def contour( defaults = {"padding": 0} if projection is not None: defaults["data_aspect"] = 1 + # HoloViews frames the projected lines by their data, which reach closer + # to the antimeridian than the shaded faces, unless the area is a limit + defaults.update(xlim=extents[0::2], ylim=extents[1::2]) else: defaults.update(width=700, height=300) kwargs = {**defaults, **kwargs}