title: What is the most efficient way to find total area covered by a bunch of polygons?
body: <p>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.</p>
  
  <p>It is OK to approximate fairly heavily, I sure have already.</p>
  
  <ol>
  <li><p>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.</p></li>
  <li><p>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.</p></li>
  <li><p>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.</p></li>
  </ol>
  
  <p>I think currently the code is not horrible, and maybe there's not <em>much</em> more I can do as far as optimization, but wanted to throw it out there in case anyone has ideas.</p>
  
  <p>As much as it has saved me some time, use of Shapely is not a requirement for this code.</p>
  
  <p>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.</p>
  
  <p>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.</p>
  
  <p>`</p>
  
  <pre><code># 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()
  </code></pre>
  
  <p>`</p>
  
  <p>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?</p>
  
  <p>EDIT: </p>
  
  <p>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:</p>
  
  <p>`</p>
  
  <pre><code># 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
  </code></pre>
  
  <p>`</p>
  
  <p>My findings (using the same timer above and the same number of points/polygons for each function):</p>
  
  <p>`</p>
  
  <pre><code>Elapsed time (s)
  1: 0.8620626926422119
  2: 0.8469529151916504
  3: 1.2330029010772705
  4: 2.49029541015625
  </code></pre>
  
  <p>`</p>