# Prototype Solar Mini Grid Layout (1st Implementation)

 Pay Notebook Creator: Roy Hyunjin Han 250 Set Container: Numerical CPU with TINY Memory for 10 Minutes 0 Total 0

# Vision¶

Inspire a community of professionals skilled in spatiotemporal analysis for community health and safety.

# Mission¶

Prototype an algorithm that places distribution poles to connect service drop poles.

Roy Hyunjin Han

# Context¶

The next step is to connect the service drop poles in a mini grid that minimizes the length of distribution line.

# Timeframe¶

20171004-2100 - 20171018-2100: 2 weeks estimated

# Objectives¶

Use brute force kruskal
Update brute force kruskal with kdtree

# Log¶

20171004-2100 - 20171004-2200: 60 minutes

Case 1: Without roads, without unpassable obstacles

• Get all possible lines connecting pairs of nodes
• Use minimum spanning tree

Case 2: Without roads, with unpassable obstacles

• Get all possible lines connecting pairs of nodes
• If a line intersects with an unpassable obstacle, modify line to follow perimeter of buffer around obstacle

Case 3: With roads, without unpassable obstacles

• Run minimum spanning tree on network

Case 4: With roads, with unpassable obstacles

• Run as in case 3 but wrap lines around perimeter of a buffer around an obstacle

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.

In :
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.

In :
import numpy as np
from itertools import combinations
from scipy.spatial.distance import euclidean as get_distance

points = np.random.rand(10, 2)

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

In :
# Cut line into separated points (convert into dotted line)
from shapely.geometry import LineString
line = LineString([(0, 0), (1, 1), (2, 0)])
line

Out:
<shapely.geometry.linestring.LineString at 0x7f0a4990ef28>
In :
list(line.coords)

Out:
[(0.0, 0.0), (1.0, 1.0), (2.0, 0.0)]
In :
print(line.interpolate(0.5))

POINT (0.3535533905932737 0.3535533905932737)

In :
print(line.interpolate(1))

POINT (0.7071067811865475 0.7071067811865475)

In :
print(line.interpolate(2))

POINT (1.414213562373095 0.5857864376269051)

In :
line.length

Out:
2.8284271247461903
In :
print(line.interpolate(line.length))

POINT (2 0)

In :
print(line.interpolate(line.length - 0.1))

POINT (1.929289321881345 0.07071067811865483)

In :
line = LineString([(0, 0), (2, 2), (2, 0)])
line

Out:
<shapely.geometry.linestring.LineString at 0x7f0a498ac160>
In :
from shapely.geometry import MultiPoint
MultiPoint(line.coords)

Out:
<shapely.geometry.multipoint.MultiPoint at 0x7f0a498ac908>
In :
maximum_distribution_length = 1
maximum_distribution_length

Out:
1
In :
line.coords[0:2]

Out:
[(0.0, 0.0), (2.0, 2.0)]
In :
LineString(line.coords[0:2]).length

Out:
2.8284271247461903
In :
LineString(line.coords[0:2]).length / maximum_distribution_length

Out:
2.8284271247461903
In :
int(LineString(line.coords[0:2]).length / maximum_distribution_length) + 1

Out:
3
In :
LineString(line.coords[0:2]).length / 3.

Out:
0.9428090415820635
In :
from shapely.geometry import MultiPoint
MultiPoint([
line.coords,
line.interpolate(0.9428090415820635),
line.interpolate(0.9428090415820635 * 2),
line.coords,
])

Out:
<shapely.geometry.multipoint.MultiPoint at 0x7f0a4990eb38>

20171016-0800 - 20171016-0930: 90 minutes (Salah, Roy)

In :
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])

Out:
<shapely.geometry.collection.GeometryCollection at 0x7f1ce0615be0>
In :
list(line.coords)

Out:
[(-1.0, 0.0), (8.0, 3.0)]
In :
list(line.interpolate(0).coords)

Out:
[(-1.0, 0.0)]
In :
list(line.interpolate(1).coords)

Out:
[(-0.05131670194948623, 0.31622776601683794)]
In :
list(line.interpolate(5).coords)

Out:
[(3.743416490252569, 1.5811388300841898)]
In :
list(line.interpolate(line.length).coords)

