content: from toolz import memoize
  import numpy as np
  
  from datashader.glyphs.line import _build_map_onto_pixel_for_line
  from datashader.glyphs.points import _GeometryLike
  from datashader.utils import ngjit
  
  try:
      import spatialpandas
  except Exception:
      spatialpandas = None
  
  
  class GeopandasPolygonGeom(_GeometryLike):
      @property
      def geom_dtypes(self):
          from geopandas.array import GeometryDtype
          return (GeometryDtype,)
  
      @memoize
      def _build_extend(self, x_mapper, y_mapper, info, append, _antialias_stage_2,
                        _antialias_stage_2_funcs):
          expand_aggs_and_cols = self.expand_aggs_and_cols(append)
          map_onto_pixel = _build_map_onto_pixel_for_line(x_mapper, y_mapper)
          draw_polygon = _build_draw_polygon(
              append, map_onto_pixel, x_mapper, y_mapper, expand_aggs_and_cols
          )
  
          perform_extend_cpu = _build_extend_geopandas_polygon_geometry(
              draw_polygon, expand_aggs_and_cols
          )
          geom_name = self.geometry
  
          def extend(aggs, df, vt, bounds, plot_start=True):
              sx, tx, sy, ty = vt
              xmin, xmax, ymin, ymax = bounds
              aggs_and_cols = aggs + info(df, aggs[0].shape[:2])
              geom_array = df[geom_name].array
  
              perform_extend_cpu(
                  sx, tx, sy, ty,
                  xmin, xmax, ymin, ymax,
                  geom_array, *aggs_and_cols
              )
  
          return extend
  
  
  class PolygonGeom(_GeometryLike):
      # spatialpandas must be available if a PolygonGeom object is created.
  
      @property
      def geom_dtypes(self):
          from spatialpandas.geometry import PolygonDtype, MultiPolygonDtype
          return PolygonDtype, MultiPolygonDtype
  
      @memoize
      def _build_extend(self, x_mapper, y_mapper, info, append, _antialias_stage_2,
                        _antialias_stage_2_funcs):
          expand_aggs_and_cols = self.expand_aggs_and_cols(append)
          map_onto_pixel = _build_map_onto_pixel_for_line(x_mapper, y_mapper)
          draw_polygon = _build_draw_polygon(
              append, map_onto_pixel, x_mapper, y_mapper, expand_aggs_and_cols
          )
  
          perform_extend_cpu = _build_extend_polygon_geometry(
              draw_polygon, expand_aggs_and_cols
          )
          geom_name = self.geometry
  
          def extend(aggs, df, vt, bounds, plot_start=True):
              sx, tx, sy, ty = vt
              xmin, xmax, ymin, ymax = bounds
              aggs_and_cols = aggs + info(df, aggs[0].shape[:2])
              geom_array = df[geom_name].array
  
              perform_extend_cpu(
                  sx, tx, sy, ty,
                  xmin, xmax, ymin, ymax,
                  geom_array, *aggs_and_cols
              )
  
          return extend
  
  
  def _build_draw_polygon(append, map_onto_pixel, x_mapper, y_mapper, expand_aggs_and_cols):
      @ngjit
      @expand_aggs_and_cols
      def draw_polygon(
              i, sx, tx, sy, ty, xmin, xmax, ymin, ymax,
              offsets, offset_multiplier, values, xs, ys, yincreasing, eligible,
              *aggs_and_cols
      ):
          """Draw a polygon using a winding-number scan-line algorithm
          """
          # Initialize values of pre-allocated buffers
          xs.fill(np.nan)
          ys.fill(np.nan)
          yincreasing.fill(0)
          eligible.fill(1)
  
          # First pass, compute bounding box of polygon vertices in data coordinates
          start_index = offset_multiplier*offsets[0]
          stop_index = offset_multiplier*offsets[-1]
          # num_edges = stop_index - start_index - 2
          poly_xmin = np.min(values[start_index:stop_index:2])
          poly_ymin = np.min(values[start_index + 1:stop_index:2])
          poly_xmax = np.max(values[start_index:stop_index:2])
          poly_ymax = np.max(values[start_index + 1:stop_index:2])
  
          # skip polygon if outside viewport
          if (poly_xmax < xmin or poly_xmin > xmax
                  or poly_ymax < ymin or poly_ymin > ymax):
              return
  
          # Compute pixel bounds for polygon
          startxi, startyi = map_onto_pixel(
              sx, tx, sy, ty, xmin, xmax, ymin, ymax,
              max(poly_xmin, xmin), max(poly_ymin, ymin)
          )
          stopxi, stopyi = map_onto_pixel(
              sx, tx, sy, ty, xmin, xmax, ymin, ymax,
              min(poly_xmax, xmax), min(poly_ymax, ymax)
          )
          stopxi += 1
          stopyi += 1
  
          # Handle subpixel polygons (pixel width and/or height of polygon is 1)
          if (stopxi - startxi) == 1 and (stopyi - startyi) == 1:
              append(i, startxi, startyi, *aggs_and_cols)
              return
          elif (stopxi - startxi) == 1:
              for yi in range(min(startyi, stopyi), max(startyi, stopyi)):
                  append(i, startxi, yi, *aggs_and_cols)
              return
          elif (stopyi - startyi) == 1:
              for xi in range(min(startxi, stopxi), max(startxi, stopxi)):
                  append(i, xi, startyi, *aggs_and_cols)
              return
  
          # Build arrays of edges in canvas coordinates
          ei = 0
          for j in range(len(offsets) - 1):
              start = offset_multiplier*offsets[j]
              stop = offset_multiplier*offsets[j + 1]
              for k in range(start, stop - 2, 2):
                  x0 = values[k]
                  y0 = values[k + 1]
                  x1 = values[k + 2]
                  y1 = values[k + 3]
  
                  # Map to canvas coordinates without rounding
                  x0c = x_mapper(x0) * sx + tx - 0.5
                  y0c = y_mapper(y0) * sy + ty - 0.5
                  x1c = x_mapper(x1) * sx + tx - 0.5
                  y1c = y_mapper(y1) * sy + ty - 0.5
  
                  if y1c > y0c:
                      xs[ei, 0] = x0c
                      ys[ei, 0] = y0c
                      xs[ei, 1] = x1c
                      ys[ei, 1] = y1c
                      yincreasing[ei] = 1
                  elif y1c < y0c:
                      xs[ei, 1] = x0c
                      ys[ei, 1] = y0c
                      xs[ei, 0] = x1c
                      ys[ei, 0] = y1c
                      yincreasing[ei] = -1
                  else:
                      # Skip horizontal edges
                      continue
  
                  ei += 1
  
          # Perform scan-line algorithm
          num_edges = ei
          for yi in range(startyi, stopyi):
              # All edges eligible at start of new row
              eligible.fill(1)
              for xi in range(startxi, stopxi):
                  # Init winding number
                  winding_number = 0
                  for ei in range(num_edges):
                      if eligible[ei] == 0:
                          # We've already determined that edge is above, below, or left
                          # of edge for the current pixel
                          continue
  
                      # Get edge coordinates.
                      # Note: y1c > y0c due to how xs/ys were populated
                      x0c = xs[ei, 0]
                      x1c = xs[ei, 1]
                      y0c = ys[ei, 0]
                      y1c = ys[ei, 1]
  
                      # Reject edges that are above, below, or left of current pixel.
                      # Note: Edge skipped if lower vertex overlaps,
                      #       but is kept if upper vertex overlaps
                      if (y0c >= yi or y1c < yi
                              or (x0c < xi and x1c < xi)
                      ):
                          # Edge not eligible for any remaining pixel in this row
                          eligible[ei] = 0
                          continue
  
                      if xi <= x0c and xi <= x1c:
                          # Edge is fully to the right of the pixel, so we know ray to the
                          # the right of pixel intersects edge.
                          winding_number += yincreasing[ei]
                      else:
                          # Now check if edge is to the right of pixel using cross product
                          # A is vector from pixel to first vertex
                          ax = x0c - xi
                          ay = y0c - yi
  
                          # B is vector from pixel to second vertex
                          bx = x1c - xi
                          by = y1c - yi
  
                          # Compute cross product of B and A
                          bxa = (bx * ay - by * ax)
  
                          if bxa < 0 or (bxa == 0 and yincreasing[ei]):
                              # Edge to the right
                              winding_number += yincreasing[ei]
                          else:
                              # Edge to left, not eligible for any remaining pixel in row
                              eligible[ei] = 0
                              continue
  
                  if winding_number != 0:
                      # If winding number is not zero, point
                      # is inside polygon
                      append(i, xi, yi, *aggs_and_cols)
  
      return draw_polygon
  
  
  def _build_extend_geopandas_polygon_geometry(
          draw_polygon, expand_aggs_and_cols
  ):
      # Lazy import shapely. Cannot get here if geopandas and shapely are not available.
      import shapely
  
      def extend_cpu(
              sx, tx, sy, ty, xmin, xmax, ymin, ymax, geometry, *aggs_and_cols
      ):
          ragged = shapely.to_ragged_array(geometry)
          geometry_type = ragged[0]
          if geometry_type not in (shapely.GeometryType.POLYGON, shapely.GeometryType.MULTIPOLYGON):
              raise ValueError(
                  "Canvas.polygons supports GeoPandas geometry types of POLYGON and MULTIPOLYGON, "
                  f"not {repr(geometry_type)}")
  
          coords = ragged[1].ravel()
          if geometry_type == shapely.GeometryType.MULTIPOLYGON:
              offsets0, offsets1, offsets2 = ragged[2]
          else:  # POLYGON
              offsets0, offsets1 = ragged[2]
              offsets2 = np.arange(len(offsets1))
  
          extend_cpu_numba(
              sx, tx, sy, ty, xmin, xmax, ymin, ymax,
              coords, offsets0, offsets1, offsets2, *aggs_and_cols
          )
  
      @ngjit
      @expand_aggs_and_cols
      def extend_cpu_numba(sx, tx, sy, ty, xmin, xmax, ymin, ymax,
                           coords, offsets0, offsets1, offsets2, *aggs_and_cols):
          # Pre-allocate temp arrays
          max_edges = 0
          n_multipolygons = len(offsets2) - 1
          for i in range(n_multipolygons):
              polygon_inds = offsets1[offsets2[i]:offsets2[i + 1] + 1]
              for j in range(len(polygon_inds) - 1):
                  start = offsets0[polygon_inds[j]]
                  stop = offsets0[polygon_inds[j + 1]]
                  max_edges = max(max_edges, stop - start - 1)
  
          xs = np.full((max_edges, 2), np.nan, dtype=np.float32)
          ys = np.full((max_edges, 2), np.nan, dtype=np.float32)
          yincreasing = np.zeros(max_edges, dtype=np.int8)
  
          # Initialize array indicating which edges are still eligible for processing
          eligible = np.ones(max_edges, dtype=np.int8)
  
          for i in range(n_multipolygons):
              polygon_inds = offsets1[offsets2[i]:offsets2[i + 1] + 1]
              for j in range(len(polygon_inds) - 1):
                  start = polygon_inds[j]
                  stop = polygon_inds[j + 1]
                  draw_polygon(i, sx, tx, sy, ty, xmin, xmax, ymin, ymax,
                              offsets0[start:stop + 1], 2,
                              coords, xs, ys, yincreasing, eligible, *aggs_and_cols)
  
      return extend_cpu
  
  
  def _build_extend_polygon_geometry(
          draw_polygon, expand_aggs_and_cols
  ):
      def extend_cpu(
              sx, tx, sy, ty, xmin, xmax, ymin, ymax, geometry, *aggs_and_cols
      ):
          values = geometry.buffer_values
          missing = geometry.isna()
          offsets = geometry.buffer_offsets
  
          if geometry._sindex is not None:
              # Compute indices of potentially intersecting polygons using
              # geometry's R-tree if there is one
              eligible_inds = geometry.sindex.intersects((xmin, ymin, xmax, ymax))
          else:
              # Otherwise, process all indices
              eligible_inds = np.arange(0, len(geometry), dtype='uint32')
  
          if len(offsets) == 3:
              # MultiPolygonArray
              offsets0, offsets1, offsets2 = offsets
          else:
              # PolygonArray
              offsets1, offsets2 = offsets
              offsets0 = np.arange(len(offsets1))
  
          extend_cpu_numba(
              sx, tx, sy, ty, xmin, xmax, ymin, ymax,
              values, missing, offsets0, offsets1, offsets2, eligible_inds, *aggs_and_cols
          )
  
      @ngjit
      @expand_aggs_and_cols
      def extend_cpu_numba(
              sx, tx, sy, ty, xmin, xmax, ymin, ymax,
              values, missing, offsets0, offsets1, offsets2,
              eligible_inds, *aggs_and_cols
      ):
          # Pre-allocate temp arrays
          max_edges = 0
          if len(offsets0) > 1:
              for i in eligible_inds:
                  if missing[i]:
                      continue
  
                  polygon_inds = offsets1[offsets0[i]:offsets0[i + 1] + 1]
                  for j in range(len(polygon_inds) - 1):
                      start = offsets2[polygon_inds[j]]
                      stop = offsets2[polygon_inds[j + 1]]
                      max_edges = max(max_edges, (stop - start - 2) // 2)
  
          xs = np.full((max_edges, 2), np.nan, dtype=np.float32)
          ys = np.full((max_edges, 2), np.nan, dtype=np.float32)
          yincreasing = np.zeros(max_edges, dtype=np.int8)
  
          # Initialize array indicating which edges are still eligible for processing
          eligible = np.ones(max_edges, dtype=np.int8)
  
          for i in eligible_inds:
              if missing[i]:
                  continue
  
              polygon_inds = offsets1[offsets0[i]:offsets0[i + 1] + 1]
              for j in range(len(polygon_inds) - 1):
                  start = polygon_inds[j]
                  stop = polygon_inds[j + 1]
  
                  draw_polygon(i, sx, tx, sy, ty, xmin, xmax, ymin, ymax,
                               offsets2[start:stop + 1], 1, values,
                               xs, ys, yincreasing, eligible, *aggs_and_cols)
  
      return extend_cpu
  