fix(network_cost): NaN for unreachable cells, value tests, example 04…

disscube/operators/proximity.py

@@ -159,6 +159,9 @@ class NetworkCostOperator(Operator):
     Replicates and improves upon TerraME's GPM (Generalized Proximity Matrix)
     Network connectivity algorithm (used in LuccME-BR for e_connmkt and e_connport).
 
+    Cells that cannot reach any target through the network (for example, a road
+    component that holds no target) are NaN, with a RuntimeWarning.
+
     Parameters in var.params:
     -------------------------
     targets : str or list
@@ -202,9 +205,9 @@ class NetworkCostOperator(Operator):
     def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray:
         if isinstance(data, gpd.GeoDataFrame):
             import shapely
-            from shapely.geometry import Point, LineString
             from scipy.sparse import csr_matrix
             from scipy.sparse.csgraph import dijkstra
+            from shapely.geometry import LineString, Point
 
             crs = var.params.get("crs")
             if crs is not None:
@@ -228,6 +231,7 @@ def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray:
             target_points = []
             if isinstance(targets_param, str):
                 import os
+
                 import pooch
                 target_path = None
                 if os.path.exists(targets_param):
@@ -262,7 +266,7 @@ def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray:
                 return xr.DataArray(dist, dims=("y", "x"), coords={"y": grid.ys, "x": grid.xs})
 
             # 2. Construir o grafo da rede viária
-            node_coords = {}
+            node_coords: dict[tuple[float, float], int] = {}
             node_list = []
             edges = []
             segment_geoms = []
@@ -381,6 +385,19 @@ def get_node(pt):
                                    node_min_cost[v] + d_v * factor * unit_scale)
                     flat_costs[c_idx] = offroad_dists[i] * (outside_factor * unit_scale) + net_cost
 
+            # Cells that cannot reach any target (a road component with no target)
+            # come out of Dijkstra as inf. Report them as NaN, like every other
+            # operator does for "no value", so they do not leak into Zarr and exports.
+            unreachable = np.isinf(flat_costs)
+            if unreachable.any():
+                warnings.warn(
+                    f"network_cost: {int(unreachable.sum())} cell(s) of '{var.name}' cannot "
+                    "reach any target through the network; set to NaN",
+                    RuntimeWarning,
+                    stacklevel=2,
+                )
+                flat_costs = np.where(unreachable, np.nan, flat_costs)
+
             out_dist = flat_costs.reshape((grid.rows, grid.cols))
             return xr.DataArray(out_dist, dims=("y", "x"), coords={"y": grid.ys, "x": grid.xs})
 

examples/04_network_connectivity.py

@@ -11,8 +11,11 @@
   - Computing multi-target shortest paths with Dijkstra graph optimization.
   - Accounting for off-road Euclidean access distance from cell centres to network.
 
