From 5134f285e07a5c8a53f51e53e52839f60468ec30 Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Mon, 21 Sep 2026 18:19:15 -0500 Subject: [PATCH 1/4] Wrap longitudes elementwise so Grid construction stops computing _set_desired_longitude_range decided whether to wrap by asking lon.max() > 180. On a dask-backed coordinate that reduction is a compute, and Grid.__init__ calls it -- so opening a chunked grid read and reduced every longitude array before the caller had asked for anything, and did it again on every isel and every copy(), each of which builds a new Grid. The wrap is now xr.where((lon > 180) | (lon < -180), (lon + 180) % 360 - 180, lon): elementwise, lazy, chunk-parallel, no reduction. Measured on a synthetic 4M-node UGRID file, chunked at 500k nodes, best of 5: dask computes wall tracemalloc peak open_grid (before) 3 23.3ms 12.2 MB open_grid (after) 2 16.2ms 1.0 MB isel (before) 2 isel (after) 1 The computes that remain are a separate site on connectivity rather than coordinates: _standardize_connectivity's conn.isnull().any() in io/_ugrid.py, reached twice on the read path, and _slice_face_indices in grid/slice.py materializing the connectivity it slices by. Neither is touched here. Three behavioural differences, all from doing this per element rather than per array. * In-range longitudes are now left exactly alone. (lon + 180) - 180 does not round-trip, so wrapping the whole array perturbed values that were already in range by up to 3e-14 degrees -- in outCSne30, 2.1182935e-14 became 2.8421709e-14. On the elements that do need wrapping the two forms are bit-identical. * Longitudes below -180 are normalized. The old test was on the maximum alone, so it reached the negative tail only when the same array also held a value above 180. * Both endpoints are kept, so the interval is the closed [-180, 180]. Folding 180.0 to -180.0 would match _xyz_to_lonlat_deg, which wraps unconditionally into the half-open interval, but it breaks antimeridian_face_indices: that reads a face as crossing from the span of its longitudes, and a face with one vertex at 180 and the rest near -170 goes from a span of 350 to a span of 10 and disappears. Caught by test_antimeridian_point_on and test_to_geodataframe_preserves_antimeridian_faces. Each variable is wrapped at most once, keyed on the xr.Variable object. Without that, edge_lat -- which calls this on every property access, outside its populate guard -- would stack a where layer onto the graph per access. Keying on the Variable makes the memo self-invalidating: assigning into _ds replaces that object, so a repopulated or user-assigned coordinate is wrapped again. The two Exodus round-trip tests compared a grid against its own lon -> xyz -> lon reload with assert_allclose(rtol=1e-8), and passed only because the old whole-array wrap applied to the original the identical perturbation the reload applies. With the original left alone, the reload's own error is exposed, and rtol is the wrong instrument for it twice over: a longitude near zero has no magnitude for a relative tolerance to measure against (outCSne30 nodes 4372, 4749, 7e-15 degrees apart), and longitude is periodic, so a node on or one ulp short of the antimeridian reads 180.0 on the original and -180.0 on the reload -- the same meridian, scored as a 360-degree error (179 nodes of outRLL1deg; outCSne30 nodes 3966, 5155). Those two assertions now compare the difference modulo 360 to an absolute tolerance, still ERROR_TOLERANCE. The helper still catches a 1e-7 shift and still rejects an antipode. Tier 0.2 of the chunked refactor plan. Test suite: 962 passed, 1 skipped. test_plot_with_features fails identically before and after (matplotlib figure size, unrelated). Co-Authored-By: Claude Opus 5 --- benchmarks/lazy_grid_construction.py | 72 ++++++++++++ test/grid/grid/test_io.py | 33 +++++- test/grid/grid/test_lon_range.py | 169 +++++++++++++++++++++++++++ uxarray/grid/coordinates.py | 50 +++++++- uxarray/grid/grid.py | 3 + 5 files changed, 319 insertions(+), 8 deletions(-) create mode 100644 benchmarks/lazy_grid_construction.py create mode 100644 test/grid/grid/test_lon_range.py diff --git a/benchmarks/lazy_grid_construction.py b/benchmarks/lazy_grid_construction.py new file mode 100644 index 000000000..f8f6f4b80 --- /dev/null +++ b/benchmarks/lazy_grid_construction.py @@ -0,0 +1,72 @@ +"""What ``Grid`` construction pulls off disk before the caller asks for it. + +A chunked ``open_grid`` is supposed to cost metadata and nothing else. It did +not: ``_set_desired_longitude_range`` decided whether to wrap by asking +``lon.max() > 180``, and on a dask-backed coordinate that reduction is a +compute -- one whole longitude array read and reduced per constructor call, +and the constructor runs again on every ``isel`` and every ``copy()``. + +Wall time will not show this on a test-sized mesh, and it is the wrong +instrument anyway: the quantity that changed is discrete. So the number here +is a count of dask graph executions, which is exact, has no variance, and +moves by one the moment the reduction comes back. +""" + +import os +from pathlib import Path + +from dask.callbacks import Callback + +import uxarray as ux + +current_path = Path(os.path.dirname(os.path.realpath(__file__))).parents[0] + +grid_path = current_path / "test" / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug" + +#: Small enough to leave several chunks per coordinate, which is what makes a +#: reduction over one visible as a scheduled graph rather than a fused no-op. +CHUNKS = {"n_node": 1000} + + +class _CountComputes(Callback): + def __init__(self): + self.n = 0 + + def _start(self, dsk): + self.n += 1 + + +class LazyGridConstruction: + def setup(self): + # Opening a grid for the first time in a process pulls in xarray's + # backend machinery and the netCDF library + ux.open_grid(grid_path, chunks=CHUNKS) + + def track_computes_open_grid_chunked(self): + """Dask graph executions during a chunked ``open_grid``. + + Three before the longitude wrap went elementwise, two after. The two + that remain are ``_standardize_connectivity``'s ``conn.isnull().any()`` + in ``io/_ugrid.py``, reached once from ``match_chunks_to_ugrid`` and + once from ``Grid.from_dataset`` -- a separate site, on connectivity + rather than coordinates, and not addressed here. + """ + counter = _CountComputes() + with counter: + ux.open_grid(grid_path, chunks=CHUNKS) + return counter.n + + def track_computes_isel(self): + """Dask graph executions during a subset of an already-open grid. + + ``isel`` builds a new ``Grid``, so it paid the constructor's reduction + on every call. Two before, one after. What remains is + ``_slice_face_indices`` (``grid/slice.py``) materializing the + face-node connectivity it slices by -- again connectivity, not + coordinates. + """ + uxgrid = ux.open_grid(grid_path, chunks=CHUNKS) + counter = _CountComputes() + with counter: + uxgrid.isel(n_face=slice(0, 100)) + return counter.n diff --git a/test/grid/grid/test_io.py b/test/grid/grid/test_io.py index db6ce6a1d..1c8b084c3 100644 --- a/test/grid/grid/test_io.py +++ b/test/grid/grid/test_io.py @@ -6,6 +6,33 @@ from uxarray.constants import ERROR_TOLERANCE +def _assert_lon_close(actual, desired, err_msg, atol=ERROR_TOLERANCE): + """Compare longitudes as directions on a circle, to an absolute tolerance. + + Used for the Exodus round-trips, which are the ones that go + lon/lat -> xyz -> lon/lat and so come back through ``_xyz_to_lonlat_deg``. + ``assert_allclose(..., rtol=...)`` is the wrong instrument for that on two + counts. A longitude near zero has no magnitude for a relative tolerance to + be measured against -- 2.1e-14 vs 2.8e-14 is 7e-15 degrees apart and fails + at rtol=1e-8 (``outCSne30`` nodes 4372 and 4749). And longitude is + periodic: ``_xyz_to_lonlat_deg`` wraps into the half-open [-180, 180), + while ``_set_desired_longitude_range`` keeps both endpoints, so a node on + the antimeridian reads 180.0 on the original and -180.0 on the reload -- + the same meridian, scored as a 360-degree error (179 nodes of + ``outRLL1deg``, and ``outCSne30`` nodes 3966 and 5155, which sit one ulp + short of 180 and land on -180.0 once rounded through Cartesian). + + All of it was masked before. The old ``_set_desired_longitude_range`` + wrapped the whole array whenever any element exceeded 180, applying to the + original grid the identical ``(lon + 180) % 360 - 180`` that the reload + applies -- so both sides carried the same perturbation and the same + endpoint fold. The wrap is elementwise now and leaves in-range longitudes + alone, which leaves the round-trip's own error exposed. + """ + diff = (np.asarray(actual) - np.asarray(desired) + 180.0) % 360.0 - 180.0 + np.testing.assert_allclose(diff, 0.0, atol=atol, err_msg=err_msg) + + def test_normalize_existing_coordinates_non_norm_initial(gridpath): from uxarray.grid.validation import _check_normalization uxgrid = ux.open_grid(gridpath("mpas", "QU", "mesh.QU.1920km.151026.nc")) @@ -145,11 +172,10 @@ def test_grid_ugrid_exodus_roundtrip(gridpath): err_msg=f"UGRID longitude mismatch for {grid_name}", rtol=ERROR_TOLERANCE ) - np.testing.assert_allclose( + _assert_lon_close( original_grid.node_lon.values, reloaded_exodus.node_lon.values, err_msg=f"Exodus longitude mismatch for {grid_name}", - rtol=ERROR_TOLERANCE ) np.testing.assert_allclose( original_grid.node_lat.values, @@ -185,11 +211,10 @@ def test_exodus_roundtrip_rll1deg_node_lonlat(gridpath, tmp_path): reloaded_exodus = ux.open_grid(exodus_filepath) - np.testing.assert_allclose( + _assert_lon_close( grid.node_lon.values, reloaded_exodus.node_lon.values, err_msg="Exodus longitude mismatch for RLL1deg", - rtol=ERROR_TOLERANCE, ) np.testing.assert_allclose( grid.node_lat.values, diff --git a/test/grid/grid/test_lon_range.py b/test/grid/grid/test_lon_range.py new file mode 100644 index 000000000..93d1c3c94 --- /dev/null +++ b/test/grid/grid/test_lon_range.py @@ -0,0 +1,169 @@ +"""Guards on the [-180, 180] longitude wrap. + +``_set_desired_longitude_range`` used to decide whether to wrap by asking +``lon.max() > 180``. On a dask-backed grid that reduction is a compute, and +``Grid.__init__`` calls it, so opening a grid pulled its longitude +coordinates into memory. The wrap is now elementwise. + +Two things have to hold for that to be an improvement rather than a trade: +the construction must stay lazy, and repeated calls must not pile ``where`` +layers onto the graph -- ``edge_lat`` invokes it on every property access. +""" + +import numpy as np +import numpy.testing as nt +import pytest +import xarray as xr +from dask.callbacks import Callback + +import uxarray as ux +from uxarray.grid.coordinates import _set_desired_longitude_range + + +class _CountComputes(Callback): + """Counts dask graph executions inside the ``with`` block.""" + + def __init__(self): + self.n = 0 + + def _start(self, dsk): + self.n += 1 + + +class _Shim: + """The only attribute ``_set_desired_longitude_range`` needs of a Grid.""" + + def __init__(self, ds): + self._ds = ds + + +def _shim(values, name="node_lon", **attrs): + return _Shim(xr.Dataset({name: (f"n_{name.split('_')[0]}", values, attrs)})) + + +def test_wrap_is_elementwise_and_leaves_in_range_values_exact(): + """Outside [-180, 180] wraps; inside is passed through bit-for-bit. + + The old whole-array form perturbed in-range longitudes by up to ~3e-14 + degrees, because ``(lon + 180) - 180`` does not round-trip exactly. The + assertion here is ``assert_array_equal``, not ``allclose``, on purpose. + """ + rng = np.random.default_rng(0) + values = np.concatenate( + [rng.uniform(0.0, 180.0, 500), rng.uniform(180.0, 360.0, 500)] + ) + + grid = _shim(values.copy()) + _set_desired_longitude_range(grid) + wrapped = grid._ds["node_lon"].values + + in_range = values <= 180.0 + nt.assert_array_equal(wrapped[in_range], values[in_range]) + nt.assert_array_equal( + wrapped[~in_range], (values[~in_range] + 180.0) % 360.0 - 180.0 + ) + assert wrapped.max() <= 180.0 and wrapped.min() >= -180.0 + + +def test_both_endpoints_of_the_antimeridian_are_left_alone(): + """180.0 stays 180.0 and -180.0 stays -180.0, however the array looks. + + Not a detail. ``antimeridian_face_indices`` reads a face as crossing from + the span of its longitudes, so a face with one vertex at exactly 180 and + the rest near -170 is only visible while that vertex reads as 180 -- fold + it to -180 and the span drops from 350 to 10. Under the old reduction the + endpoint survived only by accident: an array whose maximum was exactly + 180 was not wrapped at all, so nothing touched it. + """ + for companion in (10.0, 270.0, -170.0): + grid = _shim(np.array([180.0, companion])) + _set_desired_longitude_range(grid) + assert grid._ds["node_lon"].values[0] == 180.0 + + grid = _shim(np.array([-180.0, 0.0, 179.5])) + _set_desired_longitude_range(grid) + nt.assert_array_equal(grid._ds["node_lon"].values, [-180.0, 0.0, 179.5]) + + +def test_wrap_normalizes_both_tails_and_is_idempotent(): + """One pass lands everything in [-180, 180); a second pass changes nothing. + + Both tails matter. The reduction tested ``lon.max() > 180``, so it reached + longitudes below -180 only when the same array happened to hold one above + 180, and left them alone otherwise. + """ + rng = np.random.default_rng(1) + grid = _shim(rng.uniform(-720.0, 720.0, 1000)) + + _set_desired_longitude_range(grid) + once = grid._ds["node_lon"].values.copy() + assert once.min() >= -180.0 and once.max() <= 180.0 + + grid._wrapped_lon_vars = {} + _set_desired_longitude_range(grid) + nt.assert_array_equal(grid._ds["node_lon"].values, once) + + +@pytest.mark.parametrize("dtype", [np.float64, np.float32]) +def test_wrap_preserves_dtype_name_and_attrs(dtype): + grid = _shim(np.array([0.0, 180.0, 270.0], dtype=dtype), units="degrees_east") + _set_desired_longitude_range(grid) + out = grid._ds["node_lon"] + + assert out.dtype == dtype + assert out.name == "node_lon" + assert out.attrs == {"units": "degrees_east"} + + +def test_nan_longitudes_survive_the_wrap(): + """``nan > 180`` is False, so NaN takes the pass-through branch.""" + grid = _shim(np.array([np.nan, 270.0, 10.0])) + _set_desired_longitude_range(grid) + + nt.assert_array_equal(grid._ds["node_lon"].values, [np.nan, -90.0, 10.0]) + + +def test_open_grid_does_not_compute_longitudes(gridpath): + """The regression this change exists to fix. + + Counts every dask execution during ``open_grid``, then asserts the wrap + contributed none of them by re-running the wrap on the constructed grid + under its own counter. + """ + path = gridpath("ugrid", "outCSne30", "outCSne30.ug") + + grid = ux.open_grid(path, chunks={"n_node": 1000}) + assert grid._ds["node_lon"].chunks is not None, "fixture is not chunked" + + counter = _CountComputes() + with counter: + _set_desired_longitude_range(grid) + assert counter.n == 0 + + +def test_repeated_calls_do_not_grow_the_graph(gridpath): + """``edge_lat`` calls the wrap on every access, outside its populate guard. + + Without the memo, an unconditional elementwise wrap would add a ``where`` + layer per call -- lazy, so nothing would fail, just an ever-deepening + graph. + """ + path = gridpath("ugrid", "outCSne30", "outCSne30.ug") + grid = ux.open_grid(path, chunks={"n_node": 1000}) + + before = len(grid._ds["node_lon"].data.dask) + for _ in range(20): + _set_desired_longitude_range(grid) + assert len(grid._ds["node_lon"].data.dask) == before + + +def test_memo_reopens_when_the_variable_is_replaced(gridpath): + """The memo keys on the ``xr.Variable``, so assignment invalidates it.""" + path = gridpath("ugrid", "outCSne30", "outCSne30.ug") + grid = ux.open_grid(path) + + grid._ds["node_lon"] = grid._ds["node_lon"] + 200.0 + assert float(grid._ds["node_lon"].max()) > 180.0 + + _set_desired_longitude_range(grid) + assert float(grid._ds["node_lon"].max()) <= 180.0 diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 73ad3aed6..a6e0ac576 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -727,6 +727,41 @@ def _is_projected_grid(uxgrid) -> bool: def _set_desired_longitude_range(uxgrid): """Sets the longitude range to [-180, 180] for all longitude variables. + The wrap is elementwise rather than guarded by ``lon.max() > 180``. The + guard was a reduction, and a reduction on a dask-backed coordinate is a + compute -- so ``Grid.__init__``, which calls this, pulled every longitude + array into memory before the caller had asked for anything. The + elementwise form stays in the graph and splits over chunks. + + On the elements that need wrapping the two are bit-identical -- same + expression, same order. They differ on the elements that do not: the + reduction ran ``(lon + 180) % 360 - 180`` over the whole array once any + value exceeded 180, and that round-trip is not exact, so it perturbed + in-range longitudes by up to ~3e-14 degrees. ``xr.where`` passes those + through untouched. + + An element is wrapped when it falls outside [-180, 180] on *either* side. + The negative half of that predicate is what the reduction was missing: it + tested the maximum alone, so a longitude below -180 was normalized only + when the same array happened to also hold one above 180, and was left + where it was otherwise. + + Both endpoints are kept, rather than folding 180.0 onto -180.0 for a + half-open [-180, 180). ``_xyz_to_lonlat_deg`` does produce the half-open + interval, so a grid reloaded through Cartesian coordinates disagrees with + its original at a node sitting exactly on the antimeridian -- but that + node is the *point* of ``antimeridian_face_indices``, which reads a face + as crossing from the span of its longitudes. Folding 180.0 to -180.0 + collapses the span of a face that touches the antimeridian from the west + and hides it. The round-trip comparison is the cheaper of the two to make + periodic, and the Exodus tests do that. + + Each variable is wrapped at most once. ``edge_lat`` calls this on every + access, outside its populate guard, so without the memo an unconditional + wrap would stack a ``where`` layer onto the graph per property access. Keying on the ``xr.Variable`` object makes the memo + self-invalidating: assigning into ``_ds`` replaces that object, so a + repopulated or user-assigned coordinate is wrapped again. + Skipped entirely for projected grids: wrapping meter-scale coordinates as if they were degrees silently corrupts the geometry. A ``UserWarning`` is issued once per Grid instance so users know which operations are invalid. @@ -747,16 +782,23 @@ def _set_desired_longitude_range(uxgrid): uxgrid._projected_warning_issued = True return + memo = getattr(uxgrid, "_wrapped_lon_vars", None) + if memo is None: + memo = uxgrid._wrapped_lon_vars = {} + with xr.set_options(keep_attrs=True): for lon_name in ["node_lon", "edge_lon", "face_lon"]: if lon_name in uxgrid._ds: da = uxgrid._ds[lon_name] if da.size == 0: continue - if da.max() > 180: - wrapped = (uxgrid._ds[lon_name] + 180) % 360 - 180 - wrapped.name = da.name - uxgrid._ds[lon_name] = wrapped + if memo.get(lon_name) is da.variable: + continue + out_of_range = (da > 180) | (da < -180) + wrapped = xr.where(out_of_range, (da + 180) % 360 - 180, da) + wrapped.name = da.name + uxgrid._ds[lon_name] = wrapped + memo[lon_name] = uxgrid._ds[lon_name].variable def prepare_points(points, normalize): diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 2d0d3321d..c74b4f4f0 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -245,6 +245,9 @@ def __init__( # flag to ensure projected-grid warning fires only once per instance self._projected_warning_issued = False + # longitude Variables already wrapped into [-180, 180], by name + self._wrapped_lon_vars = {} + # set desired longitude range to [-180, 180] _set_desired_longitude_range(self) From 76c90548bc65f9cb5964841507378fef1b8af062 Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Wed, 23 Sep 2026 13:52:41 -0500 Subject: [PATCH 2/4] Keep the eager path off the elementwise longitude wrap The elementwise wrap is the right shape for a dask-backed coordinate and the wrong one for an array already in memory. There was never a compute to defer on that path, only a scan, and elementwise costs a pass plus ~18 bytes per node of temporaries -- two bool masks, the arithmetic, the result -- where a reduction costs a pass and allocates nothing. That matters because every open_grid/open_dataset call in the benchmark suite is eager; none pass chunks=. So the previous commit, measured there, was a regression and nothing else. Measured on a synthetic 4M-node UGRID file, best of 5. The wrap in isolation: in range 0..360 old reduction 4.6ms 0MB 36.8ms 64MB elementwise only 20.9ms 72MB 38.1ms 72MB elementwise + guard 9.3ms 0MB 42.7ms 72MB and through eager open_grid, where the file read dominates and the peak does not move at all (180.0 MB in every arm): in range 0..360 base (cmd/nogil) 145.5ms 181.9ms elementwise only 164.4ms 186.1ms elementwise + guard 147.9ms 185.9ms _lon_within_range is a guard, not a decision: when it is true the wrap is the identity on every element, so skipping it cannot change a value. test_eager_fast_path_agrees_with_the_where_element_for_element asserts that directly against the unguarded expression rather than assuming it, over five inputs including both endpoints, the negative tail and NaN. The reductions are the thing that made the old code compute, so they are allowed only where there is nothing to defer -- da.chunks is None. A dask-backed array skips the guard entirely, which test_guard_is_skipped_for_dask_backed_arrays pins by asserting the guard does compute when handed one. The chunked numbers are unchanged: 2 graph executions, 16.0ms, 1.0 MB peak. `and` short-circuits, so an array that does need wrapping usually pays a single max -- the same reduction the old code paid -- before falling through. The 0..360 column above is that extra max: ~4.6ms on a 190ms open_grid. Test suite: 973 passed, 1 skipped. test_plot_with_features fails identically before and after (matplotlib figure size, unrelated). Co-Authored-By: Claude Opus 5 --- test/grid/grid/test_lon_range.py | 51 +++++++++++++++++++++++++++++++- uxarray/grid/coordinates.py | 31 +++++++++++++++++++ 2 files changed, 81 insertions(+), 1 deletion(-) diff --git a/test/grid/grid/test_lon_range.py b/test/grid/grid/test_lon_range.py index 93d1c3c94..258815a76 100644 --- a/test/grid/grid/test_lon_range.py +++ b/test/grid/grid/test_lon_range.py @@ -17,7 +17,10 @@ from dask.callbacks import Callback import uxarray as ux -from uxarray.grid.coordinates import _set_desired_longitude_range +from uxarray.grid.coordinates import ( + _lon_within_range, + _set_desired_longitude_range, +) class _CountComputes(Callback): @@ -167,3 +170,49 @@ def test_memo_reopens_when_the_variable_is_replaced(gridpath): _set_desired_longitude_range(grid) assert float(grid._ds["node_lon"].max()) <= 180.0 + + +def test_eager_fast_path_agrees_with_the_where_element_for_element(): + """The in-memory guard must be a pure optimization. + + ``_lon_within_range`` lets the eager path skip the wrap entirely. That is + only sound if the wrap would have been the identity, so the two are + compared directly rather than the skip being assumed correct. + """ + rng = np.random.default_rng(2) + cases = { + "strictly inside": rng.uniform(-179.0, 179.0, 500), + "on both endpoints": np.array([-180.0, 180.0, 0.0, 179.5, -179.5]), + "needs wrapping": rng.uniform(0.0, 360.0, 500), + "negative tail": rng.uniform(-540.0, -180.5, 500), + "with nan": np.array([np.nan, 10.0, 200.0]), + } + + for label, values in cases.items(): + da = xr.DataArray(values, dims="n_node", name="node_lon") + unguarded = xr.where( + (da > 180) | (da < -180), (da + 180) % 360 - 180, da + ).values + + grid = _shim(values.copy()) + _set_desired_longitude_range(grid) + + nt.assert_array_equal(grid._ds["node_lon"].values, unguarded, err_msg=label) + + +def test_guard_is_skipped_for_dask_backed_arrays(gridpath): + """The guard is two reductions, which is exactly what must not run lazily. + + ``test_open_grid_does_not_compute_longitudes`` would catch this too, but + only as a compute count; this names the reason. + """ + grid = ux.open_grid( + gridpath("ugrid", "outCSne30", "outCSne30.ug"), chunks={"n_node": 1000} + ) + da = grid._ds["node_lon"] + assert da.chunks is not None + + counter = _CountComputes() + with counter: + _lon_within_range(da) + assert counter.n > 0, "guard is lazy here; skipping it would be pointless" diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index a6e0ac576..ed4bc92c6 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -724,6 +724,29 @@ def _is_projected_grid(uxgrid) -> bool: return False +def _lon_within_range(da) -> bool: + """Whether every element of ``da`` already lies in [-180, 180]. + + A guard, not a decision. When this is true the wrap below is the identity + on every element, so skipping it cannot change a value. + + It exists to keep the eager path off the ``where``. The wrap has to be + elementwise to stay out of a dask compute, but elementwise costs a pass + plus ~18 bytes per node of temporaries -- two bool masks, the arithmetic, + the result -- where a reduction costs a pass and allocates nothing. On an + array that is already in range, measured at 4M nodes, that is 20.9ms and + 72 MB to hand back the input unchanged, against 4.6ms and nothing. + Reductions are what made the old code compute, so this one is allowed + only where there is nothing to defer: ``da.chunks is None``. A + dask-backed array skips it and goes straight to the ``where``. + + ``and`` short-circuits, so an array that does need wrapping usually pays + a single ``max`` -- the same reduction the old code paid -- before + falling through. + """ + return bool(da.max() <= 180 and da.min() >= -180) + + def _set_desired_longitude_range(uxgrid): """Sets the longitude range to [-180, 180] for all longitude variables. @@ -756,6 +779,11 @@ def _set_desired_longitude_range(uxgrid): and hides it. The round-trip comparison is the cheaper of the two to make periodic, and the Exodus tests do that. + None of that applies to a grid whose coordinates are already in memory -- + there was never a compute to defer there, only a scan -- and elementwise + is the slower shape for it. ``_lon_within_range`` keeps the eager path off + the ``where`` when the wrap would be the identity anyway. + Each variable is wrapped at most once. ``edge_lat`` calls this on every access, outside its populate guard, so without the memo an unconditional wrap would stack a ``where`` layer onto the graph per property access. Keying on the ``xr.Variable`` object makes the memo @@ -794,6 +822,9 @@ def _set_desired_longitude_range(uxgrid): continue if memo.get(lon_name) is da.variable: continue + if da.chunks is None and _lon_within_range(da): + memo[lon_name] = da.variable + continue out_of_range = (da > 180) | (da < -180) wrapped = xr.where(out_of_range, (da + 180) % 360 - 180, da) wrapped.name = da.name From 2fe2897c7e8dc5d81f9bfd76bde6030a400dc827 Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Wed, 23 Sep 2026 14:05:37 -0500 Subject: [PATCH 3/4] Benchmark a chunked open_grid across both oQU resolutions Every other open_grid in the suite is eager, so none of it can see what Grid construction does to a dask-backed grid -- which is the path the lazy longitude wrap changed, and the path the rest of the chunked refactor will keep changing. OpenGridChunked adds time_open_grid and track_peakmem_open_grid, parametrized over the oQU 480km and 120km meshes already registered in helpers/_fixtures.py. Chunks are held at N_CHUNKS=4 per grid dimension, so the graph is the same shape at both resolutions and only the data under it grows. The file is read directly rather than through CachedFixtures, because reading it is the subject. On MPAS this is where the lazy wrap matters most. node_lon, edge_lon and face_lon all exist at construction, so the old max() > 180 check ran three computes per open; the branch runs none. Measured with this benchmark's own setup, best of 15, against HEAD's tree with coordinates.py taken from cmd/nogil: time tracemalloc peak base branch base branch 480km 60.5ms 56.5ms 3.27MB 3.27MB 120km 42.3ms 38.0ms 3.53MB 4.02MB The 120km peak reads higher on the branch, and it is not data. With gc disabled, building the where graph allocates ~1.9 MB of transient objects, nearly all in inspect.signature via dask/xarray op dispatch -- the same at both resolutions and with a single chunk, so it does not scale with the grid. At 480km it sits under the HDF5 read's own high-water mark; at 120km the netCDF3 read is cheap enough that it becomes the peak. After a gc.collect() the branch retains ~30 kB more than base, which is the extra graph layers. Worth knowing for later steps: this benchmark sees graph-construction overhead, not only bytes read. Three things the numbers above depend on: * Compare each resolution to its own history, not to the other. The two files are different formats -- oQU480.grid.nc is netCDF4/HDF5, oQU120.grid.nc is netCDF3 -- and the HDF5 open costs more, so 480km reads slower than 120km despite a sixteenth of the data. * Each open gets a fresh copy of the chunks dict. match_chunks_to_ugrid (core/utils.py) writes the source-format dimension names into the dict it is handed, so a reused one gives every sample after the first a different argument. * The warning "The specified chunks separate the stored chunks" is filtered. oQU480 stores layerThickness, ssh and zMid as one chunk of all 1791 cells, so any n_face chunking splits them and xarray warns once per open -- about data variables the grid reader drops. Checked by calling the class the way asv does (setup(param), then the time_/track_ methods) under -W error::UserWarning. Not run through asv itself: its discovery subprocess cannot import uxarray from the uxarray conda env, where the package is not installed. Co-Authored-By: Claude Opus 5.5 --- benchmarks/lazy_grid_construction.py | 69 ++++++++++++++++++++++++++-- 1 file changed, 65 insertions(+), 4 deletions(-) diff --git a/benchmarks/lazy_grid_construction.py b/benchmarks/lazy_grid_construction.py index f8f6f4b80..230f94612 100644 --- a/benchmarks/lazy_grid_construction.py +++ b/benchmarks/lazy_grid_construction.py @@ -6,19 +6,25 @@ compute -- one whole longitude array read and reduced per constructor call, and the constructor runs again on every ``isel`` and every ``copy()``. -Wall time will not show this on a test-sized mesh, and it is the wrong -instrument anyway: the quantity that changed is discrete. So the number here -is a count of dask graph executions, which is exact, has no variance, and -moves by one the moment the reduction comes back. +Two instruments. ``LazyGridConstruction`` counts dask graph executions on a +test-sized mesh: exact, no variance, and it moves by one the moment a +reduction comes back. ``OpenGridChunked`` is the wall-time and peak-memory +view of the same thing across two resolutions, which is what the rest of the +suite -- every other ``open_grid`` in it is eager -- cannot see. """ +import math import os +import warnings from pathlib import Path from dask.callbacks import Callback import uxarray as ux +from .helpers._fixtures import OQU_GRIDS, OQU_RESOLUTIONS +from .helpers._peakmem import peak_allocated + current_path = Path(os.path.dirname(os.path.realpath(__file__))).parents[0] grid_path = current_path / "test" / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug" @@ -70,3 +76,58 @@ def track_computes_isel(self): with counter: uxgrid.isel(n_face=slice(0, 100)) return counter.n + + +#: Chunks per grid dimension, held fixed across resolutions so the graph is +#: the same shape at both and only the data under it grows. +N_CHUNKS = 4 + + +class OpenGridChunked: + """``open_grid`` with ``chunks=``, across both oQU resolutions. + + Reads the file directly rather than through ``CachedFixtures``, because + reading the file is the subject here. The 120km mesh is ~16x the 480km one + (59,329 nodes against 3,947). + + Track each resolution against its own history, not against the other. The + two files are different formats -- oQU480.grid.nc is netCDF4/HDF5, + oQU120.grid.nc is netCDF3 -- and the HDF5 open costs more, so 480km reads + *slower* than 120km despite a sixteenth of the data. + """ + + param_names = ["resolution"] + params = [OQU_RESOLUTIONS] + + def setup(self, resolution): + self.grid_path = OQU_GRIDS[resolution] + # An eager open for the sizes, which also pays the first-open cost of + # the netCDF backend here rather than in the first sample. + uxgrid = ux.open_grid(self.grid_path) + self.chunks = { + "n_node": math.ceil(uxgrid.n_node / N_CHUNKS), + "n_face": math.ceil(uxgrid.n_face / N_CHUNKS), + } + self._open() + + def _open(self): + with warnings.catch_warnings(): + # oQU480 stores layerThickness, ssh and zMid as one chunk of all + # 1791 cells, so any n_face chunking splits them and xarray warns + # once per open. They are data variables the grid reader drops. + warnings.filterwarnings( + "ignore", message="The specified chunks separate the stored chunks" + ) + # A copy: open_grid adds the source-format dimension names to the + # dict it is handed (match_chunks_to_ugrid), so reusing one would + # have every sample after the first open with a different argument. + return ux.open_grid(self.grid_path, chunks=dict(self.chunks)) + + def time_open_grid(self, resolution): + self._open() + + def track_peakmem_open_grid(self, resolution): + """Transient high-water allocation of a chunked ``open_grid``.""" + return peak_allocated(self._open) + + track_peakmem_open_grid.unit = "bytes" From c92cd2291430dafb6c904d64c76c059349f0d33d Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Fri, 25 Sep 2026 19:22:07 -0500 Subject: [PATCH 4/4] lazy lon: deslop --- benchmarks/lazy_grid_construction.py | 45 +++++-------------- test/grid/grid/test_io.py | 26 +++-------- test/grid/grid/test_lon_range.py | 65 ++++++--------------------- uxarray/grid/coordinates.py | 67 +++++----------------------- 4 files changed, 40 insertions(+), 163 deletions(-) diff --git a/benchmarks/lazy_grid_construction.py b/benchmarks/lazy_grid_construction.py index 230f94612..7406f9437 100644 --- a/benchmarks/lazy_grid_construction.py +++ b/benchmarks/lazy_grid_construction.py @@ -1,16 +1,8 @@ -"""What ``Grid`` construction pulls off disk before the caller asks for it. - -A chunked ``open_grid`` is supposed to cost metadata and nothing else. It did -not: ``_set_desired_longitude_range`` decided whether to wrap by asking -``lon.max() > 180``, and on a dask-backed coordinate that reduction is a -compute -- one whole longitude array read and reduced per constructor call, -and the constructor runs again on every ``isel`` and every ``copy()``. - -Two instruments. ``LazyGridConstruction`` counts dask graph executions on a -test-sized mesh: exact, no variance, and it moves by one the moment a -reduction comes back. ``OpenGridChunked`` is the wall-time and peak-memory -view of the same thing across two resolutions, which is what the rest of the -suite -- every other ``open_grid`` in it is eager -- cannot see. +"""Cost of constructing a ``Grid`` from a chunked open. + +``LazyGridConstruction`` counts dask graph executions; ``OpenGridChunked`` times +and peak-measures a chunked ``open_grid`` at two resolutions, which the rest of +the suite, all eager, does not cover. """ import math @@ -51,11 +43,8 @@ def setup(self): def track_computes_open_grid_chunked(self): """Dask graph executions during a chunked ``open_grid``. - Three before the longitude wrap went elementwise, two after. The two - that remain are ``_standardize_connectivity``'s ``conn.isnull().any()`` - in ``io/_ugrid.py``, reached once from ``match_chunks_to_ugrid`` and - once from ``Grid.from_dataset`` -- a separate site, on connectivity - rather than coordinates, and not addressed here. + Two remain after the lazy wrap, both from ``_standardize_connectivity`` + checking connectivity for nulls. """ counter = _CountComputes() with counter: @@ -63,13 +52,9 @@ def track_computes_open_grid_chunked(self): return counter.n def track_computes_isel(self): - """Dask graph executions during a subset of an already-open grid. + """Dask graph executions during ``isel`` on a chunked grid. - ``isel`` builds a new ``Grid``, so it paid the constructor's reduction - on every call. Two before, one after. What remains is - ``_slice_face_indices`` (``grid/slice.py``) materializing the - face-node connectivity it slices by -- again connectivity, not - coordinates. + One remains, from ``_slice_face_indices`` loading the connectivity it slices by. """ uxgrid = ux.open_grid(grid_path, chunks=CHUNKS) counter = _CountComputes() @@ -84,16 +69,10 @@ def track_computes_isel(self): class OpenGridChunked: - """``open_grid`` with ``chunks=``, across both oQU resolutions. - - Reads the file directly rather than through ``CachedFixtures``, because - reading the file is the subject here. The 120km mesh is ~16x the 480km one - (59,329 nodes against 3,947). + """``open_grid`` with ``chunks=`` across both oQU resolutions. - Track each resolution against its own history, not against the other. The - two files are different formats -- oQU480.grid.nc is netCDF4/HDF5, - oQU120.grid.nc is netCDF3 -- and the HDF5 open costs more, so 480km reads - *slower* than 120km despite a sixteenth of the data. + Compare each resolution to its own history: oQU480 is netCDF4/HDF5 and oQU120 + is netCDF3, so 480km opens slower despite a sixteenth of the data. """ param_names = ["resolution"] diff --git a/test/grid/grid/test_io.py b/test/grid/grid/test_io.py index 1c8b084c3..e840d05cb 100644 --- a/test/grid/grid/test_io.py +++ b/test/grid/grid/test_io.py @@ -7,27 +7,11 @@ def _assert_lon_close(actual, desired, err_msg, atol=ERROR_TOLERANCE): - """Compare longitudes as directions on a circle, to an absolute tolerance. - - Used for the Exodus round-trips, which are the ones that go - lon/lat -> xyz -> lon/lat and so come back through ``_xyz_to_lonlat_deg``. - ``assert_allclose(..., rtol=...)`` is the wrong instrument for that on two - counts. A longitude near zero has no magnitude for a relative tolerance to - be measured against -- 2.1e-14 vs 2.8e-14 is 7e-15 degrees apart and fails - at rtol=1e-8 (``outCSne30`` nodes 4372 and 4749). And longitude is - periodic: ``_xyz_to_lonlat_deg`` wraps into the half-open [-180, 180), - while ``_set_desired_longitude_range`` keeps both endpoints, so a node on - the antimeridian reads 180.0 on the original and -180.0 on the reload -- - the same meridian, scored as a 360-degree error (179 nodes of - ``outRLL1deg``, and ``outCSne30`` nodes 3966 and 5155, which sit one ulp - short of 180 and land on -180.0 once rounded through Cartesian). - - All of it was masked before. The old ``_set_desired_longitude_range`` - wrapped the whole array whenever any element exceeded 180, applying to the - original grid the identical ``(lon + 180) % 360 - 180`` that the reload - applies -- so both sides carried the same perturbation and the same - endpoint fold. The wrap is elementwise now and leaves in-range longitudes - alone, which leaves the round-trip's own error exposed. + """Compare longitudes modulo 360, to an absolute tolerance. + + Exodus round-trips pass through ``_xyz_to_lonlat_deg``, which returns -180.0 + where the original keeps 180.0, and near-zero longitudes come back ~1e-14 off, + which a relative tolerance alone rejects. """ diff = (np.asarray(actual) - np.asarray(desired) + 180.0) % 360.0 - 180.0 np.testing.assert_allclose(diff, 0.0, atol=atol, err_msg=err_msg) diff --git a/test/grid/grid/test_lon_range.py b/test/grid/grid/test_lon_range.py index 258815a76..3915f1b5b 100644 --- a/test/grid/grid/test_lon_range.py +++ b/test/grid/grid/test_lon_range.py @@ -1,13 +1,7 @@ -"""Guards on the [-180, 180] longitude wrap. +"""Tests for the elementwise longitude wrap in ``_set_desired_longitude_range``. -``_set_desired_longitude_range`` used to decide whether to wrap by asking -``lon.max() > 180``. On a dask-backed grid that reduction is a compute, and -``Grid.__init__`` calls it, so opening a grid pulled its longitude -coordinates into memory. The wrap is now elementwise. - -Two things have to hold for that to be an improvement rather than a trade: -the construction must stay lazy, and repeated calls must not pile ``where`` -layers onto the graph -- ``edge_lat`` invokes it on every property access. +It must stay lazy on dask-backed grids and must not grow the graph when called +repeatedly. """ import numpy as np @@ -45,12 +39,7 @@ def _shim(values, name="node_lon", **attrs): def test_wrap_is_elementwise_and_leaves_in_range_values_exact(): - """Outside [-180, 180] wraps; inside is passed through bit-for-bit. - - The old whole-array form perturbed in-range longitudes by up to ~3e-14 - degrees, because ``(lon + 180) - 180`` does not round-trip exactly. The - assertion here is ``assert_array_equal``, not ``allclose``, on purpose. - """ + """Outside [-180, 180] wraps; inside passes through bit-for-bit.""" rng = np.random.default_rng(0) values = np.concatenate( [rng.uniform(0.0, 180.0, 500), rng.uniform(180.0, 360.0, 500)] @@ -69,14 +58,10 @@ def test_wrap_is_elementwise_and_leaves_in_range_values_exact(): def test_both_endpoints_of_the_antimeridian_are_left_alone(): - """180.0 stays 180.0 and -180.0 stays -180.0, however the array looks. - - Not a detail. ``antimeridian_face_indices`` reads a face as crossing from - the span of its longitudes, so a face with one vertex at exactly 180 and - the rest near -170 is only visible while that vertex reads as 180 -- fold - it to -180 and the span drops from 350 to 10. Under the old reduction the - endpoint survived only by accident: an array whose maximum was exactly - 180 was not wrapped at all, so nothing touched it. + """180.0 and -180.0 are both kept, whatever else the array holds. + + ``antimeridian_face_indices`` detects a crossing from a face's longitude span, + which folding 180 to -180 would collapse. """ for companion in (10.0, 270.0, -170.0): grid = _shim(np.array([180.0, companion])) @@ -89,12 +74,7 @@ def test_both_endpoints_of_the_antimeridian_are_left_alone(): def test_wrap_normalizes_both_tails_and_is_idempotent(): - """One pass lands everything in [-180, 180); a second pass changes nothing. - - Both tails matter. The reduction tested ``lon.max() > 180``, so it reached - longitudes below -180 only when the same array happened to hold one above - 180, and left them alone otherwise. - """ + """One pass lands everything in [-180, 180], below -180 too; a second is a no-op.""" rng = np.random.default_rng(1) grid = _shim(rng.uniform(-720.0, 720.0, 1000)) @@ -127,12 +107,7 @@ def test_nan_longitudes_survive_the_wrap(): def test_open_grid_does_not_compute_longitudes(gridpath): - """The regression this change exists to fix. - - Counts every dask execution during ``open_grid``, then asserts the wrap - contributed none of them by re-running the wrap on the constructed grid - under its own counter. - """ + """The wrap runs no dask compute on a chunked grid, the regression fixed here.""" path = gridpath("ugrid", "outCSne30", "outCSne30.ug") grid = ux.open_grid(path, chunks={"n_node": 1000}) @@ -145,12 +120,7 @@ def test_open_grid_does_not_compute_longitudes(gridpath): def test_repeated_calls_do_not_grow_the_graph(gridpath): - """``edge_lat`` calls the wrap on every access, outside its populate guard. - - Without the memo, an unconditional elementwise wrap would add a ``where`` - layer per call -- lazy, so nothing would fail, just an ever-deepening - graph. - """ + """Repeated calls, as ``edge_lat`` makes, add no ``where`` layers to the graph.""" path = gridpath("ugrid", "outCSne30", "outCSne30.ug") grid = ux.open_grid(path, chunks={"n_node": 1000}) @@ -173,12 +143,7 @@ def test_memo_reopens_when_the_variable_is_replaced(gridpath): def test_eager_fast_path_agrees_with_the_where_element_for_element(): - """The in-memory guard must be a pure optimization. - - ``_lon_within_range`` lets the eager path skip the wrap entirely. That is - only sound if the wrap would have been the identity, so the two are - compared directly rather than the skip being assumed correct. - """ + """The in-memory guard must give exactly what the unguarded ``where`` gives.""" rng = np.random.default_rng(2) cases = { "strictly inside": rng.uniform(-179.0, 179.0, 500), @@ -201,11 +166,7 @@ def test_eager_fast_path_agrees_with_the_where_element_for_element(): def test_guard_is_skipped_for_dask_backed_arrays(gridpath): - """The guard is two reductions, which is exactly what must not run lazily. - - ``test_open_grid_does_not_compute_longitudes`` would catch this too, but - only as a compute count; this names the reason. - """ + """``_lon_within_range`` computes on dask, so the wrap only runs it eagerly.""" grid = ux.open_grid( gridpath("ugrid", "outCSne30", "outCSne30.ug"), chunks={"n_node": 1000} ) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index ed4bc92c6..c14f818d1 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -725,24 +725,11 @@ def _is_projected_grid(uxgrid) -> bool: def _lon_within_range(da) -> bool: - """Whether every element of ``da`` already lies in [-180, 180]. - - A guard, not a decision. When this is true the wrap below is the identity - on every element, so skipping it cannot change a value. - - It exists to keep the eager path off the ``where``. The wrap has to be - elementwise to stay out of a dask compute, but elementwise costs a pass - plus ~18 bytes per node of temporaries -- two bool masks, the arithmetic, - the result -- where a reduction costs a pass and allocates nothing. On an - array that is already in range, measured at 4M nodes, that is 20.9ms and - 72 MB to hand back the input unchanged, against 4.6ms and nothing. - Reductions are what made the old code compute, so this one is allowed - only where there is nothing to defer: ``da.chunks is None``. A - dask-backed array skips it and goes straight to the ``where``. - - ``and`` short-circuits, so an array that does need wrapping usually pays - a single ``max`` -- the same reduction the old code paid -- before - falling through. + """Whether all of ``da`` lies in [-180, 180], making the wrap the identity. + + Lets in-memory arrays skip the elementwise wrap, which costs ~5x a reduction + plus ~18 bytes per node of temporaries. Only call it when ``da.chunks is + None``: on a dask array these reductions are the compute being avoided. """ return bool(da.max() <= 180 and da.min() >= -180) @@ -750,45 +737,11 @@ def _lon_within_range(da) -> bool: def _set_desired_longitude_range(uxgrid): """Sets the longitude range to [-180, 180] for all longitude variables. - The wrap is elementwise rather than guarded by ``lon.max() > 180``. The - guard was a reduction, and a reduction on a dask-backed coordinate is a - compute -- so ``Grid.__init__``, which calls this, pulled every longitude - array into memory before the caller had asked for anything. The - elementwise form stays in the graph and splits over chunks. - - On the elements that need wrapping the two are bit-identical -- same - expression, same order. They differ on the elements that do not: the - reduction ran ``(lon + 180) % 360 - 180`` over the whole array once any - value exceeded 180, and that round-trip is not exact, so it perturbed - in-range longitudes by up to ~3e-14 degrees. ``xr.where`` passes those - through untouched. - - An element is wrapped when it falls outside [-180, 180] on *either* side. - The negative half of that predicate is what the reduction was missing: it - tested the maximum alone, so a longitude below -180 was normalized only - when the same array happened to also hold one above 180, and was left - where it was otherwise. - - Both endpoints are kept, rather than folding 180.0 onto -180.0 for a - half-open [-180, 180). ``_xyz_to_lonlat_deg`` does produce the half-open - interval, so a grid reloaded through Cartesian coordinates disagrees with - its original at a node sitting exactly on the antimeridian -- but that - node is the *point* of ``antimeridian_face_indices``, which reads a face - as crossing from the span of its longitudes. Folding 180.0 to -180.0 - collapses the span of a face that touches the antimeridian from the west - and hides it. The round-trip comparison is the cheaper of the two to make - periodic, and the Exodus tests do that. - - None of that applies to a grid whose coordinates are already in memory -- - there was never a compute to defer there, only a scan -- and elementwise - is the slower shape for it. ``_lon_within_range`` keeps the eager path off - the ``where`` when the wrap would be the identity anyway. - - Each variable is wrapped at most once. ``edge_lat`` calls this on every - access, outside its populate guard, so without the memo an unconditional - wrap would stack a ``where`` layer onto the graph per property access. Keying on the ``xr.Variable`` object makes the memo - self-invalidating: assigning into ``_ds`` replaces that object, so a - repopulated or user-assigned coordinate is wrapped again. + Wraps elementwise with ``xr.where``, so dask-backed coordinates stay lazy and + in-range values pass through bit-for-bit; both endpoints are kept because + ``antimeridian_face_indices`` relies on a vertex at 180. Each variable is + wrapped once, memoized on its ``xr.Variable``, since ``edge_lat`` calls this + on every access. Skipped entirely for projected grids: wrapping meter-scale coordinates as if they were degrees silently corrupts the geometry. A ``UserWarning``