0

I am trying to efficiently "partition" GeoJSON objects into aligned square tiles of any size, such that:

  1. The entire polygon or multi-polygon is covered (there is no area in the polygon that doesn't have a tile covering it).
  2. There is no tile that doesn't cover any area of the polygon.

example

I have tried using libraries like Turfjs's square grid (https://turfjs.org/docs/#squareGrid), but the mask operation is unreasonably slow - it takes minutes for large areas with a low square size. It also doesn't cover the entire area of the polygon - only the interior.

I am trying to use the scanline fill algorithm to do this, since my research has shown it is incredibly fast and catered to this problem, but it doesn't always cover the entire area - it only covers the interior, and will sometimes leave corners out:

missing corners

Or leave entire horizontal areas out:

missing horizontal area

My scanline works like a normal scanline fill, but works based on floats instead of integers. It floors the initial x position such that it will be on a grid line (Math.floor((x - polygonLeftBoundary) / cellWidth) * cellWidth)

Here is my implementation (using TypeScript):

public partitionIntoGrid(polygon, cellSizeKm) {
    const bounds = this.getBoundingBox(polygon);
    const left = bounds[0];
    const cellDistance = this.distanceToDegrees(cellSizeKm);

    const grid = [];

    if (polygon.geometry.type === 'Polygon') {
      polygon = [polygon.geometry.coordinates];
    } else {
      polygon = polygon.geometry.coordinates;
    }

    polygon.forEach(poly => {
      poly = poly[0];

      const edges = poly.reduce((acc, vertex, vertexIndex) => {
        acc.push([vertex, poly[vertexIndex === poly.length - 1 ? 0 : vertexIndex + 1]]);
        return acc;
      }, []);

      let edgeTable = edges
        .map(edge => {
          const x = Math.floor((edge[0][1] < edge[1][1] ? edge[0][0] : edge[1][0]) / cellDistance) * cellDistance;
          return {
            yMin: Math.min(edge[0][1], edge[1][1]), // minumum y coordinate
            yMax: Math.max(edge[0][1], edge[1][1]), // maximum y coordinate
            originalX: x,
            x, // lower coordinate's x
            w: (edge[0][0] - edge[1][0]) / (edge[0][1] - edge[1][1]), // change of edge per y
          };
        })
        .filter(edge => !isNaN(edge.w))
        .sort((a, b) => a.yMax - b.yMax);

      let activeEdgeTable = [];

      let y = edgeTable[0].yMin;
      const originalY = y;

      while (edgeTable.length || activeEdgeTable.length) {
        // move edges from edge table under y into the active edge table
        edgeTable = edgeTable.filter(edge => {
          if (edge.yMin <= y) {
            activeEdgeTable.push(edge);
            return false;
          }
          return true;
        });

        // remove edges from the active edge table whose yMax is smaller than y
        activeEdgeTable = activeEdgeTable.filter(edge => edge.yMax >= y);
        // sort active edge table by x
        activeEdgeTable = activeEdgeTable.sort((a, b) => a.x - b.x);

        // fill pixels between even and odd adjacent pairs of intersections in active edge table
        const pairs = Array.from({ length: activeEdgeTable.length / 2 }, (v, k) => [
          activeEdgeTable[k * 2], activeEdgeTable[k * 2 + 1],
        ]);
        pairs.forEach(pair => {
          for (let x = pair[0].x; x <= pair[1].x; x += cellDistance) {
            grid.push(bboxPolygon([
              x, // minX
              y, // minY
              x + cellDistance, // maxX
              y + cellDistance, // maxY
            ]));
          }
        });

        // increment y
        y += cellDistance;

        // update x for all edges in active edge table
        activeEdgeTable.forEach(edge => edge.x += edge.w * cellDistance);
        // activeEdgeTable.forEach(edge => edge.x = edge.originalX + Math.floor((edge.w * (y - originalY) / cellDistance)) * cellDistance);
      }
    });

    return grid;
  }

I have been attempting to make the addition of the gradient work, but not getting there yet, hence the line is commented:

activeEdgeTable.forEach(edge => edge.x = edge.originalX + Math.floor((edge.w * (y - originalY) / cellDistance)) * cellDistance);
mevens
  • 341
  • 3
  • 8
  • If you know the bounding box of the polygon, then it might be simpler to employ the standard technique for determining whether a point is inside a polygon. ( See https://stackoverflow.com/questions/22521982/check-if-point-inside-a-polygon ). The idea being that you test each of the partition corners to determine if it's inside the polygon, and then if any of the four partition corners (that make up a partition) are inside, the partition is considered part of the polygon. – Trentium Jun 18 '19 at 03:44
  • An after thought... Note that testing the corners of the partitions implies the granularity of the check. Eg, the polygon could dip into the partition without actually including the partition corners, and the corners test will yield a negative. If instead checking the partition edges for intersection with the polygon, you can still run into a situation where the entire polygon is inside the partition, and again, the partition check will come up a false negative. So the size of the partition, and more importantly, which points are being sampled, dictates the accuracy of the results. – Trentium Jun 18 '19 at 12:12

0 Answers0