-Usage:
-    python examples/04_network_connectivity.py
+Everything is written to a temporary workspace (or to the directory given as
+the first argument):
+
+    python examples/04_network_connectivity.py            # temporary workspace
+    python examples/04_network_connectivity.py outputs/04 # keep the catalog and store
 """
 
 import sys
@@ -21,12 +24,12 @@
 
 import geopandas as gpd
 import numpy as np
-from shapely.geometry import LineString, Point
+from shapely.geometry import LineString
 
 from disscube import CubeClient, Derivation, GridSpec, SpatialSource
 
 
-def run_example():
+def main(workspace: Path) -> None:
     print("=" * 70)
     print(" DisSCube Example 04: Network Connectivity (GPM Network)")
     print("=" * 70)
@@ -40,70 +43,73 @@ def run_example():
         bbox=[4000000, 7000000, 4500000, 7500000],
     )
 
-    with tempfile.TemporaryDirectory() as tmp_dir:
-        workspace = Path(tmp_dir) / "workspace"
-        catalog = Path(tmp_dir) / "catalog.db"
-        cube = CubeClient(catalog=str(catalog), store=str(workspace))
-        cube.register_grid(grid)
-
-        # 2. Synthetic road network with paved trunk line and unpaved feeders
-        trunk = LineString([(4050000, 7250000), (4250000, 7250000), (4450000, 7250000)])
-        feeder_north = LineString([(4250000, 7250000), (4250000, 7450000)])
-        feeder_south = LineString([(4250000, 7250000), (4250000, 7050000)])
-
-        roads_gdf = gpd.GeoDataFrame(
-            [
-                {"geometry": trunk, "status": "paved"},
-                {"geometry": feeder_north, "status": "unpaved"},
-                {"geometry": feeder_south, "status": "unpaved"},
-            ],
+    raw = workspace / "raw"
+    raw.mkdir(parents=True, exist_ok=True)
+    cube = CubeClient(catalog=str(workspace / "catalog.db"), store=str(workspace / "store"))
+    cube.register_grid(grid)
+
+    # 2. Synthetic road network with paved trunk line and unpaved feeders
+    trunk = LineString([(4050000, 7250000), (4250000, 7250000), (4450000, 7250000)])
+    feeder_north = LineString([(4250000, 7250000), (4250000, 7450000)])
+    feeder_south = LineString([(4250000, 7250000), (4250000, 7050000)])
+
+    roads_gdf = gpd.GeoDataFrame(
+        [
+            {"geometry": trunk, "status": "paved"},
+            {"geometry": feeder_north, "status": "unpaved"},
+            {"geometry": feeder_south, "status": "unpaved"},
+        ],
+        crs="EPSG:5880",
+    )
+    roads_file = raw / "roads.geojson"
+    roads_gdf.to_file(roads_file, driver="GeoJSON")
+
+    cube.register_spatial_source(
+        SpatialSource(
+            id="roads",
+            name="Road Network",
+            format="vector",
             crs="EPSG:5880",
+            asset_url=str(roads_file),
         )
-        roads_file = Path(tmp_dir) / "roads.geojson"
-        roads_gdf.to_file(roads_file, driver="GeoJSON")
-
-        cube.register_spatial_source(
-            SpatialSource(
-                id="roads",
-                name="Road Network",
-                format="vector",
-                crs="EPSG:5880",
-                asset_url=str(roads_file),
-            )
-        )
+    )
 
-        # 3. Destination port at the eastern terminus of the highway
-        port_coords = [4450000, 7250000]
-
-        print("Deriving network connectivity with operator='network_cost'...")
-        cube.derive_declarative(
-            Derivation(
-                target="port_accessibility",
-                source_id="roads",
-                operator="network_cost",
-                params={
-                    "targets": [port_coords],
-                    "status_column": "status",
-                    "inside_paved": 1.0,
-                    "inside_unpaved": 2.5,
-                    "outside": 2.0,
-                    "crs": "EPSG:5880",
-                },
-            ),
-            grid_id="regional_grid",
-        )
+    # 3. Destination port at the eastern terminus of the highway
+    port_coords = [4450000, 7250000]
+
+    print("Deriving network connectivity with operator='network_cost'...")
+    cube.derive_declarative(
+        Derivation(
+            target="port_accessibility",
+            source_id="roads",
+            operator="network_cost",
+            params={
+                "targets": [port_coords],
+                "status_column": "status",
+                "inside_paved": 1.0,
+                "inside_unpaved": 2.5,
+                "outside": 2.0,
+                "crs": "EPSG:5880",
+            },
+        ),
+        grid_id="regional_grid",
+    )
 
-        cost_array = cube.load("port_accessibility", grid_id="regional_grid")
-        values = cost_array.values
+    cost_array = cube.load("port_accessibility", grid_id="regional_grid")
+    values = cost_array.values
 
-        print(f"Output grid dimensions: {values.shape} (rows, cols)")
-        print(f"Minimum cost (at the port): {np.nanmin(values)/1000:.1f} km")
-        print(f"Maximum cost (furthest off-road cell): {np.nanmax(values)/1000:.1f} km")
-        print(f"Mean network cost: {np.nanmean(values)/1000:.1f} km")
-        print("Network connectivity derived successfully!")
-        print("=" * 70)
+    print(f"Output grid dimensions: {values.shape} (rows, cols)")
+    print(f"Minimum cost (at the port): {np.nanmin(values)/1000:.1f} km")
+    print(f"Maximum cost (furthest off-road cell): {np.nanmax(values)/1000:.1f} km")
+    print(f"Mean network cost: {np.nanmean(values)/1000:.1f} km")
+    print("Network connectivity derived successfully!")
+    print("=" * 70)
 
 
 if __name__ == "__main__":
-    run_example()
+    if len(sys.argv) > 1:
+        main(Path(sys.argv[1]))
+    else:
+        with tempfile.TemporaryDirectory() as tmp:
+            main(Path(tmp))
 

tests/test_network_cost_operator.py

@@ -2,7 +2,6 @@
 Tests for the NetworkCostOperator and GpmNetworkOperator (TerraME GPM connectivity).
 """
 
