Inspire a community of professionals skilled in spatiotemporal analysis for community health and safety.
Prototype an algorithm that places distribution poles to connect service drop poles.
Roy Hyunjin Han
The next step is to connect the service drop poles in a mini grid that minimizes the length of distribution line.
20171004-2100 - 20171018-2100: 2 weeks estimated
Use brute force kruskal
Update brute force kruskal with kdtree20171004-2100 - 20171004-2200: 60 minutes
Case 1: Without roads, without unpassable obstacles
Case 2: Without roads, with unpassable obstacles
Case 3: With roads, without unpassable obstacles
Case 4: With roads, with unpassable obstacles
Given a road network, how do we isolate the part of the road network that is relevant to our path?
I think that we should assume the road network is fragmented.
But we still need to be able to identify the part of the road network that is relevant to our path.
I think what we can do to deal with both fragmented paths and path tracing is that we can first connect to the nearest road and then step wise move using a d-distance buffer. That means that we start with the intersection connection point of the node with the line, then search in a d-distance buffer for the next fragment or fork that is closer to our end point. We might have to do a couple of tracing paths to find the shorter candidate, as this approach might not end up with the shorter route.
Another option is to follow the road until it goes the wrong way, then do a bee line for the other point.
20171012-0915 - 20171012-0945: 30 minutes
We're going to start with the simplest, naive case of no obstacles, no roads and xy points.
import numpy as np
points = np.random.rand(10, 2)
We should implement some kind of deduplication. I think we can use maximum_drop_length as our deduplication radius. Mostly, if two poles are within a certain radius of each other, then consider them to be the same. For all the ones that are the same, get the centroid? Well in this specific case, the distribution line will have to originate from a pole. So we can compute the centroid and then get the drop pole that is closest to the centroid.
import numpy as np
from itertools import combinations
from scipy.spatial.distance import euclidean as get_distance
points = np.random.rand(10, 2)
group_radius = 0.25
indices = range(len(points))
group_by_index = {}
for index1, index2 in combinations(indices, 2):
point1, point2 = points[index1], points[index2]
get_distance(point1, point2)
# for each point, compare to another point
# if the point is within cluster radius, then put it in the same cluster
# find centroid of each cluster
# find point closest to centroid
20171012-2100 - 20171012-2130: 30 minutes
# Cut line into separated points (convert into dotted line)
from shapely.geometry import LineString
line = LineString([(0, 0), (1, 1), (2, 0)])
line
list(line.coords)
print(line.interpolate(0.5))
print(line.interpolate(1))
print(line.interpolate(2))
line.length
print(line.interpolate(line.length))
print(line.interpolate(line.length - 0.1))
line = LineString([(0, 0), (2, 2), (2, 0)])
line
from shapely.geometry import MultiPoint
MultiPoint(line.coords)
maximum_distribution_length = 1
maximum_distribution_length
line.coords[0:2]
LineString(line.coords[0:2]).length
LineString(line.coords[0:2]).length / maximum_distribution_length
int(LineString(line.coords[0:2]).length / maximum_distribution_length) + 1
LineString(line.coords[0:2]).length / 3.
from shapely.geometry import MultiPoint
MultiPoint([
line.coords[0],
line.interpolate(0.9428090415820635),
line.interpolate(0.9428090415820635 * 2),
line.coords[1],
])
20171016-0800 - 20171016-0930: 90 minutes (Salah, Roy)
from shapely.geometry import (
GeometryCollection, LineString, Point)
points = [Point(x) for x in [
(0, 3),
(1, 0),
(4, 7),
(7, 0),
(8, 5)]]
line = LineString([(-1, 0), (8, 3)])
GeometryCollection(points + [line])
list(line.coords)
list(line.interpolate(0).coords)
list(line.interpolate(1).coords)
list(line.interpolate(5).coords)
list(line.interpolate(line.length).coords)
l = LineString([(0, 0), (10, 10), (0, 0)])
slice_count = 4
for index in range(slice_count + 1):
print(tuple(l.interpolate(
index * l.length / slice_count).coords[0]))
GeometryCollection(points + [line])
GeometryCollection(points + [
line.interpolate(line.project(points[0])),
])
import networkx as nx
from itertools import combinations
from scipy.spatial.distance import euclidean as get_distance
g = nx.Graph()
indices = range(len(points))
for index1, index2 in combinations(indices, 2):
point1, point2 = points[index1], points[index2]
distance = get_distance(point1, point2)
g.add_edge(index1, index2, weight=distance)
g.edges(data=True)
def get_closest_point(line, point):
distance = line.project(point)
point_on_line = line.interpolate(distance)
return point_on_line
list(get_closest_point(line, points[0]).coords)
20171018-0830 - 20171018-0900: 30 minutes (Salah, Roy)
vim drop-pole-points.csv
wkt
POINT (0 3)
POINT (1 0)
POINT (4 7)
POINT (7 0)
POINT (8 5)
vim road-single.csv
wkt
"LINESTRING (-1 0, 8 3)"from pandas import DataFrame, read_csv
from shapely import wkt
def get_geometries_from_table(table):
geometry_wkts = table['wkt']
return [wkt.loads(x) for x in geometry_wkts]
def get_table_from_geometries(geometries):
return DataFrame({'wkt': [x.wkt for x in geometries]})
def load_geometry_table(source_path):
table = read_csv(source_path)
table.columns = table.columns.str.lower()
return get_geometries_from_table(table)
def save_geometry_table(target_path, geometries):
table = get_table_from_geometries(geometries)
table.to_csv(target_path, index=False)
import networkx as nx
from argparse import ArgumentParser
from itertools import combinations
from pprint import pprint
from scipy.spatial.distance import euclidean as get_distance
def get_closest_point(line, point):
distance = line.project(point)
point_on_line = line.interpolate(distance)
return point_on_line
def get_min_tree(nodes, roads, onroadcost, offroadcost):
# indices = range(len(nodes))
g = nx.Graph()
for point1, point2 in combinations(nodes, 2):
distance = get_distance(point1, point2)
g.add_edge(tuple(point1.coords), tuple(point2.coords), weight=(distance * offroadcost))
print(len(list(combinations(ds, 2))))
road_nodes = []
# roads should have discounted weight
for n in nodes:
xy = get_closest_point(roads, n)
road_nodes.append(xy)
distance = get_distance(xy, n)
g.add_edge(tuple(n.coords), tuple(xy.coords), weight=(distance * offroadcost))
print(len(g.edges(data=True)))
pprint(g.edges(data=True))
# get edges between nodes on road
road_nodes = sorted(road_nodes, key=lambda point: roads.project(point))
for i, j in zip(road_nodes[:-1], road_nodes[1:]):
g.add_edge(tuple(i.coords), tuple(j.coords), weight=(get_distance(i, j) * onroadcost))
# run path finding algorithm on graph
# nx.shortest_paths.dijkstra_path(g)
# problem, only want k-minimum spanning tree
mst = nx.minimum_spanning_tree(g)
return mst
ds = load_geometry_table('drop-pole-points.csv')
l = load_geometry_table('road-single.csv')[0]
distribution_line_on_road_cost_per_meter = 1
distribution_line_off_road_cost_per_meter = 2
mst = get_min_tree(ds, l, onroadcost=1, offroadcost=2)
pprint(mst)
+ Make example drop pole points in CSV format
+ Make example single road in CSV format
+ Write code to load from CSV
+ Write code to save to CSV
+ Turn drop pole points and single road into script20171018-2115 - 20171018-2145: 30 minutes
Working with Salah for the past week, I think we have enough material to put together a small tool.
We worked on the road problem, but ran into an issue where the road points were forcibly being included in the final network.
So actually the road problem has more complexity.
Then we should really just work on the basic version first.
We can generalize to drop pole polygons using projection, I think.
drop_pole_table_path = 'drop-pole-points.csv'
distribution_pole_interval_in_meters = 1
from pandas import read_csv
drop_pole_table = read_csv(drop_pole_table_path)
drop_pole_table
drop_pole_table.columns
from shapely import wkt
from shapely.geometry import MultiPoint
drop_pole_points = [wkt.loads(x) for x in drop_pole_table['wkt']]
MultiPoint(drop_pole_points)
import networkx as nx
from itertools import combinations
from scipy.spatial.distance import euclidean as get_distance
graph = nx.Graph()
for p1, p2 in combinations(drop_pole_points, 2):
distance = get_distance(p1, p2)
graph.add_edge(
tuple(p1.coords[0]),
tuple(p2.coords[0]), weight=distance)
graph.edges(data=True)
tree = nx.minimum_spanning_tree(graph)
tree.edges(data=True)
# Convert edges into poles
from shapely.geometry import LineString, MultiPoint, Point
distribution_pole_points = []
for p1_coords, p2_coords in tree.edges():
line_length_in_meters = get_distance(p1_coords, p2_coords)
distribution_pole_count = int(
line_length_in_meters / distribution_pole_interval_in_meters)
segment_count = distribution_pole_count + 1
adjusted_interval_in_meters = line_length_in_meters / float(
segment_count)
line = LineString([p1_coords, p2_coords])
print(p1_coords, p2_coords)
for segment_index in range(segment_count):
distribution_pole_point = line.interpolate(
adjusted_interval_in_meters * segment_index)
print(distribution_pole_point)
distribution_pole_points.append(distribution_pole_point)
print(p2_coords)
distribution_pole_points.append(Point(p2_coords))
print('***')
MultiPoint(distribution_pole_points)
+ Prototype distribution pole tool
+ Finish line splicing algorithm
20171107-2230 - 20171107-0000: 90 minutes (Salah)
from geotable import ColorfulGeometryCollection
from shapely.geometry import GeometryCollection, MultiPolygon, Polygon, Point
source_point = Point(0, 0)
target_point = Point(9, 9)
obstacles = [
Polygon([(1, 1), (1, 3), (4, 3), (4, 1), (1, 1)]),
Polygon([(0, 5), (0, 7), (2, 7), (2, 5), (0, 5)]),
Polygon([(4, 6), (4, 8), (9, 8), (9, 6), (4, 6)]),
]
ColorfulGeometryCollection([
MultiPolygon(obstacles),
source_point,
target_point,
], colors=[
'gray', 'red', 'blue',
])
MultiPolygon(obstacles).bounds
import numpy as np
np.arange(0.5, 5.5, 1)
import numpy as np
from itertools import chain, product
from shapely.geometry import Point
obstacle_bounds = MultiPolygon(obstacles).bounds
xs = [source_point.x, target_point.x] + [
obstacle_bounds[0], obstacle_bounds[2]]
ys = [source_point.y, target_point.y] + [
obstacle_bounds[1], obstacle_bounds[3]]
mesh_xys = product(
np.arange(min(xs), max(xs) + 1),
np.arange(min(ys), max(ys) + 1))
mesh_points = [Point(xy) for xy in mesh_xys]
GeometryCollection(mesh_points)
eligible_points = []
for point in mesh_points:
for obstacle in obstacles:
if obstacle.contains(point):
break
else:
eligible_points.append(point)
GeometryCollection(eligible_points)
import math
import networkx as nx
from shapely.geometry import LineString
graph = nx.Graph()
xys = [(_.x, _.y) for _ in eligible_points]
def is_valid_line(xy1, xy2):
if xy1 not in xys:
return False
if xy2 not in xys:
return False
for obstacle in obstacles:
line = LineString([xy1, xy2])
if line.intersects(obstacle):
return False
return True
for point in eligible_points:
x = point.x
y = point.y
adjacent_xys = [
(x + 1, y),
(x, y + 1),
(x - 1, y),
(x, y - 1),
]
for xy in filter(lambda _: _ in xys, adjacent_xys):
graph.add_edge((x, y), xy, weight=1)
diagonal_xys = [
(x + 1, y + 1),
(x - 1, y + 1),
(x - 1, y - 1),
(x + 1, y - 1),
]
for xy in filter(lambda _: is_valid_line(_, xy), diagonal_xys):
graph.add_edge((x, y), xy, weight=math.sqrt(2))
source_xy = source_point.x, source_point.y
target_xy = target_point.x, target_point.y
path = nx.algorithms.dijkstra_path(graph, source_xy, target_xy)
line = LineString(path)
line
ColorfulGeometryCollection([
GeometryCollection(obstacles),
line,
source_point,
target_point,
], colors=['gray', 'black', 'red', 'blue'])
20171115-1630 - 20171115-1700: 30 minutes
We should integrate the command line logs into this notebook.
+ Integrate command line logs into this notebook
_ Apply to case of latitude longitude
_ Use drop pole polygons instead of points
_ Finish point deduplication algorithm
+ Ask if we should assume that the road network is in one piece or if it is fragmented
+ Process tasksWrite function that takes node, integer k
Return a set of line segments connecting that node to other k other nodes
Write function that takes set of nodes, integer k
Returns a set of line segments that connects each node to the k nearest nodes
Write function that takes line, obstacle, buffer distance
Return modified line that goes around the obstacle
Write function that takes two nodes and a road network and returns a line that follows the road
Look into standard pathfinding algorithms
Look into any angle pathfinding algorithms
Compare on road to off road path when considering each pair of nodes for candidate edge
Add candidate edges
Use multiple roads instead of single road
Handle case if road is fragmented
Add pole cost if not close to a road
Use boruvka
Consider using generator to return geometries to be better about memory use