I am trying to efficiently "partition" GeoJSON objects into aligned square tiles of any size, such that:
- The entire polygon or multi-polygon is covered (there is no area in the polygon that doesn't have a tile covering it).
- There is no tile that doesn't cover any area of the polygon.
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:
Or leave entire horizontal areas out:
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);