-import warnings
 import geopandas as gpd
 import numpy as np
 import pytest
@@ -131,4 +130,162 @@ def test_network_cost_multiple_destinations(cube, tmp_path):
     row_centre = result.shape[0] // 2
     col_centre = result.shape[1] // 2
     assert result[row_centre, col_centre] <= half_len * 1.5
-
+
+
+def test_network_cost_unreachable_cells_are_nan_not_inf(cube, tmp_path):
+    """A road component without a target must give NaN (with a warning), never inf."""
+    y1 = BBOX[1] + (BBOX[3] - BBOX[1]) * 0.25
+    y2 = BBOX[1] + (BBOX[3] - BBOX[1]) * 0.75
+    reachable = LineString([(BBOX[0], y1), (BBOX[2], y1)])
+    isolated = LineString([(BBOX[0], y2), (BBOX[2], y2)])
+    _roads(cube, tmp_path, "roads_split", [reachable, isolated])
+
+    with pytest.warns(RuntimeWarning, match="cannot reach any target"):
+        cube.derive_declarative(
+            Derivation(
+                target="conn_split",
+                source_id="roads_split",
+                operator="network_cost",
+                params={"targets": [[BBOX[0], y1]]},
+            ),
+            grid_id="g",
+        )
+    result = cube.load("conn_split", grid_id="g").values
+    assert not np.isinf(result).any()
+    # grids are north-up: row 0 is nearest the isolated (northern) road, the last row the target road
+    assert np.isnan(result[0]).all()
+    assert np.isfinite(result[-1]).all()
+
+
+# ---------------------------------------------------------------------------
+# Value tests on a metric grid: 10 x 10 cells of 1 km, so every expected number
+# can be worked out by hand. Row 0 is north; the road runs along y = 5000 m, the
+# line between rows 4 and 5, so cells of rows 4 and 5 are 500 m from it.
+# ---------------------------------------------------------------------------
+
+PCRS = "EPSG:31983"
+PGRID = GridSpec(id="p", type="local", crs=PCRS, resolution=1000.0, bbox=[0.0, 0.0, 10000.0, 10000.0])
+ROAD = LineString([(0, 5000), (10000, 5000)])
+PORT = [0, 5000]
+
+
+@pytest.fixture
+def pcube(tmp_path):
+    c = CubeClient(catalog=str(tmp_path / "pc.db"), store=str(tmp_path / "ps"))
+    c.register_grid(PGRID)
+    return c
+
+
+def _metric_roads(cube, tmp_path, sid, lines, **columns):
+    path = tmp_path / f"{sid}.geojson"
+    gpd.GeoDataFrame({"geometry": lines, **columns}, crs=PCRS).to_file(path, driver="GeoJSON")
+    cube.register_spatial_source(SpatialSource(id=sid, name=sid, format="vector", crs=PCRS,
+                                               asset_url=str(path)))
+
+
+def _derive(cube, sid, **params):
+    cube.derive_declarative(
+        Derivation(target="t", source_id=sid, operator="network_cost", params=params),
+        grid_id="p",
+    )
+    return cube.load("t", grid_id="p").values
+
+
+def test_corridor_exact_values(pcube, tmp_path):
+    """Cell at column c: 500 m off the road + (c * 1000 + 500) m along it."""
+    _metric_roads(pcube, tmp_path, "r", [ROAD])
+    out = _derive(pcube, "r", targets=[PORT])
+    expected = 500.0 + (np.arange(10) * 1000.0 + 500.0)
+    np.testing.assert_allclose(out[4], expected)
+    np.testing.assert_allclose(out[5], expected)  # symmetric about the road
+
+
+def test_outside_factor_scales_only_the_off_road_leg(pcube, tmp_path):
+    _metric_roads(pcube, tmp_path, "r", [ROAD])
+    out = _derive(pcube, "r", targets=[PORT], outside=2.0)
+    np.testing.assert_allclose(out[4], 2 * 500.0 + (np.arange(10) * 1000.0 + 500.0))
+
+
+def test_unit_scale_converts_units(pcube, tmp_path):
+    _metric_roads(pcube, tmp_path, "r", [ROAD])
+    metres = _derive(pcube, "r", targets=[PORT])
+    km = _derive(pcube, "r", targets=[PORT], unit_scale=1e-3)
+    np.testing.assert_allclose(km, metres / 1000.0)
+
+
+def test_cost_column_overrides_the_status_factor(pcube, tmp_path):
+    """A per-feature multiplier doubles the on-road leg and leaves the off-road one alone."""
+    _metric_roads(pcube, tmp_path, "r", [ROAD], custo_ajus=[2.0], status=["paved"])
+    out = _derive(pcube, "r", targets=[PORT], cost_column="custo_ajus",
+                  status_column="status", inside_paved=1.0)
+    np.testing.assert_allclose(out[4], 500.0 + 2.0 * (np.arange(10) * 1000.0 + 500.0))
+
+
+def test_status_column_unpaved_factor(pcube, tmp_path):
+    _metric_roads(pcube, tmp_path, "r", [ROAD], status=["unpaved"])
+    out = _derive(pcube, "r", targets=[PORT], status_column="status",
+                  inside_paved=1.0, inside_unpaved=3.0)
+    np.testing.assert_allclose(out[4], 500.0 + 3.0 * (np.arange(10) * 1000.0 + 500.0))
+
+
+def test_entrance_segment_vs_vertex(pcube, tmp_path):
+    """The road has only two vertices, at x = 0 and x = 10000.
+
+    'segment' joins a cell at the foot of the perpendicular, 'vertex' (the TerraME rule)
+    at the nearest vertex, so a cell in the middle goes to the far end and back.
+    """
+    _metric_roads(pcube, tmp_path, "r", [ROAD])
+    seg = _derive(pcube, "r", targets=[PORT], entrance="segment")
+    vert = _derive(pcube, "r", targets=[PORT], entrance="vertex")
+    assert seg[4, 5] == pytest.approx(500.0 + 5500.0)
+    # nearest vertex of the centre (5500, 4500) is (10000, 5000): 4500 across, 500 up
+    # (off-road), then 10000 m back along the road to the port at x = 0
+    assert vert[4, 5] == pytest.approx(np.hypot(4500.0, 500.0) + 10000.0)
+    assert vert[4, 5] > seg[4, 5]  # mid-road cell: the vertex rule detours to the far end
+
+
+def test_multilinestring_is_one_connected_road(pcube, tmp_path):
+    from shapely.geometry import MultiLineString
+
+    halves = MultiLineString([[(0, 5000), (5000, 5000)], [(5000, 5000), (10000, 5000)]])
+    _metric_roads(pcube, tmp_path, "r", [halves])
+    out = _derive(pcube, "r", targets=[PORT])
+    np.testing.assert_allclose(out[4], 500.0 + (np.arange(10) * 1000.0 + 500.0))
+
+
+def test_nearest_target_wins(pcube, tmp_path):
+    _metric_roads(pcube, tmp_path, "r", [ROAD])
+    out = _derive(pcube, "r", targets=[[0, 5000], [10000, 5000]])
+    # symmetric: the two ends see the same cost, and the middle is the costliest column
+    np.testing.assert_allclose(out[4], out[4][::-1])
+    assert out[4].argmax() in (4, 5)
+
+
+def test_no_targets_gives_nan_and_warns(pcube, tmp_path):
+    _metric_roads(pcube, tmp_path, "r", [ROAD])
+    with pytest.warns(RuntimeWarning, match="no valid targets"):
+        out = _derive(pcube, "r", targets=[])
+    assert np.isnan(out).all()
+
+
+def test_targets_from_a_vector_file(pcube, tmp_path):
+    _metric_roads(pcube, tmp_path, "r", [ROAD])
+    port_file = tmp_path / "port.geojson"
+    gpd.GeoDataFrame({"geometry": [Point(*PORT)]}, crs=PCRS).to_file(port_file, driver="GeoJSON")
+    from_file = _derive(pcube, "r", targets=str(port_file))
+    from_list = _derive(pcube, "r", targets=[PORT])
+    np.testing.assert_allclose(from_file, from_list)
+
+
+def test_network_cost_rejects_a_raster_source(pcube, tmp_path):
+    import rasterio
+    from rasterio.transform import from_origin
+
+    path = tmp_path / "r.tif"
+    with rasterio.open(path, "w", driver="GTiff", height=10, width=10, count=1, dtype="float32",
+                       crs=PCRS, transform=from_origin(0, 10000, 1000, 1000)) as dst:
+        dst.write(np.ones((1, 10, 10), dtype="float32"))
+    pcube.register_spatial_source(SpatialSource(id="ras", name="ras", format="raster", crs=PCRS,
+                                                asset_url=str(path)))
+    with pytest.raises(TypeError, match="vector"):
+        _derive(pcube, "ras", targets=[PORT])