docs: def _build_voronoi_polygons(df):
      """
      Given a GeoDataFrame of point geometries and pre-computed plot extrema, build Voronoi simplexes for the given
      points in the given space and returns them.
  
      Voronoi simplexes which are located on the edges of the graph may extend into infinity in some direction. In
      other words, the set of points nearest the given point does not necessarily have to be a closed polygon. We force
      these non-hermetic spaces into polygons using a subroutine.
  
      Parameters
      ----------
      df : GeoDataFrame instance
          The `GeoDataFrame` of points being partitioned.
  
      Returns
      -------
      polygons : list of shapely.geometry.Polygon objects
          The Voronoi polygon output.
      """
      from scipy.spatial import Voronoi
      geom = np.array(df.geometry.map(lambda p: [p.x, p.y]).tolist())
      vor = Voronoi(geom)
  
      polygons = []
  
      for idx_point, _ in enumerate(vor.points):
          idx_point_region = vor.point_region[idx_point]
          idxs_vertices = np.array(vor.regions[idx_point_region])
  
          is_finite = not np.any(idxs_vertices == -1)
  
          if is_finite:
              # Easy case, the region is closed. Make a polygon out of the Voronoi ridge points.
              idx_point_region = vor.point_region[idx_point]
              idxs_vertices = np.array(vor.regions[idx_point_region])
              region_vertices = vor.vertices[idxs_vertices]
              region_poly = shapely.geometry.Polygon(region_vertices)
  
              polygons.append(region_poly)
  
          else:
              # Hard case, the region is open. Project new edges out to the margins of the plot.
              # See `scipy.spatial.voronoi_plot_2d` for the source of this calculation.
              point_idx_ridges_idx = np.where((vor.ridge_points == idx_point).any(axis=1))[0]
  
              # TODO: why does this happen?
              if len(point_idx_ridges_idx) == 0:
                  continue
  
              ptp_bound = vor.points.ptp(axis=0)
              center = vor.points.mean(axis=0)
  
              finite_segments = []
              infinite_segments = []
  
              pointwise_ridge_points = vor.ridge_points[point_idx_ridges_idx]
              pointwise_ridge_vertices = np.asarray(vor.ridge_vertices)[point_idx_ridges_idx]
  
              for pointidx, simplex in zip(pointwise_ridge_points, pointwise_ridge_vertices):
                  simplex = np.asarray(simplex)
  
                  if np.all(simplex >= 0):
                      finite_segments.append(vor.vertices[simplex])
  
                  else:
                      i = simplex[simplex >= 0][0]  # finite end Voronoi vertex
  
                      t = vor.points[pointidx[1]] - vor.points[pointidx[0]]  # tangent
                      t /= np.linalg.norm(t)
                      n = np.array([-t[1], t[0]])  # normal
  
                      midpoint = vor.points[pointidx].mean(axis=0)
                      direction = np.sign(np.dot(midpoint - center, n)) * n
                      far_point = vor.vertices[i] + direction * ptp_bound.max()
  
                      infinite_segments.append(np.asarray([vor.vertices[i], far_point]))
  
              finite_segments = finite_segments if finite_segments else np.zeros(shape=(0,2,2))
              ls = np.vstack([np.asarray(infinite_segments), np.asarray(finite_segments)])
  
              # We have to trivially sort the line segments into polygonal order. The algorithm that follows is
              # inefficient, being O(n^2), but "good enough" for this use-case.
              ls_sorted = []
  
              while len(ls_sorted) < len(ls):
                  l1 = ls[0] if len(ls_sorted) == 0 else ls_sorted[-1]
                  matches = []
  
                  for l2 in [l for l in ls if not (l == l1).all()]:
                      if np.any(l1 == l2):
                          matches.append(l2)
                      elif np.any(l1 == l2[::-1]):
                          l2 = l2[::-1]
                          matches.append(l2)
  
                  if len(ls_sorted) == 0:
                      ls_sorted.append(l1)
  
                  for match in matches:
                      # in list sytax this would be "if match not in ls_sorted"
                      # in numpy things are more complicated...
                      if not any((match == ls_sort).all() for ls_sort in ls_sorted):
                          ls_sorted.append(match)
                          break
  
              # Build and return the final polygon.
              polyline = np.vstack(ls_sorted)
              geom = shapely.geometry.Polygon(polyline).convex_hull
              polygons.append(geom)
  
      return polygons