Fixes centroid connector creation

aequilibrae/project/network/connector_creation.py

@@ -1,122 +1,67 @@
-from math import pi, sqrt
 from sqlite3 import Connection
 from typing import Optional
 
-import numpy as np
-from scipy.cluster.vq import kmeans2, whiten
-from scipy.spatial.distance import cdist
-import shapely.wkb
-from shapely.geometry import LineString
+from shapely.geometry import LineString, Polygon
 from aequilibrae.utils.db_utils import commit_and_close
 
 INFINITE_CAPACITY = 99999
 
 
 def connector_creation(
-    geo, zone_id: int, srid: int, mode_id: str, network, link_types="", connectors=1, conn_: Optional[Connection] = None
+    zone_id: int,
+    mode_id: str,
+    network,
+    link_types="",
+    connectors=1,
+    conn_: Optional[Connection] = None,
+    delimiting_area: Polygon = None,
 ):
     if len(mode_id) > 1:
         raise Exception("We can only add centroid connectors for one mode at a time")
 
-    with conn_ or commit_and_close(network.project.connect()) as conn:
+    with conn_ or commit_and_close(network.project.path_to_file) as conn:
         logger = network.project.logger
-        if conn.execute("select count(*) from nodes where node_id=?", [zone_id]).fetchone() is None:
+        if sum(conn.execute("select count(*) from nodes where node_id=?", [zone_id]).fetchone()) == 0:
             logger.warning("This centroid does not exist. Please create it first")
             return
 
-        proj_nodes = network.nodes
-        node = proj_nodes.get(zone_id)
         sql = "select count(*) from links where a_node=? and instr(modes,?) > 0"
         if conn.execute(sql, [zone_id, mode_id]).fetchone()[0] > 0:
             logger.warning("Mode is already connected")
             return
 
-        if len(link_types) > 0:
-            lt = f"*[{link_types}]*"
-        else:
-            lt = "".join([x[0] for x in conn.execute("Select link_type_id from link_types").fetchall()])
-            lt = f"*[{lt}]*"
+    proj_nodes = network.nodes.data
 
-        sql = """select node_id, ST_asBinary(geometry), modes, link_types from nodes where ST_Within(geometry, GeomFromWKB(?, ?)) and
-                        (nodes.rowid in (select rowid from SpatialIndex where f_table_name = 'nodes' and
-                        search_frame = GeomFromWKB(?, ?)))
-                and link_types glob ? and instr(modes, ?)>0"""
+    centroid = proj_nodes.query("node_id == @zone_id")  # type: gpd.GeoDataFrame
+    centroid = centroid.rename(columns={"node_id": "zone_id"})[["zone_id", "geometry"]]
 
-        # We expand the area by its average radius until it is 20 times
-        # beginning with a strict search within the zone
-        buffer = 0
-        increase = sqrt(geo.area / pi)
-        dt = []
-        while dt == [] and buffer <= increase * 10:
-            wkb = geo.buffer(buffer).wkb
-            dt = conn.execute(sql, [wkb, srid, wkb, srid, lt, mode_id]).fetchall()
-            buffer += increase
+    nodes = proj_nodes.query("is_centroid != 1 and modes.str.contains(@mode_id)", engine="python")
 
-        if buffer > increase:
-            msg = f"Could not find node inside zone {zone_id}. Search area was expanded until we found a suitable node"
-            logger.warning(msg)
-        if dt == []:
-            logger.warning(
-                f"FAILED! Could not find suitable nodes to connect within 5 times the diameter of zone {zone_id}."
-            )
-            return
-
-        coords = []
-        nodes = []
-        for node_id, wkb, modes, link_types in dt:
-            geo = shapely.wkb.loads(wkb)
-            coords.append([geo.x, geo.y])
-            nodes.append(node_id)
-
-        num_connectors = connectors
-        if len(nodes) == 0:
-            raise Exception("We could not find any candidate nodes that satisfied your criteria")
-        elif len(nodes) < connectors:
-            logger.warning(
-                f"We have fewer possible nodes than required connectors for zone {zone_id}. Will connect all of them."
-            )
-            num_connectors = len(nodes)
-
-        if num_connectors == len(coords):
-            all_nodes = nodes
-        else:
-            features = np.array(coords)
-            whitened = whiten(features)
-            centroids, allocation = kmeans2(whitened, num_connectors)
+    if len(link_types) > 0:
+        nodes = nodes[nodes.link_types.str.contains("|".join(list(link_types)))]
 
