docs: def intersect_index_quadrats(gdf, geometry, quadrat_width=0.05, min_num=3, buffer_amount=1e-9):
      """
      Intersect points with a polygon, using an r-tree spatial index and cutting
      the polygon up into smaller sub-polygons for r-tree acceleration.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          the set of points to intersect
      geometry : shapely Polygon or MultiPolygon
          the geometry to intersect with the points
      quadrat_width : numeric
          the linear length (in degrees) of the quadrats with which to cut up the
          geometry (default = 0.05, approx 4km at NYC's latitude)
      min_num : int
          the minimum number of linear quadrat lines (e.g., min_num=3 would
          produce a quadrat grid of 4 squares)
      buffer_amount : numeric
          buffer the quadrat grid lines by quadrat_width times buffer_amount
  
      Returns
      -------
      GeoDataFrame
      """
  
      # create an empty dataframe to append matches to
      points_within_geometry = pd.DataFrame()
  
      # cut the geometry into chunks for r-tree spatial index intersecting
      multipoly = quadrat_cut_geometry(geometry, quadrat_width=quadrat_width, buffer_amount=buffer_amount, min_num=min_num)
  
      # create an r-tree spatial index for the nodes (ie, points)
      start_time = time.time()
      sindex = gdf['geometry'].sindex
      log('Created r-tree spatial index for {:,} points in {:,.2f} seconds'.format(len(gdf), time.time()-start_time))
  
      # loop through each chunk of the geometry to find approximate and then
      # precisely intersecting points
      start_time = time.time()
      for poly in multipoly:
  
          # buffer by the tiny distance to account for any space lost in the
          # quadrat cutting, otherwise may miss point(s) that lay directly on
          # quadrat line
          buffer_size = quadrat_width * buffer_amount
          poly = poly.buffer(buffer_size).buffer(0)
  
          # find approximate matches with r-tree, then precise matches from those
          # approximate ones
          if poly.is_valid and poly.area > 0:
              possible_matches_index = list(sindex.intersection(poly.bounds))
              possible_matches = gdf.iloc[possible_matches_index]
              precise_matches = possible_matches[possible_matches.intersects(poly)]
              points_within_geometry = points_within_geometry.append(precise_matches)
  
      if len(points_within_geometry) > 0:
          # drop duplicate points, if buffered poly caused an overlap on point(s)
          # that lay directly on a quadrat line
          points_within_geometry = points_within_geometry.drop_duplicates(subset='node')
      else:
          # after simplifying the graph, and given the requested network type,
          # there are no nodes inside the polygon - can't create graph from that
          # so throw error
          raise Exception('There are no nodes within the requested geometry')
  
      log('Identified {:,} nodes inside polygon in {:,.2f} seconds'.format(len(points_within_geometry), time.time()-start_time))
      return points_within_geometry