Styjun

I'm working on this project in Python where I have a bunch of polygons representing the clear-sky area on the Earth. I also have some polygons representing spots where a satellite looks (they're essentially squares). I'm trying to figure out how much clear area has been covered by the satellite, and I can do that, but it takes forever. Now I'm wondering how to speed up the calculation process.

It is OK to approximate fairly heavily, I sure have already.

1. Initially, I tried to make all the geometry functions myself, (e.g. intersection of a polygon with another, area of polygon) and that was generally going okay, but my code was slow.




2. Discovered the Shapely library. It's pretty great! Thanks to the people involved with that. Now I use Shapely shapes for the polygons, and finding intersections and areas etc. is super easy. However, it's not always super fast.




3. Because the area calculation was taking so long, I have given up on trying to make it accurate, now I am just going for something proportional to area. I am now using the number of points (which are centroids of the initial polygons) covered by a polygon as my value for area. If you can get this to calculate real area in a good time, I will give you extra imaginary points.

I think currently the code is not horrible, and maybe there's not much more I can do as far as optimization, but wanted to throw it out there in case anyone has ideas.

As much as it has saved me some time, use of Shapely is not a requirement for this code.

One idea I have had was that I should take all the small polygons and join them into a mega-polygon and calculate the overlap area with that, but I'm not sure how to do that.

Running on Ubuntu (not that it really matters here, right?), and using Python 3.7. If you'd like any other details I'd be happy to give more info.

`

# for the area stuff we care about
import numpy as np
from shapely.geometry import Polygon, Point

# for timing things
from functools import wraps
from time import time

def timing(f):
 """Prints elapsed time after a function runs"""
 # Thanks, яүυк!
 # https://codereview.stackexchange.com/questions/169870/decorator-to-measure-execution-time-of-a-function
 # I'm aware the timing won't be perfect, this is just to get an idea

 @wraps(f)
 def wrapper(*args, **kwargs):
 start = time()
 result = f(*args, **kwargs)
 end = time()
 print('Elapsed time: '.format(end-start))
 return result
 return wrapper

@timing
def generate_points_polys():
 """I do this in a different way in the actual code, but have simplified
 the approach for this question"""
 # domain of interest
 x_low = -1000
 x_high = 1000
 y_low = -1000
 y_high = 1000
 resolution = 0.1

 # these are the points I want to cover
 x_points = np.arange(x_low, x_high, resolution)
 y_points = np.arange(y_low, y_high, resolution)

 # convert them into shapely points for easier geometrical operations
 points = zip(x_points, y_points)
 points = np.array([Point(x, y) for x, y in points], dtype=object)

 # generate polygons
 polys_per_col = 10
 polys_per_row = 10
 polys = np.empty((polys_per_col, polys_per_row), dtype=object)

 for index, _ in np.ndenumerate(polys):
 # let's say they're 5x5 for the sake of the question
 # the shapes are dynamic in real life
 x = np.random.uniform(low=x_low, high=x_high)
 y = np.random.uniform(low=y_low, high=y_high)
 vertices = [[x,y],[x+5,y],[x+5,y+5],[x,y+5]]
 polys[index] = Polygon(vertices)
 return points, polys

@timing
def calculate_area(points, polys):
 """This is the function we're trying to optimize"""
 areas = np.empty(polys.shape, dtype=float)
 for index, poly in np.ndenumerate(polys):
 # get number of polygons at least partially covered
 num_covered = len([1 for point in points if poly.intersects(point)])
 # area is not equal to this but I figure it's proportional
 # and calculating the real area would be expensive
 areas[index] = num_covered
 return areas

@timing
def do_the_things():
 points, polys = generate_points_polys()
 calculate_area(points, polys)

do_the_things()

`

I would like to be able to calculate the total area quickly, but it is taking a long time using the methods I've tried. Any ideas of code or different approximations I should use?

EDIT:

I have tried one of the suggestions in the comments, the one about using a MultiPolygon. Here are a couple trials involving that kind of thing, using the same number of points:

`

# Function 1, finds how many points are covered
def calculate_area(points, polys):
 areas = np.empty(polys.shape, dtype=float)
 for index, poly in np.ndenumerate(polys):
 num_covered = len([1 for point in points if poly.intersects(point)])
 areas[index] = num_covered
 return areas

# Function 2, number of background polygons (b_polys) covered
def calculate_area2(b_polys, polys):
 areas = np.empty(polys.shape, dtype=float)
 for index, poly in np.ndenumerate(polys):
 close_polys = [p for p in b_polys if poly.intersects(p)]
 areas[index] = len(close_polys)
 return areas

# Function 3, Function 2, but with actual areas
def calculate_area3(b_polys, polys):
 areas = np.empty(polys.shape, dtype=float)
 for index, poly in np.ndenumerate(polys):
 close_polys = [p for p in b_polys if poly.intersects(p)]
 total_area = 0
 for p in close_polys:
 total_area += p.area
 areas[index] = total_area
 return areas

# Function 4, calculates overlap area with cascaded_union
def calculate_area4(b_polys, polys):
 areas = np.empty(polys.shape, dtype=float)
 for index, poly in np.ndenumerate(polys):
 area = shapely.ops.cascaded_union(list(b_polys)+[poly]).area
 areas[index] = area
 return areas

`

My findings (using the same timer above and the same number of points/polygons for each function):

`

Elapsed time (s)
1: 0.8620626926422119
2: 0.8469529151916504
3: 1.2330029010772705
4: 2.49029541015625

`