Out:
[(8.0, 3.0)]
In :
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.0, 0.0)
(5.0, 5.0)
(10.0, 10.0)
(5.0, 5.0)
(0.0, 0.0)

In :
GeometryCollection(points + [line])

Out:
<shapely.geometry.collection.GeometryCollection at 0x7f1cc77be470>
In :
GeometryCollection(points + [
line.interpolate(line.project(points)),
])

Out:
<shapely.geometry.collection.GeometryCollection at 0x7f1cc77b9978>
In :
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.edges(data=True)

Out:
EdgeDataView([(0, 1, {'weight': 3.1622776601683795}), (0, 2, {'weight': 5.656854249492381}), (0, 3, {'weight': 7.615773105863909}), (0, 4, {'weight': 8.246211251235321}), (1, 2, {'weight': 7.615773105863909}), (1, 3, {'weight': 6.0}), (1, 4, {'weight': 8.602325267042627}), (2, 3, {'weight': 7.615773105863909}), (2, 4, {'weight': 4.47213595499958}), (3, 4, {'weight': 5.0990195135927845})])
In :
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).coords)

Out:
[(0.8, 0.6000000000000001)]

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)

wkt
"LINESTRING (-1 0, 8 3)"
In :
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]})

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)

In :
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

# indices = range(len(nodes))
g = nx.Graph()
for point1, point2 in combinations(nodes, 2):
distance = get_distance(point1, point2)
print(len(list(combinations(ds, 2))))
# roads should have discounted weight
for n in nodes:
distance = get_distance(xy, n)

print(len(g.edges(data=True)))
pprint(g.edges(data=True))

# get edges between nodes on road

# 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

pprint(mst)

10
10
10
10
10
10
10
10
10
10
15
EdgeDataView([(((0.0, 3.0),), ((1.0, 0.0),), {'weight': 6.324555320336759}), (((0.0, 3.0),), ((4.0, 7.0),), {'weight': 11.313708498984761}), (((0.0, 3.0),), ((7.0, 0.0),), {'weight': 15.231546211727817}), (((0.0, 3.0),), ((8.0, 5.0),), {'weight': 16.492422502470642}), (((0.0, 3.0),), ((0.8, 0.6000000000000001),), {'weight': 5.059644256269407}), (((1.0, 0.0),), ((4.0, 7.0),), {'weight': 15.231546211727817}), (((1.0, 0.0),), ((7.0, 0.0),), {'weight': 12.0}), (((1.0, 0.0),), ((8.0, 5.0),), {'weight': 17.204650534085253}), (((1.0, 0.0),), ((0.8, 0.6000000000000001),), {'weight': 1.2649110640673518}), (((4.0, 7.0),), ((7.0, 0.0),), {'weight': 15.231546211727817}), (((4.0, 7.0),), ((8.0, 5.0),), {'weight': 8.94427190999916}), (((4.0, 7.0),), ((5.6, 2.1999999999999997),), {'weight': 10.119288512538814}), (((7.0, 0.0),), ((8.0, 5.0),), {'weight': 10.198039027185569}), (((7.0, 0.0),), ((6.2, 2.4000000000000004),), {'weight': 5.059644256269407}), (((8.0, 5.0),), ((8.0, 3.0),), {'weight': 4.0})])
<networkx.classes.graph.Graph object at 0x7f1cb3e6fb70>

+ 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 script

20171018-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.

In :
drop_pole_table_path = 'drop-pole-points.csv'
distribution_pole_interval_in_meters = 1

In :
from pandas import read_csv
drop_pole_table

Out:
<style> .dataframe thead tr:only-child th { text-align: right; } .dataframe thead th { text-align: left; } .dataframe tbody tr th { vertical-align: top; } </style>
wkt
0 POINT (0 3)
1 POINT (1 0)
2 POINT (4 7)
3 POINT (7 0)
4 POINT (8 5)
In :
drop_pole_table.columns

Out:
Index(['wkt'], dtype='object')
In :
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)

Out:
<shapely.geometry.multipoint.MultiPoint at 0x7f1cb3e5db00>
In :
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)
tuple(p1.coords),
tuple(p2.coords), weight=distance)
graph.edges(data=True)

