Raise the underlying cause when geospatial data cannot be read

ci/environment.yml

@@ -32,6 +32,6 @@ dependencies:
   - shapely
   - spatialpandas
   - pooch
-  - geopandas
+  - geopandas>=1.0
   - xarray
   - asv

pyproject.toml

@@ -40,7 +40,7 @@ dependencies = [
   "scipy",
   "shapely",
   "spatialpandas",
-  "geopandas",
+  "geopandas>=1.0",  # 1.0 is the first release requiring pyogrio, the file-IO backend read_file needs. See #1693.
   "xarray",
   "hvplot",
   "healpix",

test/core/test_api.py

@@ -218,6 +218,21 @@ def mock_open_dataset(*args, **kwargs):
             os.unlink(tmp_path)
 
 
+def test_open_dataset_with_fallback_chains_both_engine_errors(tmp_path):
+    """When both engines fail, the fallback error must be chained onto the
+    default engine's error rather than replacing it."""
+
+    not_netcdf = tmp_path / "not_netcdf.nc"
+    not_netcdf.write_text("this is not a netcdf file")
+
+    with pytest.raises(Exception) as excinfo:
+        _open_dataset_with_fallback(str(not_netcdf))
+
+    assert excinfo.value.__cause__ is not None, (
+        "the default engine's error was discarded"
+    )
+
+
 def test_list_grid_names_multigrid(gridpath):
     """List grids from an OASIS-style multi-grid file."""
     grid_file = gridpath("scrip", "oasis", "grids.nc")

test/io/test_fesom.py

@@ -4,6 +4,7 @@
 import pytest
 
 import uxarray as ux
+from uxarray.io._fesom2 import _parse_elem2d, _parse_nod2d
 
 
 
@@ -62,3 +63,15 @@ def test_open_mfdataset_pi_path(test_data_dir):
     assert "n_node" in uxds.dims
     assert "n_face" in uxds.dims
     assert len(uxds) == 3
+
+
+def test_parse_nod2d_missing_file_names_the_file(tmp_path):
+    """A missing 'nod2d.out' must say which file is missing and where."""
+    with pytest.raises(FileNotFoundError, match="nod2d.out"):
+        _parse_nod2d(str(tmp_path))
+
+
+def test_parse_elem2d_missing_file_names_the_file(tmp_path):
+    """A missing 'elem2d.out' must say which file is missing and where."""
+    with pytest.raises(FileNotFoundError, match="elem2d.out"):
+        _parse_elem2d(str(tmp_path))

test/io/test_geopandas.py

@@ -1,5 +1,10 @@
+import geopandas as gpd
 import numpy as np
+import pytest
+from shapely.geometry import Point, Polygon
+
 import uxarray as ux
+from uxarray.io._geopandas import _extract_geometry_info, _gpd_read, _set_crs
 
 def test_read_shpfile(test_data_dir):
     """Read a shapefile."""
@@ -40,3 +45,45 @@ def test_load_xarray_with_from_file(gridpath):
     nc_filename = gridpath("scrip", "outCSne8", "outCSne8.nc")
     uxgrid = ux.Grid.from_file(nc_filename, backend="xarray")
     uxgrid.validate()
+
+
+def test_read_failure_raises(tmp_path):
+    """A read failure must surface the backend's own error rather than being
+    printed and swallowed into an UnboundLocalError.
+
+    Regression test for issue #1693.
+    """
+    not_geospatial = tmp_path / "not_geospatial.shp"
+    not_geospatial.write_text("this is not a shapefile")
+
+    with pytest.raises(Exception) as excinfo:
+        _gpd_read(str(not_geospatial))
+
+    assert not isinstance(excinfo.value, UnboundLocalError)
+
+
+def test_set_crs_warns_when_crs_is_missing():
+    """Assuming WGS84 for CRS-less data is a guess and must be announced."""
+    gdf = gpd.GeoDataFrame(
+        geometry=[Polygon([(0, 0), (1, 0), (1, 1)])], crs=None
+    )
+
+    with pytest.warns(UserWarning, match="no CRS"):
+        out = _set_crs(gdf)
+
+    assert out.crs is not None
+
+
+def test_unsupported_geometry_is_reported():
+    """Dropping a geometry silently would yield a grid missing a face with no
+    indication that anything was skipped."""
+    gdf = gpd.GeoDataFrame(
+        geometry=[Polygon([(0, 0), (1, 0), (1, 1), (0, 0)]), Point(5, 5)],
+        crs="EPSG:4326",
+    )
+
+    with pytest.warns(UserWarning, match="unsupported geometry type"):
+        node_lon, node_lat, connectivity = _extract_geometry_info(gdf, 4)
+
+    # Only the polygon contributes a face; the point is skipped.
+    assert connectivity.shape[0] == 1

test/io/test_utils.py

@@ -3,7 +3,7 @@
 import xarray as xr
 
 from uxarray.errors import GridInvalidError
-from uxarray.io.utils import _parse_grid_type
+from uxarray.io.utils import _is_structured, _parse_grid_type
 
 
 @pytest.mark.parametrize(
@@ -64,3 +64,42 @@ def test_parse_grid_type_detects_structured_grid():
 def test_parse_grid_type_rejects_incomplete_format_signals(dataset):
     with pytest.raises(GridInvalidError, match="Failed to parse uxgrid information from xarray.Dataset."):
         _parse_grid_type(dataset)
+
+
+def test_parse_grid_type_is_quiet_for_non_structured_grids(capsys):
+    """`_is_structured` runs before every other format check, so an unstructured
+    grid carrying lat/lon coordinates must not produce spurious output."""
+    lat = xr.DataArray(
+        np.array([[0.0, 1.0], [2.0, 3.0]]),
+        dims=["y", "x"],
+        attrs={"standard_name": "latitude"},
+    )
+    lon = xr.DataArray(
+        np.array([[0.0, 1.0], [2.0, 3.0]]),
+        dims=["y", "x"],
+        attrs={"standard_name": "longitude"},
+    )
+    ds = xr.Dataset(coords={"lat": lat, "lon": lon})
+
+    with pytest.raises(GridInvalidError):
+        _parse_grid_type(ds)
+
+    assert capsys.readouterr().out == ""
+
+
+def test_parse_grid_type_is_quiet_for_irregular_spacing():
+    """Irregularly spaced coordinates are simply 'not structured', not an event
+    worth reporting to stdout."""
+    lon = xr.DataArray(
+        np.array([0.0, 1.0, 4.0]),
+        dims=["lon"],
+        attrs={"standard_name": "longitude"},
+    )
+    lat = xr.DataArray(
+        np.array([-1.0, 0.0, 1.0]),
+        dims=["lat"],
+        attrs={"standard_name": "latitude"},
+    )
+    structured, _, _ = _is_structured(xr.Dataset(coords={"lon": lon, "lat": lat}))
+
+    assert not structured

uxarray/core/utils.py

@@ -32,11 +32,18 @@ def _open_dataset_with_fallback(filename_or_obj, chunks=None, **kwargs):
     try:
         # Try opening with xarray's default read engine
         return xr.open_dataset(filename_or_obj, chunks=chunks, **kwargs)
-    except Exception:
+    except Exception as default_engine_error:
         # If it fails, use the "netcdf4" engine as backup
         # Extract engine from kwargs to prevent duplicate parameter error
         engine = kwargs.pop("engine", "netcdf4")
-        return xr.open_dataset(filename_or_obj, engine=engine, chunks=chunks, **kwargs)
+        try:
+            return xr.open_dataset(
+                filename_or_obj, engine=engine, chunks=chunks, **kwargs
+            )
+        except Exception as fallback_error:
+            # Chain the fallback onto the original so both engines' reasons are
+            # visible; otherwise the default engine's error is lost entirely.
+            raise fallback_error from default_engine_error
 
 
 def _map_dims_to_ugrid(

uxarray/io/_fesom2.py

@@ -89,7 +89,10 @@ def _parse_nod2d(grid_path):
 
     file_path = os.path.join(grid_path, "nod2d.out")
     if not os.path.isfile(file_path):
-        raise FileNotFoundError("TODO: ")
+        raise FileNotFoundError(
+            f"Expected a FESOM2 ASCII grid directory containing 'nod2d.out', "
+            f"but no such file exists under {os.path.abspath(grid_path)!r}."
+        )
 
     nodes = pd.read_csv(
         file_path,
@@ -120,7 +123,10 @@ def _parse_elem2d(grid_path):
     """
     file_path = os.path.join(grid_path, "elem2d.out")
     if not os.path.isfile(file_path):
-        raise FileNotFoundError("TODO: ")
+        raise FileNotFoundError(
+            f"Expected a FESOM2 ASCII grid directory containing 'elem2d.out', "
+            f"but no such file exists under {os.path.abspath(grid_path)!r}."
+        )
 
     file_content = pd.read_csv(
         file_path,

uxarray/io/_geopandas.py

@@ -1,3 +1,5 @@
+import warnings
+
 import numpy as np
 import xarray as xr
 
@@ -66,11 +68,8 @@ def _gpd_read(filepath, driver=None, **kwargs):
 
     import geopandas as gpd
 
-    try:
-        gdf = gpd.read_file(filepath, driver=driver, **kwargs)
-        gdf = _set_crs(gdf)
-    except Exception as e:
-        print(f"An error occurred while reading the geospatial data: {e}")
+    gdf = gpd.read_file(filepath, driver=driver, **kwargs)
+    gdf = _set_crs(gdf)
 
     max_polygon_nodes = gdf["geometry"].apply(_get_num_nodes).max()
 
@@ -92,11 +91,14 @@ def _set_crs(gdf):
     """
     if gdf.crs is None:
         gdf = gdf.set_crs(WGS84_CRS)
-        print("Original CRS: None\nAssigned CRS:", gdf.crs)
+        warnings.warn(
+            f"The geospatial data declares no CRS; assuming {WGS84_CRS}. "
+            f"Coordinates will be wrong if the source uses a different CRS.",
+            stacklevel=2,
+        )
 
     if gdf.crs != WGS84_CRS:
         gdf = gdf.to_crs(WGS84_CRS)
-        print("Transformed CRS:", gdf.crs)
 
     return gdf
 
@@ -137,7 +139,14 @@ def _extract_geometry_info(gdf, max_coord_size):
                 geometry, node_lat_list, node_lon_list, connectivity, node_index
             )
         else:
-            print(f"Unsupported geometry type: {geometry.geom_type}")
+            # Skipping a geometry silently would yield a grid that is missing
+            # faces without any indication that data was dropped.
+            warnings.warn(
+                f"Skipping unsupported geometry type {geometry.geom_type!r}; "
+                f"only Polygon and MultiPolygon are read. The resulting grid "
+                f"will not contain a face for this geometry.",
+                stacklevel=2,
+            )
 
     # Convert lists to numpy arrays at the end
     node_lon = np.array(node_lon_list)

uxarray/io/utils.py

@@ -146,6 +146,12 @@ def _is_structured(dataset: xr.Dataset, tol: float = 1e-5) -> bool:
     bool
         True if the dataset is structured with regularly spaced latitude and longitude,
         False otherwise.
+
+    Note
+    ----
+    ``_parse_grid_type`` calls this speculatively for every dataset, so a
+    negative result is the normal case for all other grid formats. It must stay
+    quiet rather than reporting why the dataset is not structured.
     """
     # Extract all 'standard_name' attributes in lower case
     standard_names = [
@@ -176,7 +182,6 @@ def _is_structured(dataset: xr.Dataset, tol: float = 1e-5) -> bool:
 
     # Ensure that latitude and longitude are one-dimensional
     if lat.ndim != 1 or lon.ndim != 1:
-        print("Latitude and/or longitude coordinates are not one-dimensional.")
         return False, None, None
 
     # Calculate the differences between consecutive latitude and longitude values
@@ -187,11 +192,6 @@ def _is_structured(dataset: xr.Dataset, tol: float = 1e-5) -> bool:
     lat_regular = np.all(np.abs(lat_diffs - lat_diffs[0]) <= tol)
     lon_regular = np.all(np.abs(lon_diffs - lon_diffs[0]) <= tol)
 
-    if not lat_regular:
-        print("Latitude coordinates are not regularly spaced.")
-    if not lon_regular:
-        print("Longitude coordinates are not regularly spaced.")
-
     return lat_regular and lon_regular, lon_name, lat_name