From a7f7f8ed811425e511bfc88f029b547a865bcd05 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 4 Oct 2026 17:39:52 +0000 Subject: [PATCH 1/2] fix(network_cost): NaN for unreachable cells, value tests, example 04 workspace, lint - network_cost reports cells that cannot reach any target as NaN (RuntimeWarning), not inf - examples/04_network_connectivity.py takes the workspace argument like the others, which fixes tests/test_examples.py - tests/test_network_cost_operator.py: value tests with hand-computed expectations (corridor, outside factor, unit_scale, cost_column, status_column, entrance segment vs vertex, MultiLineString, several targets, empty targets, targets file, raster source rejected) - ruff fixes in the touched files; proximity.py back to mode 644 --- disscube/operators/proximity.py | 19 +++- examples/04_network_connectivity.py | 130 +++++++++++----------- tests/test_network_cost_operator.py | 161 +++++++++++++++++++++++++++- 3 files changed, 245 insertions(+), 65 deletions(-) mode change 100755 => 100644 disscube/operators/proximity.py diff --git a/disscube/operators/proximity.py b/disscube/operators/proximity.py old mode 100755 new mode 100644 index ef8bc8b..0aad68c --- a/disscube/operators/proximity.py +++ b/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): @@ -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}) diff --git a/examples/04_network_connectivity.py b/examples/04_network_connectivity.py index 1261952..11dc1bb 100644 --- a/examples/04_network_connectivity.py +++ b/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)) \ No newline at end of file diff --git a/tests/test_network_cost_operator.py b/tests/test_network_cost_operator.py index b583965..bbb01dc 100644 --- a/tests/test_network_cost_operator.py +++ b/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 - \ No newline at end of file + + +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]) From 180cfd5b745c450d010be0df5c18b0ebba1abfc6 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Sun, 4 Oct 2026 15:16:34 -0300 Subject: [PATCH 2/2] fix(network_cost): annotate node_coords for mypy --- disscube/operators/proximity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/disscube/operators/proximity.py b/disscube/operators/proximity.py index 0aad68c..4a3f419 100644 --- a/disscube/operators/proximity.py +++ b/disscube/operators/proximity.py @@ -266,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 = []