Feat/network cost operator

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
 # ---------------------------------------------------------------------------

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()
+