Out:
EdgeDataView([((0.0, 3.0), (1.0, 0.0), {'weight': 3.1622776601683795}), ((0.0, 3.0), (4.0, 7.0), {'weight': 5.656854249492381}), ((0.0, 3.0), (7.0, 0.0), {'weight': 7.615773105863909}), ((0.0, 3.0), (8.0, 5.0), {'weight': 8.246211251235321}), ((1.0, 0.0), (4.0, 7.0), {'weight': 7.615773105863909}), ((1.0, 0.0), (7.0, 0.0), {'weight': 6.0}), ((1.0, 0.0), (8.0, 5.0), {'weight': 8.602325267042627}), ((4.0, 7.0), (7.0, 0.0), {'weight': 7.615773105863909}), ((4.0, 7.0), (8.0, 5.0), {'weight': 4.47213595499958}), ((7.0, 0.0), (8.0, 5.0), {'weight': 5.0990195135927845})])
In :
tree = nx.minimum_spanning_tree(graph)
tree.edges(data=True)

Out:
EdgeDataView([((0.0, 3.0), (1.0, 0.0), {'weight': 3.1622776601683795}), ((0.0, 3.0), (4.0, 7.0), {'weight': 5.656854249492381}), ((4.0, 7.0), (8.0, 5.0), {'weight': 4.47213595499958}), ((7.0, 0.0), (8.0, 5.0), {'weight': 5.0990195135927845})])
In :
# 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
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(
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)

(0.0, 3.0) (1.0, 0.0)
POINT (0 3)
POINT (0.25 2.25)
POINT (0.5 1.5)
POINT (0.75 0.75)
(1.0, 0.0)
***
(0.0, 3.0) (4.0, 7.0)
POINT (0 3)
POINT (0.6666666666666667 3.666666666666667)
POINT (1.333333333333333 4.333333333333334)
POINT (2 5)
POINT (2.666666666666667 5.666666666666667)
POINT (3.333333333333333 6.333333333333333)
(4.0, 7.0)
***
(4.0, 7.0) (8.0, 5.0)
POINT (4 7)
POINT (4.8 6.6)
POINT (5.6 6.2)
POINT (6.4 5.8)
POINT (7.2 5.4)
(8.0, 5.0)
***
(7.0, 0.0) (8.0, 5.0)
POINT (7 0)
POINT (7.166666666666667 0.8333333333333333)
POINT (7.333333333333333 1.666666666666667)
POINT (7.5 2.5)
POINT (7.666666666666667 3.333333333333333)
POINT (7.833333333333333 4.166666666666666)
(8.0, 5.0)
***

Out:
<shapely.geometry.multipoint.MultiPoint at 0x7f1cb3e64390>
+ Prototype distribution pole tool
+ Finish line splicing algorithm



20171107-2230 - 20171107-0000: 90 minutes (Salah)

In :
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',
])

Out:
<geotable.ColorfulGeometryCollection at 0x7f5b1406fd68>
In :
MultiPolygon(obstacles).bounds

Out:
(0.0, 1.0, 9.0, 8.0)
In :
import numpy as np
np.arange(0.5, 5.5, 1)

Out:
array([ 0.5,  1.5,  2.5,  3.5,  4.5])
In :
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, obstacle_bounds]
ys = [source_point.y, target_point.y] + [
obstacle_bounds, obstacle_bounds]
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)

Out:
<shapely.geometry.collection.GeometryCollection at 0x7f5b01dc30b8>
In :
eligible_points = []
for point in mesh_points:
for obstacle in obstacles:
if obstacle.contains(point):
break
else:
eligible_points.append(point)
GeometryCollection(eligible_points)

Out:
<shapely.geometry.collection.GeometryCollection at 0x7f5b01d12b70>
In :
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
(x + 1, y),
(x, y + 1),
(x - 1, y),
(x, y - 1),
]
for xy in filter(lambda _: _ in xys, adjacent_xys):
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):

In :
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

Out:
<shapely.geometry.linestring.LineString at 0x7f5af7a62908>
In :
ColorfulGeometryCollection([
GeometryCollection(obstacles),
line,
source_point,
target_point,
], colors=['gray', 'black', 'red', 'blue'])

Out:
<geotable.ColorfulGeometryCollection at 0x7f5af6d60438>

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 tasks

## Future¶

Write 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
Consider using generator to return geometries to be better about memory use