Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 58 additions & 0 deletions docs/user-guide/grid-formats.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
====

Expand Down
66 changes: 66 additions & 0 deletions test/grid/grid/test_projected.py
Original file line number Diff line number Diff line change
@@ -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
72 changes: 71 additions & 1 deletion uxarray/grid/coordinates.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"}
Comment on lines +714 to +716
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
Comment on lines +740 to +754

with xr.set_options(keep_attrs=True):
for lon_name in ["node_lon", "edge_lon", "face_lon"]:
if lon_name in uxgrid._ds:
Expand Down
Loading