From 52e767e96cd5f5edc99fc75b8b02eba7f55f91af Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Fri, 12 Jun 2026 16:23:21 -0500 Subject: [PATCH 1/2] Fix silent lon wrapping for projected (non-spherical) UGRID grids MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Fixes #1524. _is_projected_grid() in coordinates.py detects projected coordinates from CF standard_name ("projection_x_coordinate"), length units (m, km, ...), or a non-latlon grid_mapping variable — in that priority order, falling back to geographic (backward-compatible). _set_desired_longitude_range skips wrapping and emits a UserWarning when projected coordinates are detected, naming which spherical operations are invalid. Coordinates are preserved as-is. 4 tests: no-wrap on load, units detection, grid_mapping detection, geographic wrapping still works. --- test/grid/grid/test_projected.py | 66 +++++++++++++++++++++++++++++ uxarray/grid/coordinates.py | 72 +++++++++++++++++++++++++++++++- 2 files changed, 137 insertions(+), 1 deletion(-) create mode 100644 test/grid/grid/test_projected.py diff --git a/test/grid/grid/test_projected.py b/test/grid/grid/test_projected.py new file mode 100644 index 000000000..d9172d239 --- /dev/null +++ b/test/grid/grid/test_projected.py @@ -0,0 +1,66 @@ +"""Tests for projected (non-spherical) coordinate detection and safe loading.""" + +import numpy as np +import pytest +import xarray as xr + +import uxarray as ux + + +def _make_grid(stdname, units, values=None): + if values is None: + values = np.array([1_500_000.0, 1_500_100.0, 1_500_100.0, 1_500_000.0]) + node_lon = xr.DataArray( + values, dims=["n_node"], attrs={"standard_name": stdname, "units": units} + ) + node_lat = xr.DataArray( + np.array([800_000.0, 800_000.0, 800_100.0, 800_100.0]), dims=["n_node"] + ) + fnc = xr.DataArray( + np.array([[0, 1, 2], [0, 2, 3]]), + dims=["n_face", "n_max_face_nodes"], + attrs={"cf_role": "face_node_connectivity", "start_index": 0, "_FillValue": -1}, + ) + return xr.Dataset( + {"node_lon": node_lon, "node_lat": node_lat, "face_node_connectivity": fnc} + ) + + +def test_projected_coordinates_not_wrapped(): + """Meter-scale coordinates with standard_name=projection_x_coordinate must not be wrapped.""" + original = np.array([1_500_000.0, 1_500_100.0, 1_500_100.0, 1_500_000.0]) + ds = _make_grid("projection_x_coordinate", "m") + with pytest.warns(UserWarning, match="Projected"): + grid = ux.Grid(ds, source_grid_spec="UGRID") + np.testing.assert_array_equal(grid.node_lon.values, original) + + +def test_projected_detected_from_units(): + """units='m' with no standard_name is sufficient to detect projected coords.""" + ds = _make_grid("", "m") + ds["node_lon"].attrs.pop("standard_name", None) + with pytest.warns(UserWarning, match="Projected"): + ux.Grid(ds, source_grid_spec="UGRID") + + +def test_projected_detected_from_grid_mapping(): + """A grid_mapping variable with a non-latlon name signals projected coords.""" + ds = _make_grid("", "") + ds["node_lon"].attrs = {"grid_mapping": "crs"} + ds["crs"] = xr.DataArray( + np.int32(0), attrs={"grid_mapping_name": "lambert_conformal_conic"} + ) + with pytest.warns(UserWarning, match="Projected"): + ux.Grid(ds, source_grid_spec="UGRID") + + +def test_geographic_coordinates_still_wrapped(): + """Geographic [0, 360] coords must still be normalized to [-180, 180].""" + ds = _make_grid( + "longitude", + "degrees_east", + values=np.array([270.0, 271.0, 271.0, 270.0]), + ) + ds["node_lat"] = xr.DataArray(np.array([43.0, 43.0, 44.0, 44.0]), dims=["n_node"]) + grid = ux.Grid(ds, source_grid_spec="UGRID") + assert grid.node_lon.values.max() <= 180.0 diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 0e97bcf5a..6870bbdcc 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -681,8 +681,78 @@ def _construct_edge_centroids(node_x, node_y, node_z, edge_node_conn): return _normalize_xyz(centroid_x, centroid_y, centroid_z) +def _is_projected_grid(uxgrid) -> bool: + """Return True if the grid uses projected (non-spherical) coordinates. + + Detection priority (mirrors CF conventions): + + 1. ``standard_name = "projection_x_coordinate"`` on ``node_lon`` — the + clearest CF signal used by regional ocean and coastal models. + 2. ``units`` on ``node_lon`` is a length unit (m, km, ft, …) rather than + angular — catches files that omit ``standard_name`` but carry units. + 3. A ``grid_mapping`` attribute on ``node_lon`` pointing to a variable + whose ``grid_mapping_name`` is not ``"latitude_longitude"`` — handles + files that embed a full CRS variable (e.g. Lambert Conformal, Albers). + + Falls back to ``False`` (geographic) when no signal is found, preserving + backward compatibility with all existing geographic UGRID files. + """ + if "node_lon" not in uxgrid._ds: + return False + + da = uxgrid._ds["node_lon"] + attrs = da.attrs + + # 1. standard_name + stdname = attrs.get("standard_name", "") + if stdname == "projection_x_coordinate": + return True + if stdname in ("longitude",): + return False + + # 2. units — angular units mean geographic + units = attrs.get("units", "").lower().strip() + _ANGULAR_UNITS = {"degrees_east", "degrees", "degree", "deg", "rad", "radians"} + _LENGTH_UNITS = {"m", "km", "meter", "meters", "metre", "metres", "ft", "feet"} + if units in _LENGTH_UNITS: + return True + if units in _ANGULAR_UNITS: + return False + + # 3. grid_mapping variable + gm_name = attrs.get("grid_mapping", "") + if gm_name and gm_name in uxgrid._ds: + gm_attrs = uxgrid._ds[gm_name].attrs + gm_name_val = gm_attrs.get("grid_mapping_name", "latitude_longitude") + if gm_name_val != "latitude_longitude": + return True + + return False + + def _set_desired_longitude_range(uxgrid): - """Sets the longitude range to [-180, 180] for all longitude variables.""" + """Sets the longitude range to [-180, 180] for all longitude variables. + + Skipped entirely for projected grids: wrapping meter-scale coordinates + as if they were degrees silently corrupts the geometry. A ``UserWarning`` + is issued on the first access so users know which operations are invalid. + """ + if _is_projected_grid(uxgrid): + import warnings + + warnings.warn( + "Projected (non-spherical) coordinates detected on this grid " + "(standard_name='projection_x_coordinate' or length units). " + "UXarray's geometry algorithms (face areas, GCA intersections, " + "zonal mean, bounds, cross-sections) assume a unit sphere and " + "will produce incorrect results on projected grids. " + "Loading, plotting, and connectivity operations are unaffected. " + "Check Grid.is_projected for the detected coordinate type.", + UserWarning, + stacklevel=3, + ) + return + with xr.set_options(keep_attrs=True): for lon_name in ["node_lon", "edge_lon", "face_lon"]: if lon_name in uxgrid._ds: From 22c826f4fbf77a57a17675657e0ee46984b9fd7a Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Fri, 12 Jun 2026 16:35:30 -0500 Subject: [PATCH 2/2] docs: document projected coordinate support in grid-formats --- docs/user-guide/grid-formats.rst | 58 ++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/docs/user-guide/grid-formats.rst b/docs/user-guide/grid-formats.rst index 05cb8eb4e..4b89844ca 100644 --- a/docs/user-guide/grid-formats.rst +++ b/docs/user-guide/grid-formats.rst @@ -46,6 +46,64 @@ References * https://cf-convention.github.io/Data/cf-conventions/cf-conventions-1.13/cf-conventions.html#ugrid-conventions * https://cf-convention.github.io/Data/cf-conventions/cf-conventions-1.13/cf-conventions.html#mesh-topology-variables +Projected (Non-Spherical) Coordinates +-------------------------------------- + +UGRID does not assume a sphere — some regional and coastal models write UGRID +files with **projected coordinates** (e.g. a Lambert Conformal or Albers Equal +Area projection) in units of meters rather than degrees. UXarray detects this +automatically from the CF metadata already present on the coordinate variables: + +.. list-table:: + :header-rows: 1 + :widths: 40 60 + + * - CF signal on ``node_lon`` + - Detected as + * - ``standard_name = "projection_x_coordinate"`` + - Projected + * - ``standard_name = "longitude"`` + - Geographic + * - ``units`` is ``m``, ``km``, ``ft``, … + - Projected + * - ``units`` is ``degrees_east``, ``deg``, … + - Geographic + * - ``grid_mapping`` variable with non-latlon ``grid_mapping_name`` + - Projected + * - No recognizable signal + - Geographic (backward-compatible default) + +When projected coordinates are detected, UXarray: + +* **Preserves coordinates as-is** — no longitude wrapping is applied. +* **Emits a** ``UserWarning`` **on load** listing which operations are invalid. + +**What works on projected grids:** loading, data access, connectivity +traversal, and plotting (coordinates are passed through as-is to the plotting +backend). + +**What does not work:** UXarray's geometry algorithms assume a unit sphere — +face areas, GCA intersections, zonal mean, face bounds, and cross-sections will +produce incorrect results and are not supported on projected grids. + +.. code-block:: python + + import xarray as xr + import uxarray as ux + + # Load a projected UGRID file — coordinates are preserved, warning is issued + import warnings + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + grid = ux.open_grid("projected_grid.nc") + if w: + print(w[0].message) + +.. note:: + ``Grid.from_topology`` accepts raw NumPy arrays with no CF attributes, so + projection cannot be auto-detected. Use ``Grid.from_dataset`` with a + properly attributed ``xr.Dataset`` for projected grids. + MPAS ====