diff --git a/disscube/operators/proximity.py b/disscube/operators/proximity.py old mode 100644 new mode 100755 index 1f9cc2a..ef8bc8b --- a/disscube/operators/proximity.py +++ b/disscube/operators/proximity.py @@ -9,6 +9,7 @@ import geopandas as gpd import numpy as np +import pandas as pd import rasterio.features import xarray as xr from rasterio.warp import Resampling @@ -147,6 +148,251 @@ def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray: ) + + +class NetworkCostOperator(Operator): + """ + Least-cost network travel distance/cost from each cell centre to target + destinations via a linear transport network (e.g. roads, railways), + accounting for off-road access and on-road impedance. + + Replicates and improves upon TerraME's GPM (Generalized Proximity Matrix) + Network connectivity algorithm (used in LuccME-BR for e_connmkt and e_connport). + + Parameters in var.params: + ------------------------- + targets : str or list + Path to vector file (e.g. shapefile/GeoJSON) of target destinations + (points or polygons), or a list of [x, y] coordinates. + cost_column : str, optional + Name of attribute column with individual road impedance/cost multiplier + (e.g. 'custo_ajus'). If present, overrides status_column. + status_column : str, optional + Name of attribute column distinguishing road status (e.g. 'status', 'paved'). + inside_paved : float, optional + Cost/impedance factor on paved lines (default 1.0). + inside_unpaved : float, optional + Cost/impedance factor on unpaved lines (default 1.0). + outside : float, optional + Cost/impedance factor outside the network (default 1.0). + unit_scale : float, optional + Multiplier applied to distances (e.g. 1e-3 to convert metres to kilometres, default 1.0). + entrance : str, optional + How cell connects to the network: 'segment' (orthogonal projection, default) + or 'vertex' (nearest line vertex, replicating TerraME). + crs : str, optional + Projected CRS for metric distance calculation (e.g. 'EPSG:5880'). + """ + + name = "network_cost" + _resampling = Resampling.nearest + clip_to_grid = False + params: ClassVar[dict[str, str]] = { + "targets": "path to destinations vector layer or coordinate list [[x, y], ...]", + "cost_column": "column in network vector with per-segment cost multiplier (e.g. 'custo_ajus')", + "status_column": "column in network vector distinguishing road status (e.g. 'status')", + "inside_paved": "impedance/friction multiplier on paved lines (default: 1.0)", + "inside_unpaved": "impedance/friction multiplier on unpaved lines (default: 1.0)", + "outside": "impedance/friction multiplier off-road from cell to network (default: 1.0)", + "unit_scale": "scale factor for output units, e.g. 1e-3 for km (default: 1.0)", + "entrance": "connection rule: 'segment' (nearest edge) or 'vertex' (nearest vertex, TerraME)", + "crs": "projected CRS to measure distances in metres (e.g. 'EPSG:5880')", + } + + 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 + + crs = var.params.get("crs") + if crs is not None: + data = data.to_crs(crs) + + geoms = [g for g in data.geometry if g is not None and not g.is_empty] + if not geoms: + dist = np.full((grid.rows, grid.cols), np.nan) + return xr.DataArray(dist, dims=("y", "x"), coords={"y": grid.ys, "x": grid.xs}) + + outside_factor = float(var.params.get("outside", 1.0)) + inside_paved_factor = float(var.params.get("inside_paved", 1.0)) + inside_unpaved_factor = float(var.params.get("inside_unpaved", 1.0)) + unit_scale = float(var.params.get("unit_scale", 1.0)) + entrance_mode = var.params.get("entrance", "segment").lower() + cost_col = var.params.get("cost_column") + status_col = var.params.get("status_column") + + # 1. Carregar os destinos (alvos) + targets_param = var.params.get("targets") + target_points = [] + if isinstance(targets_param, str): + import os + import pooch + target_path = None + if os.path.exists(targets_param): + target_path = targets_param + elif (pooch.os_cache("disscube") / "raw" / targets_param).exists(): + target_path = str(pooch.os_cache("disscube") / "raw" / targets_param) + else: + matches = list((pooch.os_cache("disscube") / "raw").glob(f"**/{os.path.basename(targets_param)}")) + if matches: + target_path = str(matches[0]) + if target_path and os.path.exists(target_path): + tgdf = gpd.read_file(target_path) + if crs is not None and tgdf.crs is not None: + tgdf = tgdf.to_crs(crs) + for g in tgdf.geometry: + if g is not None and not g.is_empty: + target_points.append(g if g.geom_type == "Point" else g.centroid) + elif isinstance(targets_param, (list, tuple)): + for pt in targets_param: + if isinstance(pt, (list, tuple)) and len(pt) >= 2: + target_points.append(Point(pt[0], pt[1])) + elif isinstance(pt, Point): + target_points.append(pt) + + if not target_points: + warnings.warn( + f"network_cost: no valid targets provided for '{var.name}'; returning NaN", + RuntimeWarning, + stacklevel=2, + ) + dist = np.full((grid.rows, grid.cols), np.nan) + 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_list = [] + edges = [] + segment_geoms = [] + seg_to_edge = [] + + def get_node(pt): + key = (round(pt.x, 3), round(pt.y, 3)) + if key not in node_coords: + idx = len(node_coords) + node_coords[key] = idx + node_list.append(Point(pt.x, pt.y)) + return idx + return node_coords[key] + + for _, row in data.iterrows(): + geom = row.geometry + if geom is None or geom.is_empty: + continue + + # Determinar o fator de impedância da linha + if cost_col and cost_col in row and not pd.isna(row[cost_col]): + factor = float(row[cost_col]) + elif status_col and status_col in row: + val = str(row[status_col]).lower() + if any(term in val for term in ("unpaved", "terra", "implantada", "natural", "leito")): + factor = inside_unpaved_factor + else: + factor = inside_paved_factor + else: + factor = inside_paved_factor + + lines = geom.geoms if geom.geom_type == "MultiLineString" else [geom] + for line in lines: + coords = list(line.coords) + for i in range(len(coords) - 1): + p1 = Point(coords[i]) + p2 = Point(coords[i + 1]) + seg_dist = p1.distance(p2) + if seg_dist == 0: + continue + u = get_node(p1) + v = get_node(p2) + seg_cost = seg_dist * factor * unit_scale + edges.append((u, v, seg_cost)) + segment_geoms.append(LineString([coords[i], coords[i + 1]])) + seg_to_edge.append((u, v, factor)) + + num_nodes = len(node_list) + if num_nodes == 0 or not segment_geoms: + dist = np.full((grid.rows, grid.cols), np.nan) + return xr.DataArray(dist, dims=("y", "x"), coords={"y": grid.ys, "x": grid.xs}) + + row_idx, col_idx, cost_data = [], [], [] + for u, v, cost in edges: + row_idx.extend([u, v]) + col_idx.extend([v, u]) + cost_data.extend([cost, cost]) + + adj_matrix = csr_matrix((cost_data, (row_idx, col_idx)), shape=(num_nodes, num_nodes)) + node_tree = shapely.STRtree(node_list) + + # 3. Mapear alvos na rede com custo de acesso inicial (TerraME GPM) + target_indices = [] + target_access_penalties = [] + for t in target_points: + nearest_node_idx = int(np.atleast_1d(node_tree.query_nearest(t))[0]) + nearest_pt = node_list[nearest_node_idx] + target_indices.append(nearest_node_idx) + # Adiciona o custo de acesso euclidiano fora da rede do porto até a linha + access_dist = t.distance(nearest_pt) + target_access_penalties.append(access_dist * outside_factor * unit_scale) + + # 4. Resolver caminhos mínimos via Dijkstra + dist_matrix = dijkstra(adj_matrix, indices=target_indices, directed=False) + + # Somar a penalidade de acesso do alvo ao vetor de distâncias do grafo + if len(target_indices) == 1: + node_min_cost = dist_matrix[0] if dist_matrix.ndim == 2 else dist_matrix + node_min_cost = node_min_cost + target_access_penalties[0] + else: + penalties = np.array(target_access_penalties)[:, np.newaxis] + node_min_cost = np.min(dist_matrix + penalties, axis=0) + + # 5. Interpolar para as células da grade + xx, yy = np.meshgrid(grid.xs, grid.ys) + if crs is not None: + from pyproj import Transformer + xx, yy = Transformer.from_crs(grid.crs, crs, always_xy=True).transform(xx, yy) + + centres = shapely.points(xx.ravel(), yy.ravel()) + + if entrance_mode == "vertex": + # Regra do vértice mais próximo (reproduz a lógica discreta do TerraME) + indices, offroad_dists = node_tree.query_nearest(centres, return_distance=True, all_matches=False) + c_idx, v_idx = indices[0], indices[1] + flat_costs = offroad_dists * (outside_factor * unit_scale) + node_min_cost[v_idx] + else: + # Regra do segmento contínuo (ortogonal, mais precisa) + seg_tree = shapely.STRtree(segment_geoms) + indices, offroad_dists = seg_tree.query_nearest(centres, return_distance=True, all_matches=False) + pt_idx, seg_idx = indices[0], indices[1] + + flat_costs = np.full(len(centres), np.nan) + for i in range(len(pt_idx)): + s_idx = seg_idx[i] + c_idx = pt_idx[i] + u, v, factor = seg_to_edge[s_idx] + seg_geom = segment_geoms[s_idx] + c_pt = centres[c_idx] + + dist_along = seg_geom.project(c_pt) + p_proj = seg_geom.interpolate(dist_along) + d_u = p_proj.distance(node_list[u]) + d_v = p_proj.distance(node_list[v]) + net_cost = min(node_min_cost[u] + d_u * factor * unit_scale, + node_min_cost[v] + d_v * factor * unit_scale) + flat_costs[c_idx] = offroad_dists[i] * (outside_factor * unit_scale) + net_cost + + out_dist = flat_costs.reshape((grid.rows, grid.cols)) + return xr.DataArray(out_dist, dims=("y", "x"), coords={"y": grid.ys, "x": grid.xs}) + + raise TypeError(f"'network_cost' expects a vector source, got {type(data).__name__}") + +class GpmNetworkOperator(NetworkCostOperator): + """ + Alias for NetworkCostOperator implementing TerraME GPM Network connectivity. + """ + name = "gpm_network" + + # --------------------------------------------------------------------------- # Legacy compatibility shim # --------------------------------------------------------------------------- diff --git a/examples/04_network_connectivity.py b/examples/04_network_connectivity.py new file mode 100644 index 0000000..1261952 --- /dev/null +++ b/examples/04_network_connectivity.py @@ -0,0 +1,109 @@ +""" +04 — Network Connectivity (TerraME GPM Network replication): +Least-cost path travel distance to destinations via a transport network. + +Replicates the behavior of TerraME's GPM (Generalized Proximity Matrix) Network +connectivity algorithm, used in LuccME-BR to derive market and port accessibility +(e_connmkt and e_connport). + +Demonstrates: + - Defining linear transport infrastructure (roads/railways) with surface impedance. + - 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 +""" + +import sys +import tempfile +from pathlib import Path + +import geopandas as gpd +import numpy as np +from shapely.geometry import LineString, Point + +from disscube import CubeClient, Derivation, GridSpec, SpatialSource + + +def run_example(): + print("=" * 70) + print(" DisSCube Example 04: Network Connectivity (GPM Network)") + print("=" * 70) + + # 1. Grid specification (50 km cells in SIRGAS 2000 Polyconic / EPSG:5880) + grid = GridSpec( + id="regional_grid", + type="local", + crs="EPSG:5880", + resolution=50000, + 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"}, + ], + crs="EPSG:5880", + ) + 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", + ) + + 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) + + +if __name__ == "__main__": + run_example() + \ No newline at end of file diff --git a/tests/test_network_cost_operator.py b/tests/test_network_cost_operator.py new file mode 100644 index 0000000..b583965 --- /dev/null +++ b/tests/test_network_cost_operator.py @@ -0,0 +1,134 @@ +""" +Tests for the NetworkCostOperator and GpmNetworkOperator (TerraME GPM connectivity). +""" + +import warnings +import geopandas as gpd +import numpy as np +import pytest +from shapely.geometry import LineString, Point + +from disscube import CubeClient, Derivation, GridSpec, SpatialSource + +CRS = "EPSG:4618" +RES = 0.0045 +BBOX = [-54.8415, -3.5865, -54.459, -3.168] +GRID = GridSpec(id="g", type="local", crs=CRS, resolution=RES, bbox=BBOX) + + +@pytest.fixture +def cube(tmp_path): + c = CubeClient(catalog=str(tmp_path / "c.db"), store=str(tmp_path / "s")) + c.register_grid(GRID) + return c + + +def _roads(cube, tmp_path, sid, lines, statuses=None): + path = tmp_path / f"{sid}.geojson" + data = {"geometry": lines} + if statuses is not None: + data["status"] = statuses + gpd.GeoDataFrame(data, crs=CRS).to_file(path, driver="GeoJSON") + cube.register_spatial_source(SpatialSource(id=sid, name=sid, format="vector", crs=CRS, + asset_url=str(path))) + + +def test_network_cost_linear_corridor(cube, tmp_path): + y_road = BBOX[1] + (BBOX[3] - BBOX[1]) / 2 + road = LineString([(BBOX[0], y_road), (BBOX[2], y_road)]) + _roads(cube, tmp_path, "roads", [road], ["paved"]) + + target = [BBOX[0], y_road] + cube.derive_declarative( + Derivation( + target="conn", + source_id="roads", + operator="network_cost", + params={ + "targets": [target], + "inside_paved": 1.0, + "outside": 2.0, + }, + ), + grid_id="g", + ) + result = cube.load("conn", grid_id="g").values + assert not np.isnan(result).any() + assert (result >= 0).all() + + +def test_gpm_network_alias(cube, tmp_path): + y_road = BBOX[1] + (BBOX[3] - BBOX[1]) / 2 + road = LineString([(BBOX[0], y_road), (BBOX[2], y_road)]) + _roads(cube, tmp_path, "roads", [road]) + + target = [BBOX[0], y_road] + cube.derive_declarative( + Derivation( + target="conn_gpm", + source_id="roads", + operator="gpm_network", + params={ + "targets": [target], + "outside": 1.5, + }, + ), + grid_id="g", + ) + result = cube.load("conn_gpm", grid_id="g").values + assert not np.isnan(result).any() + + +def test_network_cost_unpaved_impedance(cube, tmp_path): + y1 = BBOX[1] + (BBOX[3] - BBOX[1]) * 0.25 + y2 = BBOX[1] + (BBOX[3] - BBOX[1]) * 0.75 + road_paved = LineString([(BBOX[0], y1), (BBOX[2], y1)]) + road_unpaved = LineString([(BBOX[0], y2), (BBOX[2], y2)]) + + _roads(cube, tmp_path, "roads_multi", [road_paved, road_unpaved], ["paved", "unpaved"]) + + targets = [[BBOX[0], y1], [BBOX[0], y2]] + cube.derive_declarative( + Derivation( + target="conn_imp", + source_id="roads_multi", + operator="network_cost", + params={ + "targets": targets, + "status_column": "status", + "inside_paved": 0.5, + "inside_unpaved": 2.0, + "outside": 1.0, + }, + ), + grid_id="g", + ) + result = cube.load("conn_imp", grid_id="g").values + assert not np.isnan(result).any() + + +def test_network_cost_multiple_destinations(cube, tmp_path): + y_road = BBOX[1] + (BBOX[3] - BBOX[1]) / 2 + road = LineString([(BBOX[0], y_road), (BBOX[2], y_road)]) + _roads(cube, tmp_path, "roads_corridor", [road], ["paved"]) + + targets = [[BBOX[0], y_road], [BBOX[2], y_road]] + cube.derive_declarative( + Derivation( + target="conn_multi_dest", + source_id="roads_corridor", + operator="network_cost", + params={ + "targets": targets, + "inside_paved": 1.0, + "outside": 1.0, + }, + ), + grid_id="g", + ) + result = cube.load("conn_multi_dest", grid_id="g").values + half_len = (BBOX[2] - BBOX[0]) / 2 + 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