Fix MPAS reader bug: Normalize after loading (#1344)

* o MPAS reader bug: Normalize after loading

* Add test for MPAS grid normalization verification

* Add face area and edge length normalization tests for unit sphere

* o Fix area normalization in MPAS reader

* o MPAS files might not have the sphere_radius attribute fix reader

* io/mpas: guard normalization of dvEdge/dcEdge by radius!=1.0 for consistency with XYZ/areas; no functional change when unit sphere. Verified by test/test_mpas.py (11 passed).

---------

Co-authored-by: Rajeev Jain <rajeeja@@gmail.com>
Co-authored-by: Philip Chmielowiec <[email protected]>

Author: 3 people

test/test_mpas.py

@@ -133,3 +133,100 @@ def test_distance_units():
 
     nt.assert_array_almost_equal(uxgrid['edge_node_distances'].values, (xrds['dvEdge'].values / xrds.attrs['sphere_radius']))
     nt.assert_array_almost_equal(uxgrid['edge_face_distances'].values, (xrds['dcEdge'].values / xrds.attrs['sphere_radius']))
+
+
+def test_ocean_mesh_normalization():
+    """Test that MPAS ocean mesh with non-unit sphere radius is properly normalized."""
+    # Ocean mesh has sphere_radius = 6371229.0 meters
+    uxgrid = ux.open_grid(mpas_ocean_mesh, use_dual=False)
+
+    # Check node coordinates are normalized to unit sphere
+    node_x = uxgrid._ds['node_x'].values
+    node_y = uxgrid._ds['node_y'].values
+    node_z = uxgrid._ds['node_z'].values
+
+    # Compute radius for each node
+    radii = np.sqrt(node_x**2 + node_y**2 + node_z**2)
+
+    # All radii should be 1.0 (unit sphere)
+    nt.assert_array_almost_equal(radii, np.ones_like(radii), decimal=10)
+
+    # Check face areas are normalized for unit sphere
+    # Ocean mesh only covers ~70% of sphere (ocean area)
+    face_areas = uxgrid.face_areas.values
+
+    # Check that all face areas are positive
+    assert np.all(face_areas > 0), "All face areas should be positive"
+
+    # Total area should be less than full sphere (4*pi) but reasonable
+    total_area = np.sum(face_areas)
+    full_sphere_area = 4.0 * np.pi
+    assert 0.5 < total_area < full_sphere_area, f"Total area {total_area} should be between 0.5 and {full_sphere_area}"
+
+    # Check that individual face areas are reasonable for unit sphere
+    max_face_area = np.max(face_areas)
+    min_face_area = np.min(face_areas)
+    assert max_face_area < 0.1 * full_sphere_area, "Maximum face area seems too large for unit sphere"
+    assert min_face_area > 1e-10, "Minimum face area seems too small"
+
+    # Check edge lengths are normalized for unit sphere
+    if "edge_node_distances" in uxgrid._ds:
+        edge_lengths = uxgrid._ds["edge_node_distances"].values
+
+        # All edge lengths should be positive and <= pi
+        assert np.all(edge_lengths > 0), "All edge lengths should be positive"
+        assert np.max(edge_lengths) <= np.pi, "Edge lengths should not exceed pi on unit sphere"
+
+
+def test_grid_normalization():
+    """Test that MPAS grid coordinates are properly normalized."""
+    uxgrid = ux.open_grid(mpas_grid_path, use_dual=False)
+
+    # Check node coordinates are normalized
+    node_lon = uxgrid._ds['node_lon'].values
+    node_lat = uxgrid._ds['node_lat'].values
+    node_x = uxgrid._ds['node_x'].values
+    node_y = uxgrid._ds['node_y'].values
+    node_z = uxgrid._ds['node_z'].values
+
+    # Compute radius for each node
+    radii = np.sqrt(node_x**2 + node_y**2 + node_z**2)
+
+    # All radii should be 1.0 (unit sphere)
+    nt.assert_array_almost_equal(radii, np.ones_like(radii), decimal=10)
+
+    # Check face areas are normalized for unit sphere
+    # Total surface area of unit sphere is 4*pi
+    expected_total_area = 4.0 * np.pi
+
+    # Get face areas
+    face_areas = uxgrid.face_areas.values
+
+    # Check that all face areas are positive
+    assert np.all(face_areas > 0), "All face areas should be positive"
+
+    # Check that total area equals surface area of unit sphere
+    total_area = np.sum(face_areas)
+    nt.assert_almost_equal(total_area, expected_total_area, decimal=7)
+
+    # Check that face areas are reasonable (not too small or too large)
+    # For a unit sphere, face areas should be much smaller than total area
+    max_face_area = np.max(face_areas)
+    min_face_area = np.min(face_areas)
+
+    assert max_face_area < 0.1 * expected_total_area, "Maximum face area seems too large for unit sphere"
+    assert min_face_area > 1e-10, "Minimum face area seems too small"
+
+    # Check edge lengths are normalized for unit sphere
+    # Edge lengths should be reasonable for a unit sphere (max possible is pi for antipodal points)
+    if "edge_node_distances" in uxgrid._ds:
+        edge_lengths = uxgrid._ds["edge_node_distances"].values
+
+        # All edge lengths should be positive
+        assert np.all(edge_lengths > 0), "All edge lengths should be positive"
+
+        # Maximum edge length on unit sphere cannot exceed pi
+        assert np.max(edge_lengths) <= np.pi, "Edge lengths should not exceed pi on unit sphere"
+
+        # Minimum edge length should be reasonable
+        assert np.min(edge_lengths) > 1e-10, "Minimum edge length seems too small"

uxarray/io/_mpas.py

@@ -190,6 +190,13 @@ def _parse_node_xyz_coords(in_ds, out_ds, mesh_type):
         node_y = node_y.rename({"nCells": ugrid.NODE_DIM})
         node_z = node_z.rename({"nCells": ugrid.NODE_DIM})
 
+    # Normalize coordinates to unit sphere if needed
+    radius = in_ds.attrs.get("sphere_radius", 1.0)
+    if radius != 1.0:
+        node_x = node_x / radius
+        node_y = node_y / radius
+        node_z = node_z / radius
+
     out_ds["node_x"] = node_x.assign_attrs(ugrid.NODE_X_ATTRS)
     out_ds["node_y"] = node_y.assign_attrs(ugrid.NODE_Y_ATTRS)
     out_ds["node_z"] = node_z.assign_attrs(ugrid.NODE_Z_ATTRS)
@@ -237,6 +244,13 @@ def _parse_face_xyz_coords(in_ds, out_ds, mesh_type):
         face_y = face_y.rename({"nVertices": ugrid.FACE_DIM})
         face_z = face_z.rename({"nVertices": ugrid.FACE_DIM})
 
+    # Normalize coordinates to unit sphere if needed
+    radius = in_ds.attrs.get("sphere_radius", 1.0)
+    if radius != 1.0:
+        face_x = face_x / radius
+        face_y = face_y / radius
+        face_z = face_z / radius
+
     out_ds["face_x"] = face_x.assign_attrs(ugrid.FACE_X_ATTRS)
     out_ds["face_y"] = face_y.assign_attrs(ugrid.FACE_Y_ATTRS)
     out_ds["face_z"] = face_z.assign_attrs(ugrid.FACE_Z_ATTRS)
@@ -266,6 +280,13 @@ def _parse_edge_xyz_coords(in_ds, out_ds, mesh_type):
     edge_y = edge_y.rename({"nEdges": ugrid.EDGE_DIM})
     edge_z = edge_z.rename({"nEdges": ugrid.EDGE_DIM})
 
+    # Normalize coordinates to unit sphere if needed
+    radius = in_ds.attrs.get("sphere_radius", 1.0)
+    if radius != 1.0:
+        edge_x = edge_x / radius
+        edge_y = edge_y / radius
+        edge_z = edge_z / radius
+
     out_ds["edge_x"] = edge_x.assign_attrs(ugrid.EDGE_X_ATTRS)
     out_ds["edge_y"] = edge_y.assign_attrs(ugrid.EDGE_Y_ATTRS)
     out_ds["edge_z"] = edge_z.assign_attrs(ugrid.EDGE_Z_ATTRS)
@@ -380,7 +401,10 @@ def _parse_face_faces(in_ds, out_ds, mesh_type):
 
 def _parse_edge_node_distances(in_ds, out_ds):
     """Parses ``edge_node_distances``"""
-    edge_node_distances = in_ds["dvEdge"] / in_ds.attrs["sphere_radius"]
+    radius = in_ds.attrs.get("sphere_radius", 1.0)
+    edge_node_distances = in_ds["dvEdge"]
+    if radius != 1.0:
+        edge_node_distances = edge_node_distances / radius
 
     out_ds["edge_node_distances"] = edge_node_distances.assign_attrs(
         descriptors.EDGE_NODE_DISTANCES_ATTRS
@@ -389,7 +413,10 @@ def _parse_edge_node_distances(in_ds, out_ds):
 
 def _parse_edge_face_distances(in_ds, out_ds):
     """Parses ``edge_face_distances``"""
-    edge_face_distances = in_ds["dcEdge"] / in_ds.attrs["sphere_radius"]
+    radius = in_ds.attrs.get("sphere_radius", 1.0)
+    edge_face_distances = in_ds["dcEdge"]
+    if radius != 1.0:
+        edge_face_distances = edge_face_distances / radius
 
     out_ds["edge_face_distances"] = edge_face_distances.assign_attrs(
         descriptors.EDGE_FACE_DISTANCES_ATTRS
@@ -408,6 +435,12 @@ def _parse_face_areas(in_ds, out_ds, mesh_type):
     else:
         face_area = in_ds["areaTriangle"]
 
+    # Normalize face areas to unit sphere if needed
+    radius = in_ds.attrs.get("sphere_radius", 1.0)
+    if radius != 1.0:
+        # Area scales with radius squared
+        face_area = face_area / (radius * radius)
+
     out_ds["face_areas"] = face_area.assign_attrs(descriptors.FACE_AREAS_ATTRS).rename(
         {face_area.dims[0]: ugrid.FACE_DIM}
     )