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

Prototype an algorithm that places service drop poles near customers.

Roy Hyunjin Han

Placing service drop poles is the first step to connecting a village. It is also possible that we can abstract this algorithm to similar problem definitions.

20170912-1715 - 20170913-1715: 1 day estimated

20170912-1715 - 20170919-2215: 1 week actual

```
+ Consider different ways to implement algorithm
+ Implement at least one version
```

20170912-2130 - 20170912-2200: 30 minutes

While walking home, I thought of a potential solution.

- Let's take each point (or geometry) that represents a customer. Then we'll draw a buffer around the geometry. That buffer represents the locus of possible locations of a service drop pole for that customer.
- Then let's take the intersection of all the buffers. If there are multiple buffers, there there are likely to be multiple areas of intersection. For each area of intersection, we'll count how many times it intersects with a buffer, which indicates the number of customers a service drop pole in that area will be able to reach.
- Then we'll start with the intersection area that has the most number of overlapping buffers and put as many service drop poles in it until it satisfies its customers.
- Then we'll go through the intersection areas in order, from the most number of overlapping buffers to the least until we've satisfied them all.
- Each time we satisfy a customer, we'll remove it from the unsatisfied customers list.
- Finally, we'll have the unsatisfied buffers, which we'll keep.
- The end result is a set of geometries, preserving a degree of freedom on where exactly to place each service drop pole for future steps.

We should refine these thoughts in more precise pseudocode.

`+ Draft potential solution in words`

- Start with a set of geometries that represent each customer.
- Draw a buffer around each geometry that represents the possible locations of a service drop pole for that customer.
- Compute the intersection of all buffers.
- Split the buffer intersection into separate geometries which we will call buffer intersections.
- For each buffer intersection, loop through the buffers and count the number of times the buffer intersection is contained in a buffer.
- Sort the buffer intersections in order of most number of overlapping buffers.
- For each buffer intersection, assign the minimum number of service drop poles that will cover the customers for that buffer intersection. We'll call these social service drop poles.
- Remove each satisfied buffer from the list of unsatisfied buffers.
- Assign a service drop pole to each unsatified buffer, which we will call lonely service drop poles.
- Count the number of service drop poles.
- Output the buffer intersections and unsatisfied buffers, which is our solution.

I think our metric for now is simply the number of service drop poles combined with the sum of distances from customer to service drop pole for each customer.

```
score = A * sum(customer_distances) - B * pole_count
where A is $/length and B is $/pole
+ Draft potential solution in pseudocode
+ Think of different ways to implement algorithm
+ Propose a metric
```

For the toy problem, we'll start with points, though we should confirm that we can get a buffer around any geometry.

In [ ]:

```
from shapely.geometry import Point
point = Point(0, 0)
point
```

In [ ]:

```
point.buffer(1)
```

In [ ]:

```
from shapely.geometry import LineString
line = LineString([(0, 0), (1, 1)])
line
```

In [ ]:

```
line.buffer(1)
```

In [ ]:

```
from shapely.geometry import Polygon
polygon = Polygon([(0, 0), (1, 1), (1, 0)])
polygon
```

In [ ]:

```
polygon.buffer(1)
```

`+ Review use of buffers`

In [ ]:

```
# Start with a set of geometries that represent each customer
from shapely.geometry import MultiPolygon, Polygon
customers = MultiPolygon([
Polygon([(0, 0), (1, 1), (1, 0)]),
Polygon([(0, 0.5), (0, 1.5), (1, 1.5)]),
Polygon([(2.5, 0.5), (2.5, 1), (3, 1), (3, 0.5)]),
])
customers
```

In [ ]:

```
# Draw a buffer around each geometry that represents
# possible locations of a service drop pole for that customer
maximum_drop_length = 0.5
customers_buffered = customers.buffer(maximum_drop_length)
customers_buffered
```

I think we want to keep each customer separate.

In [ ]:

```
buffered_customers = [x.buffer(maximum_drop_length) for x in customers]
MultiPolygon(buffered_customers)
```

In [ ]:

```
buffered_customers[0].intersection(buffered_customers[1])
```

Now we are trying to compute the intersection of all buffers. There seems to be one possible solution here: https://gis.stackexchange.com/questions/187402/how-to-find-the-intersection-areas-of-overlapping-buffer-zones-in-single-shapefi/187499#187499

In [ ]:

```
from shapely.ops import unary_union, polygonize
pieces = list(polygonize(unary_union(buffered_customers)))
print(len(pieces)) # Expect 4 pieces
MultiPolygon(pieces) # This is not what we want
```

In [ ]:

```
# Compute the intersection of all buffers
from shapely.geometry import LineString
from shapely.ops import unary_union, polygonize
rings = [LineString(list(x.exterior.coords)) for x in buffered_customers]
pieces = list(polygonize(unary_union(rings)))
print(len(pieces)) # Expect 4 pieces
MultiPolygon(pieces)
```

```
+ Start with a set of geometries that represent each customer
+ Draw a buffer around each geometry
+ Compute the intersection of all buffers
```

In [ ]:

```
buffered_customers[0]
```

In [ ]:

```
buffered_customers[1]
```

In [ ]:

```
pieces[1]
```

In [ ]:

```
pieces[2]
```

In [ ]:

```
index = 2
print([
buffered_customers[0].contains(pieces[index]),
buffered_customers[0].overlaps(pieces[index]),
buffered_customers[0].intersects(pieces[index]),
buffered_customers[0].touches(pieces[index]),
])
print([
buffered_customers[1].contains(pieces[index]),
buffered_customers[1].overlaps(pieces[index]),
buffered_customers[1].intersects(pieces[index]),
buffered_customers[1].touches(pieces[index]),
])
```

In [ ]:

```
# Isolate buffer intersections
intersections = []
buffers = buffered_customers
for piece in pieces:
membership_count = 0
for buffer in buffers:
if buffer.contains(piece):
membership_count += 1
if membership_count > 1:
piece.membership_count = membership_count
intersections.append(piece)
MultiPolygon(intersections)
```

20170919-2100 - 20170919-2200: 60 minutes

I was trying to isolate the buffer intersections last time. Maybe we can get a random point inside each intersection (such as the centroid) and use that as the test point.

In [ ]:

```
piece.centroid
```

In [ ]:

```
# Compute the intersection of all buffers
from shapely.geometry import LineString
from shapely.ops import unary_union, polygonize
rings = [LineString(list(x.exterior.coords)) for x in buffered_customers]
pieces = list(polygonize(unary_union(rings)))
intersections = []
buffers = buffered_customers
for piece in pieces:
membership_count = 0
for buffer in buffers:
if buffer.contains(piece.centroid):
membership_count += 1
if membership_count > 1:
piece.membership_count = membership_count
intersections.append(piece)
MultiPolygon(intersections)
```

We were able to isolate the areas of intersection.

```
+ Split the buffer intersection into separate geometries which we will call buffer intersections.
+ For each buffer intersection, loop through the buffers and count the number of times the buffer intersection is contained in a buffer.
```

In [ ]:

```
for intersection in intersections:
print(intersection.membership_count)
```

In [ ]:

```
sorted_intersections = sorted(intersections, key=lambda x: -x.membership_count)
sorted_intersections
```

Now we will assign service drop poles for each intersection area.

```
+ Sort the buffer intersections in order of most number of overlapping buffers.
```

The number of service drop poles is membership_count divided by maximum_customer_count_per_drop_pole, rounded up to nearest integer.

`+ For each buffer intersection, assign the minimum number of service drop poles that will cover the customers for that buffer intersection. We'll call these social service drop poles.`

In [ ]:

```
maximum_customer_count_per_drop_pole = 4
membership_count = 10
from math import ceil
ceil(membership_count / maximum_customer_count_per_drop_pole)
```

In [ ]:

```
# Compute the intersection of all buffers
from shapely.geometry import LineString
from shapely.ops import unary_union, polygonize
rings = [LineString(list(x.exterior.coords)) for x in buffered_customers]
pieces = list(polygonize(unary_union(rings)))
import copy
intersections = []
buffers = buffered_customers
lonely_buffers = copy.copy(buffers)
for piece in pieces:
piece.buffers = []
for buffer in buffers:
if buffer.contains(piece.centroid):
piece.buffers.append(buffer)
if len(piece.buffers) > 1:
for buffer in piece.buffers:
lonely_buffers.remove(buffer)
intersections.append(piece)
MultiPolygon(lonely_buffers)
```

In [ ]:

```
drop_pole_count = 0
drop_pole_count += len(lonely_buffers)
drop_pole_count
```

The next step is to process the isolated customers.

```
+ Remove each satisfied buffer from the list of unsatisfied buffers.
+ Assign a service drop pole to each unsatified buffer, which we will call lonely service drop poles.
+ Count the number of service drop poles.
+ Output the buffer intersections and unsatisfied buffers, which is our solution.
```

In [ ]:

```
# Start with a set of geometries that represent each customer
from shapely.geometry import MultiPolygon, Polygon
customers = MultiPolygon([
Polygon([(0, 0), (1, 1), (1, 0)]),
Polygon([(0, 0.5), (0, 1.5), (1, 1.5)]),
Polygon([(2.5, 0.5), (2.5, 1), (3, 1), (3, 0.5)]),
])
customers
```

