1from warnings import warn 2 3from libpysal import cg 4from libpysal.common import requires 5import numpy 6from rtree import Rtree 7 8 9try: 10 import geopandas 11 from shapely.geometry import Point, LineString 12except ImportError: 13 msg = "geopandas/shapely not available. Some functionality will be disabled." 14 warn(msg) 15 16 17def compute_length(v0, v1): 18 """Compute the euclidean distance between two points. 19 20 Parameters 21 ---------- 22 v0 : tuple 23 Coordinate sequence in the form x,y. 24 vq : tuple 25 Coordinate sequence in the form x,y. 26 27 Returns 28 ------- 29 euc_dist : float 30 Euclidean distance. 31 32 Examples 33 -------- 34 35 >>> import spaghetti 36 >>> point1, point2 = (0,0), (1,1) 37 >>> spaghetti.util.compute_length(point1, point2) 38 1.4142135623730951 39 40 """ 41 42 euc_dist = cg.standalone.get_points_dist(v0, v1) 43 44 return euc_dist 45 46 47def get_neighbor_distances(ntw, v0, l): 48 """Get distances to the nearest vertex neighbors along 49 connecting arcs. 50 51 Parameters 52 ---------- 53 ntw : spaghetti.Network 54 A spaghetti network object. 55 v0 : int 56 The vertex ID. 57 l : dict 58 The key is a tuple (start vertex, end vertex); value is ``float``. 59 Cost per arc to travel, e.g. distance. 60 61 Returns 62 ------- 63 neighbors : dict 64 The key is an integer (vertex ID); value is ``float`` (distance). 65 66 Examples 67 -------- 68 69 >>> import spaghetti 70 >>> from libpysal import examples 71 >>> ntw = spaghetti.Network(examples.get_path("streets.shp")) 72 >>> neighs = spaghetti.util.get_neighbor_distances(ntw, 0, ntw.arc_lengths) 73 >>> numpy.round(neighs[1], 10) 74 102.6235345344 75 76 """ 77 78 # fetch links associated with vertices 79 arcs = ntw.enum_links_vertex(v0) 80 81 # create neighbor distance lookup 82 neighbors = {} 83 84 # iterate over each associated link 85 for arc in arcs: 86 87 # set distance from vertex1 to vertex2 (link length) 88 if arc[0] != v0: 89 neighbors[arc[0]] = l[arc] 90 else: 91 neighbors[arc[1]] = l[arc] 92 93 return neighbors 94 95 96def generatetree(pred): 97 """Rebuild the shortest path from root origin to destination. 98 99 Parameters 100 ---------- 101 pred : list 102 List of preceding vertices for traversal route. 103 104 Returns 105 ------- 106 tree : dict 107 The key is the root origin; value is the root origin to destination. 108 109 Examples 110 -------- 111 112 >>> import spaghetti 113 >>> from libpysal import examples 114 >>> ntw = spaghetti.Network(examples.get_path("streets.shp")) 115 >>> distance, pred = spaghetti.util.dijkstra(ntw, 0) 116 >>> tree = spaghetti.util.generatetree(pred) 117 >>> tree[3] 118 [23, 22, 20, 19, 170, 2, 0] 119 120 """ 121 122 # instantiate tree lookup 123 tree = {} 124 125 # iterate over the list of predecessor vertices 126 for i, p in enumerate(pred): 127 128 # if the route begins/ends with itself set the 129 # root vertex and continue to next iteration 130 if p == -1: 131 132 # tree keyed by root vertex with root vertex as path 133 tree[i] = [i] 134 continue 135 136 # set the initial vertex `p` as `idx` 137 idx = p 138 # and add it as the first vertex in the path 139 path = [idx] 140 141 # iterate through the path until back to home vertex 142 while idx >= 0: 143 # set the next vertex on the path 144 next_vertex = pred[idx] 145 # and redeclare the current `idx` 146 idx = next_vertex 147 148 # add the vertex to path while not at home vertex 149 if idx >= 0: 150 path.append(next_vertex) 151 152 # tree keyed by root vertex with network vertices as path 153 tree[i] = path 154 155 return tree 156 157 158def dijkstra(ntw, v0, initial_dist=numpy.inf): 159 """Compute the shortest path between a start vertex and 160 all other vertices in an origin-destination matrix. 161 162 Parameters 163 ---------- 164 ntw : spaghetti.Network 165 A spaghetti network object. 166 v0 : int 167 Start vertex ID. 168 initial_dist : float 169 Integer break point to stop iteration and return n neighbors. 170 Default is ``numpy.inf``. 171 172 Returns 173 ------- 174 distance : list 175 List of distances from vertex to all other vertices. 176 pred : list 177 List of preceeding vertices for traversal route. 178 179 Notes 180 ----- 181 182 Based on :cite:`Dijkstra1959a`. 183 184 Examples 185 -------- 186 187 >>> import spaghetti 188 >>> from libpysal import examples 189 >>> ntw = spaghetti.Network(examples.get_path("streets.shp")) 190 >>> distance, pred = spaghetti.util.dijkstra(ntw, 0) 191 >>> round(distance[196], 4) 192 5505.6682 193 >>> pred[196] 194 133 195 196 """ 197 198 # cost per arc to travel, e.g. distance 199 cost = ntw.arc_lengths 200 201 # initialize travel costs as `inf` for all distances 202 distance = [initial_dist for x in ntw.vertex_list] 203 204 # label distance to self as 0 205 distance[ntw.vertex_list.index(v0)] = 0 206 207 # instantiate set of unvisited vertices 208 unvisited = set([v0]) 209 210 # initially label as predecessor vertices with -1 as path 211 pred = [-1 for x in ntw.vertex_list] 212 213 # iterate over `unvisited` until all vertices have been visited 214 while len(unvisited) > 0: 215 216 # get vertex with the lowest value from distance 217 dist = initial_dist 218 219 for vertex in unvisited: 220 if distance[vertex] < dist: 221 dist = distance[vertex] 222 current = vertex 223 224 # remove that vertex from the set 225 unvisited.remove(current) 226 227 # get the neighbors (and costs) to the current vertex 228 neighbors = get_neighbor_distances(ntw, current, cost) 229 230 # iterate over neighbors to find least cost along path 231 for v1, indiv_cost in neighbors.items(): 232 233 # if the labeled cost is greater than 234 # the currently calculated cost 235 if distance[v1] > distance[current] + indiv_cost: 236 237 # relabel to the currently calculated cost 238 distance[v1] = distance[current] + indiv_cost 239 240 # set the current vertex as a predecessor on the path 241 pred[v1] = current 242 243 # add the neighbor vertex to `unvisted` 244 unvisited.add(v1) 245 246 # cast preceding vertices list as an array of integers 247 pred = numpy.array(pred, dtype=numpy.int) 248 249 return distance, pred 250 251 252def dijkstra_mp(ntw_vertex): 253 """Compute the shortest path between a start vertex and all other 254 vertices in the matrix utilizing multiple cores upon request. 255 256 Parameters 257 ---------- 258 ntw_vertex : tuple 259 Tuple of arguments to pass into ``dijkstra()`` as 260 (1) ``ntw`` - ``spaghetti.Network object``; 261 (2) ``vertex`` - int (start node ID) 262 263 Returns 264 ------- 265 distance : list 266 List of distances from vertex to all other vertices. 267 pred : list 268 List of preceeding vertices for traversal route. 269 270 Notes 271 ----- 272 273 Based on :cite:`Dijkstra1959a`. 274 275 Examples 276 -------- 277 278 >>> import spaghetti 279 >>> from libpysal import examples 280 >>> ntw = spaghetti.Network(examples.get_path("streets.shp")) 281 >>> distance, pred = spaghetti.util.dijkstra_mp((ntw, 0)) 282 >>> round(distance[196], 4) 283 5505.6682 284 >>> pred[196] 285 133 286 287 """ 288 289 # unpack network object and source vertex 290 ntw, vertex = ntw_vertex 291 292 # calculate shortest path distances and predecessor vertices 293 distance, pred = dijkstra(ntw, vertex) 294 295 return distance, pred 296 297 298def squared_distance_point_link(point, link): 299 """Find the squared distance between a point and a link. 300 301 Parameters 302 ---------- 303 point : tuple 304 Point coordinates (x,y). 305 link : list 306 List of 2 point coordinate tuples [(x0, y0), (x1, y1)]. 307 308 Returns 309 ------- 310 sqd : float 311 The distance squared between the point and edge. 312 nearp : numpy.ndarray 313 An array of (xb, yb); the nearest point on the edge. 314 315 Examples 316 -------- 317 318 >>> import spaghetti 319 >>> point, link = (1,1), ((0,0), (2,0)) 320 >>> spaghetti.util.squared_distance_point_link(point, link) 321 (1.0, array([1., 0.])) 322 323 """ 324 325 # cast vertices comprising the network link as an array 326 p0, p1 = [numpy.array(p) for p in link] 327 328 # cast the observation point as an array 329 p = numpy.array(point) 330 331 # subtract point 0 coords from point 1 332 v = p1 - p0 333 # subtract point 0 coords from the observation coords 334 w = p - p0 335 336 # if the point 0 vertex is the closest point along the link 337 c1 = numpy.dot(w, v) 338 if c1 <= 0.0: 339 sqd = numpy.dot(w.T, w) 340 nearp = p0 341 342 return sqd, nearp 343 344 # if the point 1 vertex is the closest point along the link 345 c2 = numpy.dot(v, v) 346 if c2 <= c1: 347 dp1 = p - p1 348 sqd = numpy.dot(dp1.T, dp1) 349 nearp = p1 350 351 return sqd, nearp 352 353 # otherwise the closest point along the link lies between p0 and p1 354 b = c1 / c2 355 bv = numpy.dot(b, v) 356 pb = p0 + bv 357 d2 = p - pb 358 sqd = numpy.dot(d2, d2) 359 nearp = pb 360 361 return sqd, nearp 362 363 364def snap_points_to_links(points, links): 365 """Place points onto closest link in a set of links (arc/edges). 366 367 Parameters 368 ---------- 369 points : dict 370 Point ID as key and (x,y) coordinate as value. 371 links : list 372 Elements are of type ``libpysal.cg.shapes.Chain`` 373 ** Note ** each element is a link represented as a chain with 374 *one head and one tail vertex* in other words one link only. 375 376 Returns 377 ------- 378 point2link : dict 379 Key [point ID (see points in arguments)]; value [a 2-tuple 380 ((head, tail), point) where (head, tail) is the target link, 381 and point is the snapped location on the link. 382 383 Examples 384 -------- 385 386 >>> import spaghetti 387 >>> from libpysal.cg.shapes import Point, Chain 388 >>> points = {0: Point((1,1))} 389 >>> link = [Chain([Point((0,0)), Point((2,0))])] 390 >>> spaghetti.util.snap_points_to_links(points, link) 391 {0: ([(0.0, 0.0), (2.0, 0.0)], array([1., 0.]))} 392 393 """ 394 395 # instantiate an rtree 396 rtree = Rtree() 397 # set the smallest possible float epsilon on machine 398 SMALL = numpy.finfo(float).eps 399 400 # initialize network vertex to link lookup 401 vertex_2_link = {} 402 403 # iterate over network links 404 for i, link in enumerate(links): 405 406 # extract network link (x,y) vertex coordinates 407 head, tail = link.vertices 408 x0, y0 = head 409 x1, y1 = tail 410 411 if (x0, y0) not in vertex_2_link: 412 vertex_2_link[(x0, y0)] = [] 413 414 if (x1, y1) not in vertex_2_link: 415 vertex_2_link[(x1, y1)] = [] 416 417 vertex_2_link[(x0, y0)].append(link) 418 vertex_2_link[(x1, y1)].append(link) 419 420 # minimally increase the bounding box exterior 421 bx0, by0, bx1, by1 = link.bounding_box 422 bx0 -= SMALL 423 by0 -= SMALL 424 bx1 += SMALL 425 by1 += SMALL 426 427 # insert the network link and its associated 428 # rectangle into the rtree 429 rtree.insert(i, (bx0, by0, bx1, by1), obj=link) 430 431 # build a KDtree on link vertices 432 kdtree = cg.KDTree(list(vertex_2_link.keys())) 433 434 point2link = {} 435 436 for pt_idx, point in points.items(): 437 438 # first, find nearest neighbor link vertices for the point 439 dmin, vertex = kdtree.query(point, k=1) 440 vertex = tuple(kdtree.data[vertex]) 441 closest = vertex_2_link[vertex][0].vertices 442 443 # Use this link as the candidate closest link: closest 444 # Use the distance as the distance to beat: dmin 445 point2link[pt_idx] = (closest, numpy.array(vertex)) 446 x0 = point[0] - dmin 447 y0 = point[1] - dmin 448 x1 = point[0] + dmin 449 y1 = point[1] + dmin 450 451 # Find all links with bounding boxes that intersect 452 # a query rectangle centered on the point with sides 453 # of length dmin * dmin 454 rtree_lookup = rtree.intersection([x0, y0, x1, y1], objects=True) 455 candidates = [cand.object.vertices for cand in rtree_lookup] 456 457 # Sorting the candidate ensures reproducible results from OS to OS. 458 # See: 459 # https://github.com/pysal/spaghetti/pull/595 460 # https://github.com/pysal/spaghetti/issues/598 461 # https://github.com/pysal/spaghetti/pull/599 462 candidates.sort(reverse=True) 463 dmin += SMALL 464 dmin2 = dmin * dmin 465 466 # of the candidate arcs, find the nearest to the query point 467 for candidate in candidates: 468 dist2cand, nearp = squared_distance_point_link(point, candidate) 469 if dist2cand <= dmin2: 470 closest = candidate 471 dmin2 = dist2cand 472 point2link[pt_idx] = (closest, nearp) 473 474 return point2link 475 476 477def network_has_cycle(adjacency): 478 """Searches for a cycle in the complete network/graph. 479 480 Parameters 481 ---------- 482 adjacency : spaghetti.Network.adjacencylist 483 Vertex adjacency relationships. 484 485 Returns 486 ------- 487 network_cycle_found : bool 488 ``True`` for a cycle being found in the network/graph, 489 otherwise ``False``. 490 491 """ 492 493 def tree_has_cycle(_parent, _v): 494 """Searches for a cycle in the subtree. 495 496 Parameters 497 ---------- 498 _parent : int 499 Root vertex for the subnetwork/graph. 500 _v : int 501 Current vertex index of in the complete network. 502 503 Returns 504 ------- 505 subtree_cycle_found : bool 506 Current recursion found a cycle in the subtree. 507 508 """ 509 510 # Set current cycle tag as False 511 subtree_cycle_found = False 512 513 # Label the current network vertex as seen 514 seen[_v] = True 515 516 # Perform recursion for all adjacent network/graph vertices 517 for rv in adjacency[_v]: 518 519 # If vertex already seen, skip it 520 if not seen[rv]: 521 # Perform recursion down the depth-first search tree 522 if tree_has_cycle(_v, rv): 523 subtree_cycle_found = True 524 break 525 526 # If an adjacent vertex has not been seen and it is not the 527 # parent of current vertex, then a cycle is present 528 elif _parent != rv: 529 subtree_cycle_found = True 530 break 531 532 return subtree_cycle_found 533 534 # set initial cycle tag as False 535 network_cycle_found = False 536 537 # Label all network/graph vertices as not seen 538 vids = list(adjacency.keys()) 539 seen = {vid: False for vid in vids} 540 541 # Perform depth-first search recursion to isolate cycles 542 for idx, v in enumerate(vids): 543 # If vertex already seen, skip it 544 if not seen[v]: 545 # Perform recursion down the depth-first search tree 546 if tree_has_cycle(-1, v): 547 network_cycle_found = True 548 break 549 550 return network_cycle_found 551 552 553def chain_constr(vcoords, arcs): 554 """Create the spatial representation of a network arc. 555 556 Parameters 557 ---------- 558 vcoords : dict 559 Vertex to coordinate lookup (see ``spaghetti.Network.vertex_coords``). 560 arcs : list 561 Arcs represented as start and end vertices. 562 563 Returns 564 ------- 565 spatial_reps : list 566 Spatial representations of arcs - ``libpysal.cg.Chain`` objects. 567 568 """ 569 spatial_reps = [_chain_constr(vcoords, vs) for vs in arcs] 570 return spatial_reps 571 572 573def _chain_constr(_vcoords, _vs): 574 """Construct a libpysal.cg.Chain object. 575 576 Parameters 577 ---------- 578 _vcoords : {dict, None} 579 See ``vcoords`` in ``get_chains()``. 580 _vs : tuple 581 Start and end vertex IDs of arc. 582 583 Returns 584 ------- 585 chain : libpysal.cg.Chain 586 Spatial representation of the arc. 587 588 """ 589 590 if _vcoords: 591 chain_vtx_points = [cg.Point((_vcoords[v])) for v in _vs] 592 else: 593 chain_vtx_points = _vs 594 chain = cg.Chain(chain_vtx_points) 595 596 return chain 597 598 599def build_chains(space_h, space_v, exterior, bounds, h=True): 600 """Generate line segments for a lattice. 601 602 Parameters 603 ---------- 604 space_h : list 605 Horizontal spacing. 606 space_v : list 607 Vertical spacing. 608 exterior : bool 609 Flag for including the outer bounding box segments. 610 bounds : list 611 Area bounds in the form - <minx,miny,maxx,maxy>. 612 h : bool 613 Generate horizontal line segments. 614 Default is ``True``. ``False`` generates vertical segments. 615 616 Returns 617 ------- 618 chains : list 619 All horizontal or vertical line segments in the lattice. 620 621 """ 622 623 # Initialize starting and ending indices 624 start_h, end_h, start_v, end_v = 0, len(space_h), 0, len(space_v) 625 626 # set inital index track back to 0 627 minus_y, minus_x = 0, 0 628 629 if h: # start track back at 1 for horizontal lines 630 minus_x = 1 631 if not exterior: # do not include borders 632 start_v += 1 633 end_v -= 1 634 635 else: # start track back at 1 for vertical lines 636 minus_y = 1 637 if not exterior: # do not include borders 638 start_h += 1 639 end_h -= 1 640 641 # Create empty line list and fill 642 chains = [] 643 644 # for element in the horizontal index 645 for plus_h in range(start_h, end_h): 646 647 # for element in the vertical index 648 for plus_v in range(start_v, end_v): 649 650 # ignore if a -1 index 651 if plus_h - minus_x == -1 or plus_v - minus_y == -1: 652 continue 653 else: 654 # Point 1 (start point + previous slot in 655 # horizontal or vertical space index) 656 p1x = bounds[0] + space_h[plus_h - minus_x] 657 p1y = bounds[1] + space_v[plus_v - minus_y] 658 p1 = cg.Point((p1x, p1y)) 659 660 # Point 2 (start point + current slot in 661 # horizontal or vertical space index) 662 p2x = bounds[0] + space_h[plus_h] 663 p2y = bounds[1] + space_v[plus_v] 664 p2 = cg.Point((p2x, p2y)) 665 666 # libpysal.cg.Chain 667 chains.append(_chain_constr(None, [p1, p2])) 668 669 return chains 670 671 672@requires("geopandas", "shapely") 673def _points_as_gdf( 674 net, vertices, vertices_for_arcs, pp_name, snapped, id_col=None, geom_col=None 675): 676 """Internal function for returning a point ``geopandas.GeoDataFrame`` 677 called from within ``spaghetti.element_as_gdf()``. 678 679 Parameters 680 ---------- 681 vertices_for_arcs : bool 682 Flag for points being an object returned (``False``) or for merely 683 creating network arcs (``True``). Set from within the parent 684 function (``spaghetti.element_as_gdf()``). 685 686 Raises 687 ------ 688 689 KeyError 690 In order to extract a ``network.PointPattern`` it must already 691 be a part of the network object. This exception is raised 692 when a ``network.PointPattern`` is being extracted that does not 693 exist within the network object. 694 695 Returns 696 ------- 697 points : geopandas.GeoDataFrame 698 Network point elements (either vertices or ``network.PointPattern`` 699 points) as a simple ``geopandas.GeoDataFrame`` of 700 ``shapely.geometry.Point`` objects with an ``"id"`` column and 701 ``"geometry"`` column. 702 703 Notes 704 ----- 705 706 1. See ``spaghetti.element_as_gdf()`` for description of arguments. 707 2. This function requires ``geopandas``. 708 709 """ 710 711 # vertices / nodes 712 if vertices or vertices_for_arcs: 713 pts_dict = net.vertex_coords 714 715 if pp_name: 716 try: 717 pp = net.pointpatterns[pp_name] 718 except KeyError: 719 err_msg = "Available point patterns are {}" 720 raise KeyError(err_msg.format(list(net.pointpatterns.keys()))) 721 722 # raw point pattern 723 if not snapped: 724 pp_pts = pp.points 725 n_pp_pts = range(len(pp_pts)) 726 pts_dict = {point: pp_pts[point]["coordinates"] for point in n_pp_pts} 727 728 # snapped point pattern 729 else: 730 pts_dict = pp.snapped_coordinates 731 732 # instantiate geopandas.GeoDataFrame 733 pts_list = list(pts_dict.items()) 734 points = geopandas.GeoDataFrame(pts_list, columns=[id_col, geom_col]) 735 points.geometry = points.geometry.apply(lambda p: Point(p)) 736 737 # additional columns 738 if not pp_name: 739 ncv_tag = "network_component_vertices" 740 if hasattr(net, ncv_tag): 741 ncv = getattr(net, ncv_tag) 742 ncv_map = {v: k for k, verts in ncv.items() for v in verts} 743 points["comp_label"] = points[id_col].map(ncv_map) 744 if pp_name: 745 c2o_tag = "component_to_obs" 746 if hasattr(pp, c2o_tag): 747 c2o = getattr(pp, c2o_tag) 748 o2c_map = {o: c for c, obs in c2o.items() for o in obs} 749 points["comp_label"] = points[id_col].map(o2c_map) 750 751 return points 752 753 754@requires("geopandas", "shapely") 755def _arcs_as_gdf(net, points, id_col=None, geom_col=None): 756 """Internal function for returning an arc ``geopandas.GeoDataFrame`` 757 called from within ``spaghetti.element_as_gdf()``. 758 759 Returns 760 ------- 761 arcs : geopandas.GeoDataFrame 762 Network arc elements as a ``geopandas.GeoDataFrame`` of 763 ``shapely.geometry.LineString`` objects with an ``"id"`` 764 column and ``geometry`` column. 765 766 Notes 767 ----- 768 769 1. See ``spaghetti.element_as_gdf()`` for description of arguments. 770 2. This function requires ``geopandas``. 771 772 """ 773 774 # arcs 775 arcs = {} 776 777 # iterate over network arcs 778 for (vtx1_id, vtx2_id) in net.arcs: 779 780 # extract vertices comprising the network arc 781 vtx1 = points.loc[(points[id_col] == vtx1_id), geom_col].squeeze() 782 vtx2 = points.loc[(points[id_col] == vtx2_id), geom_col].squeeze() 783 # create a LineString for the network arc 784 arcs[(vtx1_id, vtx2_id)] = LineString((vtx1, vtx2)) 785 786 # instantiate GeoDataFrame 787 arcs = geopandas.GeoDataFrame( 788 sorted(list(arcs.items())), columns=[id_col, geom_col] 789 ) 790 791 # additional columns 792 if hasattr(net, "network_component_labels"): 793 arcs["comp_label"] = net.network_component_labels 794 795 return arcs 796 797 798@requires("geopandas", "shapely") 799def _routes_as_gdf(paths, id_col, geom_col): 800 """Internal function for returning a shortest paths 801 ``geopandas.GeoDataFrame`` called from within 802 ``spaghetti.element_as_gdf()``. 803 804 Returns 805 ------- 806 paths : geopandas.GeoDataFrame 807 Network shortest paths as a ``geopandas.GeoDataFrame`` of 808 ``shapely.geometry.LineString`` objects with an ``"O"`` (origin), 809 ``D`` (destination), and ``geometry`` column. An additional 810 column storing the ID as a tuple is available. 811 812 Notes 813 ----- 814 815 1. See ``spaghetti.element_as_gdf()`` for description of arguments. 816 2. This function requires ``geopandas``. 817 818 """ 819 820 # isolate the origins, destinations, and geometries 821 origs = [o for (o, d), g in paths] 822 dests = [d for (o, d), g in paths] 823 geoms = [LineString(g.vertices) for (o, d), g in paths] 824 825 # instantiate as a geodataframe 826 paths = geopandas.GeoDataFrame(geometry=geoms) 827 paths["O"] = origs 828 paths["D"] = dests 829 830 if id_col: 831 paths[id_col] = paths.apply(lambda x: (x["O"], x["D"]), axis=1) 832 833 return paths 834