-            all_nodes = set()
-            for i in range(num_connectors):
-                nds = [x for x, y in zip(nodes, list(allocation)) if y == i]
-                centr = centroids[i]
-                positions = [x for x, y in zip(whitened, allocation) if y == i]
-                if positions:
-                    dist = cdist(np.array([centr]), np.array(positions)).flatten()
-                    node_to_connect = nds[dist.argmin()]
-                    all_nodes.add(node_to_connect)
+    if delimiting_area is not None:
+        nodes = nodes[nodes.geometry.within(delimiting_area)]
 
-        nds = list(all_nodes)
-        data = [zone_id] + nds
-        sql = f'select b_node from links where a_node=? and b_node in ({",".join(["?"] * len(nds))})'
-        data = [x[0] for x in conn.execute(sql, data).fetchall()]
+    if nodes.empty:
+        zone_id = centroid["zone_id"].values[0]
+        logger.warning(f"No nodes found for centroid {zone_id} (mode {mode_id} and link types {link_types})")
+        return
 
-        if data:
-            qry = ",".join(["?"] * len(data))
-            dt = [mode_id, zone_id] + data
-            conn.execute(f"Update links set modes=modes || ? where a_node=? and b_node in ({qry})", dt)
-            nds = [x for x in nds if x not in data]
-            logger.warning(f"Mode {mode_id} added to {len(data)} existing centroid connectors for zone {zone_id}")
+    joined = nodes[["node_id", "geometry"]].sjoin_nearest(centroid, distance_col="distance_connector")
+    joined = joined.nsmallest(connectors, "distance_connector")
 
+    centr_geo = centroid.geometry.values[0]
     links = network.links
-    for node_to_connect in nds:
+    for _, rec in joined.iterrows():
         link = links.new()
-        node_to = proj_nodes.get(node_to_connect)
-        link.geometry = LineString([node.geometry, node_to.geometry])
+        link.geometry = LineString([centr_geo, rec.geometry])
         link.modes = mode_id
         link.direction = 0
         link.link_type = "centroid_connector"
         link.name = f"centroid connector zone {zone_id}"
         link.capacity_ab = INFINITE_CAPACITY
         link.capacity_ba = INFINITE_CAPACITY
         link.save()
-    if nds:
-        logger.warning(f"{len(nds)} new centroid connectors for mode {mode_id} added for centroid {zone_id}")
+    if not joined.empty:
+        logger.warning(f"{joined.shape[0]} new centroid connectors for mode {mode_id} added for centroid {zone_id}")

aequilibrae/project/network/node.py

@@ -113,7 +113,14 @@ def __save_existing_node(self):
         sql = f"Update Nodes set {txts}"
         return data, sql
 
-    def connect_mode(self, area: Polygon, mode_id: str, link_types="", connectors=1, conn: Optional[Connection] = None):
+    def connect_mode(
+        self,
+        mode_id: str,
+        link_types="",
+        connectors=1,
+        conn: Optional[Connection] = None,
+        area: Optional[Polygon] = None,
+    ):
         """Adds centroid connectors for the desired mode to the network file
 
         Centroid connectors are created by connecting the zone centroid to one or more nodes selected from
@@ -129,28 +136,28 @@ def connect_mode(self, area: Polygon, mode_id: str, link_types="", connectors=1,
         If fewer candidates than required connectors are found, all candidates are connected.
 
         :Arguments:
-            **area** (:obj:`Polygon`): Initial area where AequilibraE will look for nodes to connect
 
             **mode_id** (:obj:`str`): Mode ID we are trying to connect
 
             **link_types** (:obj:`str`, *Optional*): String with all the link type IDs that can
             be considered. eg: yCdR. Defaults to ALL link types
 
             **connectors** (:obj:`int`, *Optional*): Number of connectors to add. Defaults to 1
+
+            **area** (:obj:`Polygon`, *Optional*): Area limiting the search for connectors
         """
         if self.is_centroid != 1 or self.__original__["is_centroid"] != 1:
             self._logger.warning("Connecting a mode only makes sense for centroids and not for regular nodes")
             return
 
         connector_creation(
-            area,
-            self.node_id,
-            self.__srid__,
-            mode_id,
+            zone_id=self.node_id,
+            mode_id=mode_id,
             link_types=link_types,
             connectors=connectors,
             network=self.project.network,
             conn_=conn,
+            delimiting_area=area,
         )
 
     def __setattr__(self, instance, value) -> None:

aequilibrae/project/zone.py

@@ -90,18 +90,16 @@ def add_centroid(self, point: Point, robust=True) -> None:
             data = [self.zone_id, point.wkb, self.__srid__]
             conn.execute(sql, data)
 
-    def connect_mode(self, mode_id: str, link_types="", connectors=1, conn: Optional[Connection] = None) -> None:
+    def connect_mode(
+        self, mode_id: str, link_types="", connectors=1, conn: Optional[Connection] = None, limit_to_zone=True
+    ) -> None:
         """Adds centroid connectors for the desired mode to the network file
 
         Centroid connectors are created by connecting the zone centroid to one or more nodes selected from
         all those that satisfy the mode and link_types criteria and are inside the zone.
 
-        The selection of the nodes that will be connected is done simply by computing running the
-        `KMeans2 <https://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.vq.kmeans2.html>`_
-        clustering algorithm from SciPy and selecting the nodes closest to each cluster centroid.
-
-        When there are no node candidates inside the zone, the search area is progressively expanded until
-        at least one candidate is found.
+        The selection of the nodes that will be connected is done simply by searching for the node closest to the
+        zone centroid, or the N closest nodes to the centroid.
 
         If fewer candidates than required connectors are found, all candidates are connected.
 
@@ -112,16 +110,21 @@ def connect_mode(self, mode_id: str, link_types="", connectors=1, conn: Optional
             eg: yCdR. Defaults to ALL link types
 
             **connectors** (:obj:`int`, *Optional*): Number of connectors to add. Defaults to 1
+
+            **conn** (:obj:`sqlite3.Connection`, *Optional*): Connection to the database.
+
+            **limit_to_zone** (:obj:`bool`): Limits the search for nodes inside the zone. Defaults to ``True``.
         """
+
+        area = self.geometry if limit_to_zone else None
         connector_creation(
-            self.geometry,
             zone_id=self.zone_id,
-            srid=self.__srid__,
             mode_id=mode_id,
             link_types=link_types,
             connectors=connectors,
             network=self.project.network,
             conn_=conn,
+            delimiting_area=area,
         )
 
     def disconnect_mode(self, mode_id: str) -> None:

docs/source/examples/creating_models/plot_create_zoning.py

@@ -153,8 +153,7 @@
 # %%
 # When connecting a centroid not associated with a zone, we need to tell AequilibraE what is the initial area around
 # the centroid that needs to be considered when looking for candidate nodes.
-# Distance here is in degrees, so 0.01 is equivalent to roughly 1.1km
-airport.connect_mode(airport.geometry.buffer(0.01), mode_id="c", link_types="ytrusP", connectors=1)
+airport.connect_mode(mode_id="c", link_types="ytrusP", connectors=1)
 
 # %%
 project.close()

tests/aequilibrae/project/test_zone.py

@@ -154,9 +154,8 @@ def test_connect_mode(self):
         self.assertIsNot(0, curr.fetchone()[0], "failed to add connectors for mode t")
 
         # Cannot connect a centroid that does not exist
-        with self.assertRaises(ValueError):
-            zone2 = zones.get(2)
-            zone2.connect_mode("c")
+        zone2 = zones.get(2)
+        zone2.connect_mode("c", conn=self.proj.conn)
 
     def test_disconnect_mode(self):
         self.__change_project()