clean and speed up nwa

Author: mortewle

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "ssb-sgis"
-version = "1.3.7"
+version = "1.3.8"
 description = "GIS functions used at Statistics Norway."
 authors = ["Morten Letnes <morten.letnes@ssb.no>"]
 license = "MIT"

src/sgis/geopandas_tools/runners.py

@@ -239,9 +239,6 @@ def run(
                 left, right = results
                 return left, right
             return results
-            left = np.concatenate([x[0] for x in results])
-            right = np.concatenate([x[1] for x in results])
-            return left, right
         elif (
             (self.n_jobs or 1) > 1
             and len(arr2) / self.n_jobs > 10_000
@@ -264,9 +261,6 @@ def run(
                 left, right = results
                 return left, right
             return results
-            left = np.concatenate([x[0] for x in results])
-            right = np.concatenate([x[1] for x in results])
-            return left, right
 
         return _strtree_query(arr1, arr2, method=method, **kwargs)
 

src/sgis/networkanalysis/_od_cost_matrix.py

@@ -26,13 +26,11 @@ def _od_cost_matrix(
     # calculating all-to-all distances is much faster than looping rowwise,
     # so filtering to rowwise afterwards instead
     if rowwise:
-        rowwise_df = DataFrame(
-            {
-                "origin": origins.index,
-                "destination": destinations.index,
-            }
+        keys = pd.MultiIndex.from_arrays(
+            [origins.index, destinations.index],
+            names=["origin", "destination"],
         )
-        results = rowwise_df.merge(results, on=["origin", "destination"], how="left")
+        results = results.set_index(["origin", "destination"]).loc[keys].reset_index()
 
     results["geom_ori"] = results["origin"].map(origins.geometry)
     results["geom_des"] = results["destination"].map(destinations.geometry)

src/sgis/networkanalysis/_points.py

@@ -16,10 +16,7 @@
 
 
 class Points:
-    def __init__(
-        self,
-        points: GeoDataFrame,
-    ) -> None:
+    def __init__(self, points: GeoDataFrame) -> None:
         self.gdf = points.copy()
 
     def _make_temp_idx(self, start: int) -> None:

src/sgis/networkanalysis/closing_network_holes.py

@@ -79,6 +79,7 @@ def close_network_holes(
     gdf: GeoDataFrame,
     max_distance: int | float,
     max_angle: int,
+    *,
     hole_col: str | None = "hole",
 ) -> GeoDataFrame:
     """Fills network gaps with straigt lines.
@@ -282,11 +283,13 @@ def _close_holes_all_lines(
 ) -> GeoSeries:
     k = min(len(nodes), 50)
 
+    n_dict = nodes.set_index("wkt")["n"]
+
     # make points for the deadends and the other endpoint of the deadend lines
-    deadends_target = lines.loc[lines["n_target"] == 1].rename(
+    deadends_target = lines.loc[lines["target_wkt"].map(n_dict) == 1].rename(
         columns={"target_wkt": "wkt", "source_wkt": "wkt_other_end"}
     )
-    deadends_source = lines.loc[lines["n_source"] == 1].rename(
+    deadends_source = lines.loc[lines["source_wkt"].map(n_dict) == 1].rename(
         columns={"source_wkt": "wkt", "target_wkt": "wkt_other_end"}
     )
     deadends = pd.concat([deadends_source, deadends_target], ignore_index=True)
@@ -349,12 +352,6 @@ def get_angle_difference(angle1, angle2):
         to_idx = indices[condition]
         to_wkt = nodes.iloc[to_idx]["wkt"]
 
-        # all_angles = all_angles + [
-        #     diff
-        #     for f, diff in zip(from_wkt, angles_difference[condition], strict=True)
-        #     if f not in new_sources
-        # ]
-
         # now add the wkts to the lists of new sources and targets. If the source
         # is already added, the new wks will not be added again
         new_targets = new_targets + [

src/sgis/networkanalysis/finding_isolated_networks.py

@@ -57,9 +57,7 @@ def get_connected_components(gdf: GeoDataFrame) -> GeoDataFrame:
 
     gdf["connected"] = gdf.source.map(largest_component_dict).fillna(0)
 
-    gdf = gdf.drop(
-        ["source_wkt", "target_wkt", "source", "target", "n_source", "n_target"], axis=1
-    )
+    gdf = gdf.drop(["source_wkt", "target_wkt", "source", "target"], axis=1)
 
     return gdf
 
@@ -120,8 +118,6 @@ def get_component_size(gdf: GeoDataFrame) -> GeoDataFrame:
     gdf["component_index"] = gdf["source"].map(mapper["component_index"])
     gdf["component_size"] = gdf["source"].map(mapper["component_size"])
 
-    gdf = gdf.drop(
-        ["source_wkt", "target_wkt", "source", "target", "n_source", "n_target"], axis=1
-    )
+    gdf = gdf.drop(["source_wkt", "target_wkt", "source", "target"], axis=1)
 
     return gdf

src/sgis/networkanalysis/network.py

@@ -41,7 +41,6 @@ def __init__(self, gdf: GeoDataFrame) -> None:
         self.gdf = self._prepare_network(gdf)
 
         self._make_node_ids()
-        self._percent_bidirectional = self._check_percent_bidirectional()
 
     def _make_node_ids(self) -> None:
         """Gives the lines node ids and return lines (edges) and nodes.

src/sgis/networkanalysis/networkanalysis.py

@@ -1413,7 +1413,6 @@ def _get_edges_and_weights(
         """
         if self.rules.split_lines:
             self._split_lines()
-            # self.network._make_node_ids()
             self.origins._make_temp_idx(
                 start=max(self.network.nodes.node_id.astype(int)) + 1
             )
@@ -1428,6 +1427,7 @@ def _get_edges_and_weights(
 
         self.network.gdf["src_tgt_wt"] = self.network._create_edge_ids(edges, weights)
 
+        # add edges between origins+destinations to the network nodes
         edges_start, weights_start = self.origins._get_edges_and_weights(
             nodes=self.network.nodes,
             rules=self.rules,
@@ -1587,7 +1587,7 @@ def _graph_is_up_to_date(self) -> bool:
         for points in ["origins", "destinations"]:
             if self[points] is None:
                 continue
-            if points not in self.wkts:
+            if not hasattr(self, points) or self[points] is None:
                 return False
             if self._points_have_changed(self[points].gdf, what=points):
                 return False
@@ -1617,8 +1617,6 @@ def _update_wkts(self) -> None:
         """
         self.wkts = {}
 
-        self.wkts["network"] = self.network.gdf.geometry.to_wkt().values
-
         if not hasattr(self, "origins"):
             return
 

src/sgis/networkanalysis/nodes.py

@@ -47,16 +47,14 @@ def make_node_ids(
         gdf = make_edge_coords_cols(gdf)
         geomcol1, geomcol2, geomcol_final = "source_coords", "target_coords", "coords"
 
-    # remove identical lines in opposite directions
+    # remove identical lines in opposite directions in order to get n==1 for deadends
     gdf["meters_"] = gdf.length.astype(str)
-
     sources = gdf[[geomcol1, geomcol2, "meters_"]].rename(
         columns={geomcol1: geomcol_final, geomcol2: "temp"}
     )
     targets = gdf[[geomcol1, geomcol2, "meters_"]].rename(
         columns={geomcol2: geomcol_final, geomcol1: "temp"}
     )
-
     nodes = (
         pd.concat([sources, targets], axis=0, ignore_index=True)
         .drop_duplicates([geomcol_final, "temp", "meters_"])
@@ -66,18 +64,12 @@ def make_node_ids(
     gdf = gdf.drop("meters_", axis=1)
 
     nodes["n"] = nodes.assign(n=1).groupby(geomcol_final)["n"].transform("sum")
-
     nodes = nodes.drop_duplicates(subset=[geomcol_final]).reset_index(drop=True)
-
     nodes["node_id"] = nodes.index
     nodes["node_id"] = nodes["node_id"].astype(str)
 
     gdf = _map_node_ids_from_wkt(gdf, nodes, wkt=wkt)
 
-    n_dict = {geom: n for geom, n in zip(nodes[geomcol_final], nodes["n"], strict=True)}
-    gdf["n_source"] = gdf[geomcol1].map(n_dict)
-    gdf["n_target"] = gdf[geomcol2].map(n_dict)
-
     if wkt:
         nodes["geometry"] = gpd.GeoSeries.from_wkt(nodes[geomcol_final], crs=gdf.crs)
     else:

tests/test_network_analysis.py

@@ -273,6 +273,8 @@ def not_test_get_k_routes(nwa, p):
 
 def not_test_direction(roads_oslo):
     """Check that a route that should go in separate tunnels, goes in correct tunnels."""
+    m = 5
+
     vippetangen = sg.to_gdf([10.741527, 59.9040595], crs=4326).to_crs(roads_oslo.crs)
     ryen = sg.to_gdf([10.8047522, 59.8949826], crs=4326).to_crs(roads_oslo.crs)
 
@@ -283,25 +285,21 @@ def not_test_direction(roads_oslo):
     tunnel_tofrom = sg.to_gdf([10.7724645, 59.899908], crs=4326).to_crs(roads_oslo.crs)
 
     clipped = sg.clean_clip(roads_oslo, tunnel_fromto.buffer(2000))
+
     connected_roads = sg.get_connected_components(clipped).query("connected == 1")
     directed_roads = sg.make_directed_network_norway(connected_roads, dropnegative=True)
     rules = sg.NetworkAnalysisRules(directed=True, weight="minutes")
     nwa = sg.NetworkAnalysis(directed_roads, rules=rules)
-
     route_fromto = nwa.get_route(vippetangen, ryen)
     route_tofrom = nwa.get_route(ryen, vippetangen)
-
-    m = 5
     should_be_within = route_fromto.sjoin_nearest(tunnel_fromto, distance_col="dist")
     assert should_be_within["dist"].max() < m, should_be_within["dist"]
     should_be_within = route_tofrom.sjoin_nearest(tunnel_tofrom, distance_col="dist")
     assert should_be_within["dist"].max() < m, should_be_within["dist"]
-
     should_not_be_within = route_fromto.sjoin_nearest(
         tunnel_tofrom, distance_col="dist"
     )
     assert should_not_be_within["dist"].max() > m, should_not_be_within["dist"]
-
     should_not_be_within = route_tofrom.sjoin_nearest(
         tunnel_fromto, distance_col="dist"
     )