In [ ]:

```
maximum_drop_length = 0.5
# Compute buffers
buffers = [x.buffer(maximum_drop_length) for x in customers]
# Compute the intersection of all buffers
from shapely.geometry import LineString
from shapely.ops import unary_union, polygonize
rings = [LineString(list(x.exterior.coords)) for x in buffers]
pieces = list(polygonize(unary_union(rings)))
import copy
intersections = []
lonely_buffers = copy.copy(buffers)
for piece in pieces:
piece.buffers = []
for buffer in buffers:
if buffer.contains(piece.centroid):
piece.buffers.append(buffer)
if len(piece.buffers) > 1:
for buffer in piece.buffers:
lonely_buffers.remove(buffer)
intersections.append(piece)
from math import ceil
maximum_customer_count_per_drop_pole = 4
drop_pole_count = 0
for intersection in intersections:
drop_pole_count += ceil(len(intersection.buffers) / maximum_customer_count_per_drop_pole)
drop_pole_count += len(lonely_buffers)
print('drop_pole_count = %s' % drop_pole_count)
```

In [ ]:

```
MultiPolygon(intersections)
```

In [ ]:

```
MultiPolygon(lonely_buffers)
```

`+ Combine the steps into a single algorithm`

20170919-2200 - 20170919-2215: 15 minutes

In [ ]:

```
maximum_customer_count_per_drop_pole = 4
maximum_drop_length = 0.5
drop_pole_count = 0
from copy import copy
from math import ceil
from shapely.geometry import LineString
from shapely.ops import unary_union, polygonize
def get_disjoint_polygons(overlapping_polygons):
rings = [LineString(list(
x.exterior.coords)) for x in overlapping_polygons]
return list(polygonize(unary_union(rings)))
# Compute buffers
buffers = [x.buffer(maximum_drop_length) for x in customers]
lonely_buffers = copy(buffers)
# Identify social service drop poles
intersections = []
for polygon in get_disjoint_polygons(buffers):
polygon.buffers = []
for buffer in buffers:
if buffer.contains(polygon.centroid):
polygon.buffers.append(buffer)
if len(polygon.buffers) > 1:
for buffer in polygon.buffers:
lonely_buffers.remove(buffer)
intersections.append(polygon)
# Count service drop poles
for intersection in intersections:
drop_pole_count += ceil(len(
intersection.buffers) / maximum_customer_count_per_drop_pole)
drop_pole_count += len(lonely_buffers)
drop_pole_count
```

Okay, now we have one version.

```
+ Clean code
+ Implement at least one version
```

I think we can call service drop poles as drop poles for short.

20171116-0100 - 20171116-0130: 30 minutes

Deril noted that the service drop pole algorithm is not solving the triangle problem correctly. It seems like I forgot the step of sorting the intersections by connection_count.

In [13]:

```
from copy import copy
from shapely.geometry import LineString
from shapely.ops import polygonize, unary_union
def get_source_polygons(target_geometries, maximum_distance):
# Convert target_geometries into target_polygons using maximum_distance
target_polygons = [x.buffer(maximum_distance) for x in target_geometries]
# Identify overlapping areas
sliced_polygons = get_disjoint_polygons(target_polygons)
for polygon in sliced_polygons:
candidate_polygons = []
for target_polygon in target_polygons:
if target_polygon.contains(polygon.centroid):
candidate_polygons.append(target_polygon)
polygon.candidate_polygons = candidate_polygons
# Sort overlapping areas by overlap count
sorted_polygons = sorted(
sliced_polygons, key=lambda x: -len(x.candidate_polygons))
# Assign target_polygons to each sorted_polygon
social_polygons, lonely_polygons = [], []
for polygon in sorted_polygons:
connected_polygons = []
for candidate_polygon in polygon.candidate_polygons:
try:
target_polygons.remove(candidate_polygon)
except ValueError:
continue
connected_polygons.append(candidate_polygon)
connection_count = len(connected_polygons)
if connection_count > 1:
social_polygons.append(polygon)
elif connection_count == 1:
lonely_polygons.append(polygon)
polygon.connected_polygons = connected_polygons
return social_polygons + lonely_polygons
def get_disjoint_polygons(overlapping_polygons):
'Split overlapping polygons into disjoint polygons'
rings = [LineString(list(
x.exterior.coords)) for x in overlapping_polygons]
return list(polygonize(unary_union(rings)))
```

In [19]:

```
from shapely.geometry import GeometryCollection, Point
points = [
Point(-1, 0),
Point(+1, 0),
Point(0, +2),
]
source_polygons = get_source_polygons(points, 1.5)
GeometryCollection(points + source_polygons)
```

Out[19]: