Skip to content

Commit 52e767e

Browse files
committed
Fix silent lon wrapping for projected (non-spherical) UGRID grids
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.
1 parent ae51717 commit 52e767e

2 files changed

Lines changed: 137 additions & 1 deletion

File tree

test/grid/grid/test_projected.py

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
"""Tests for projected (non-spherical) coordinate detection and safe loading."""
2+
3+
import numpy as np
4+
import pytest
5+
import xarray as xr
6+
7+
import uxarray as ux
8+
9+
10+
def _make_grid(stdname, units, values=None):
11+
if values is None:
12+
values = np.array([1_500_000.0, 1_500_100.0, 1_500_100.0, 1_500_000.0])
13+
node_lon = xr.DataArray(
14+
values, dims=["n_node"], attrs={"standard_name": stdname, "units": units}
15+
)
16+
node_lat = xr.DataArray(
17+
np.array([800_000.0, 800_000.0, 800_100.0, 800_100.0]), dims=["n_node"]
18+
)
19+
fnc = xr.DataArray(
20+
np.array([[0, 1, 2], [0, 2, 3]]),
21+
dims=["n_face", "n_max_face_nodes"],
22+
attrs={"cf_role": "face_node_connectivity", "start_index": 0, "_FillValue": -1},
23+
)
24+
return xr.Dataset(
25+
{"node_lon": node_lon, "node_lat": node_lat, "face_node_connectivity": fnc}
26+
)
27+
28+
29+
def test_projected_coordinates_not_wrapped():
30+
"""Meter-scale coordinates with standard_name=projection_x_coordinate must not be wrapped."""
31+
original = np.array([1_500_000.0, 1_500_100.0, 1_500_100.0, 1_500_000.0])
32+
ds = _make_grid("projection_x_coordinate", "m")
33+
with pytest.warns(UserWarning, match="Projected"):
34+
grid = ux.Grid(ds, source_grid_spec="UGRID")
35+
np.testing.assert_array_equal(grid.node_lon.values, original)
36+
37+
38+
def test_projected_detected_from_units():
39+
"""units='m' with no standard_name is sufficient to detect projected coords."""
40+
ds = _make_grid("", "m")
41+
ds["node_lon"].attrs.pop("standard_name", None)
42+
with pytest.warns(UserWarning, match="Projected"):
43+
ux.Grid(ds, source_grid_spec="UGRID")
44+
45+
46+
def test_projected_detected_from_grid_mapping():
47+
"""A grid_mapping variable with a non-latlon name signals projected coords."""
48+
ds = _make_grid("", "")
49+
ds["node_lon"].attrs = {"grid_mapping": "crs"}
50+
ds["crs"] = xr.DataArray(
51+
np.int32(0), attrs={"grid_mapping_name": "lambert_conformal_conic"}
52+
)
53+
with pytest.warns(UserWarning, match="Projected"):
54+
ux.Grid(ds, source_grid_spec="UGRID")
55+
56+
57+
def test_geographic_coordinates_still_wrapped():
58+
"""Geographic [0, 360] coords must still be normalized to [-180, 180]."""
59+
ds = _make_grid(
60+
"longitude",
61+
"degrees_east",
62+
values=np.array([270.0, 271.0, 271.0, 270.0]),
63+
)
64+
ds["node_lat"] = xr.DataArray(np.array([43.0, 43.0, 44.0, 44.0]), dims=["n_node"])
65+
grid = ux.Grid(ds, source_grid_spec="UGRID")
66+
assert grid.node_lon.values.max() <= 180.0

uxarray/grid/coordinates.py

Lines changed: 71 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -681,8 +681,78 @@ def _construct_edge_centroids(node_x, node_y, node_z, edge_node_conn):
681681
return _normalize_xyz(centroid_x, centroid_y, centroid_z)
682682

683683

684+
def _is_projected_grid(uxgrid) -> bool:
685+
"""Return True if the grid uses projected (non-spherical) coordinates.
686+
687+
Detection priority (mirrors CF conventions):
688+
689+
1. ``standard_name = "projection_x_coordinate"`` on ``node_lon`` — the
690+
clearest CF signal used by regional ocean and coastal models.
691+
2. ``units`` on ``node_lon`` is a length unit (m, km, ft, …) rather than
692+
angular — catches files that omit ``standard_name`` but carry units.
693+
3. A ``grid_mapping`` attribute on ``node_lon`` pointing to a variable
694+
whose ``grid_mapping_name`` is not ``"latitude_longitude"`` — handles
695+
files that embed a full CRS variable (e.g. Lambert Conformal, Albers).
696+
697+
Falls back to ``False`` (geographic) when no signal is found, preserving
698+
backward compatibility with all existing geographic UGRID files.
699+
"""
700+
if "node_lon" not in uxgrid._ds:
701+
return False
702+
703+
da = uxgrid._ds["node_lon"]
704+
attrs = da.attrs
705+
706+
# 1. standard_name
707+
stdname = attrs.get("standard_name", "")
708+
if stdname == "projection_x_coordinate":
709+
return True
710+
if stdname in ("longitude",):
711+
return False
712+
713+
# 2. units — angular units mean geographic
714+
units = attrs.get("units", "").lower().strip()
715+
_ANGULAR_UNITS = {"degrees_east", "degrees", "degree", "deg", "rad", "radians"}
716+
_LENGTH_UNITS = {"m", "km", "meter", "meters", "metre", "metres", "ft", "feet"}
717+
if units in _LENGTH_UNITS:
718+
return True
719+
if units in _ANGULAR_UNITS:
720+
return False
721+
722+
# 3. grid_mapping variable
723+
gm_name = attrs.get("grid_mapping", "")
724+
if gm_name and gm_name in uxgrid._ds:
725+
gm_attrs = uxgrid._ds[gm_name].attrs
726+
gm_name_val = gm_attrs.get("grid_mapping_name", "latitude_longitude")
727+
if gm_name_val != "latitude_longitude":
728+
return True
729+
730+
return False
731+
732+
684733
def _set_desired_longitude_range(uxgrid):
685-
"""Sets the longitude range to [-180, 180] for all longitude variables."""
734+
"""Sets the longitude range to [-180, 180] for all longitude variables.
735+
736+
Skipped entirely for projected grids: wrapping meter-scale coordinates
737+
as if they were degrees silently corrupts the geometry. A ``UserWarning``
738+
is issued on the first access so users know which operations are invalid.
739+
"""
740+
if _is_projected_grid(uxgrid):
741+
import warnings
742+
743+
warnings.warn(
744+
"Projected (non-spherical) coordinates detected on this grid "
745+
"(standard_name='projection_x_coordinate' or length units). "
746+
"UXarray's geometry algorithms (face areas, GCA intersections, "
747+
"zonal mean, bounds, cross-sections) assume a unit sphere and "
748+
"will produce incorrect results on projected grids. "
749+
"Loading, plotting, and connectivity operations are unaffected. "
750+
"Check Grid.is_projected for the detected coordinate type.",
751+
UserWarning,
752+
stacklevel=3,
753+
)
754+
return
755+
686756
with xr.set_options(keep_attrs=True):
687757
for lon_name in ["node_lon", "edge_lon", "face_lon"]:
688758
if lon_name in uxgrid._ds:

0 commit comments

Comments
 (0)