diff --git a/README.md b/README.md index 246942a..dd31dc0 100644 --- a/README.md +++ b/README.md @@ -3,8 +3,28 @@ this is general utilities for usage in Map Colonies project. # included components - [http client](#http-client) + - [footprint to tile ranges](#footprint-to-tile-ranges) # usage + ## footprint to tile ranges + `footprintToTileRanges(footprint, { minZoom, maxZoom })` lazily computes the tile ranges intersecting a footprint (Polygon / MultiPolygon / Feature, holes supported) for every zoom level in the span, on the WGS84 geodetic 2:1 grid. + + it uses hierarchical descent (see `docs/adr/0001-hierarchical-precompute-in-mc-utils.md`): fully-interior subtrees are emitted as compressed ranges with no per-tile work, so memory stays flat regardless of how many tiles the footprint covers. a tile touching the footprint with zero overlap area counts as intersecting. ranges of the same zoom are disjoint; their order is deterministic but unspecified. + + typical use is tile-cache invalidation: expand each range to your cache's key format and batch-delete. + +```typescript +import { footprintToTileRanges } from '@map-colonies/mc-utils'; + +for (const range of footprintToTileRanges(footprint, { minZoom: 0, maxZoom: 18 })) { + for (let x = range.minX; x <= range.maxX; x++) { + for (let y = range.minY; y <= range.maxY; y++) { + actionOnTile(range.zoom, x, y); // e.g. batch UNLINK keys in redis + } + } +} +``` + ## http client this is abstract base class for sending http request with logging and request retries. diff --git a/src/geo/footprintToTileRanges.ts b/src/geo/footprintToTileRanges.ts new file mode 100644 index 0000000..6115f68 --- /dev/null +++ b/src/geo/footprintToTileRanges.ts @@ -0,0 +1,213 @@ +/** + * Hierarchical precompute of the tiles intersecting a footprint. + * The tile pyramid is a quadtree — tile (z, x, y) contains tiles (z+1, 2x..2x+1, 2y..2y+1) - + * so classifying one coarse tile can settle its entire subtree; the walk only recurses where + * the footprint's boundary passes, making total work O(boundary tiles), not O(bbox area). + */ +import type { Position } from 'geojson'; +import { ITileRange } from '../models/interfaces/geo/iTile'; +import { Footprint } from './geoIntersection'; +import { degreesPerTile } from './tiles'; + +interface ZoomSpan { + minZoom: number; + maxZoom: number; +} + +interface Segment { + x1: number; + y1: number; + x2: number; + y2: number; +} + +interface TileBounds { + minLon: number; + minLat: number; + maxLon: number; + maxLat: number; +} + +const MIN_SUPPORTED_ZOOM = 0; +const MAX_SUPPORTED_ZOOM = 22; +const GRID_MIN_LON = -180; +const GRID_MIN_LAT = -90; +const SUBTILES_PER_AXIS = 2; + +// eslint-disable-next-line @typescript-eslint/no-magic-numbers +const midpoint = (a: number, b: number): number => (a + b) / 2; + +// normalizes any accepted footprint into a flat list of rings +function ringsOf(footprint: Footprint): Position[][] { + const geometry = footprint.type === 'Feature' ? footprint.geometry : footprint; + if (geometry.type === 'Polygon') { + return geometry.coordinates; + } + // flattening MultiPolygon parts is safe: edge clipping and the even-odd point test are + // both ring-agnostic, so outer rings, holes and disjoint parts all work with one flat list + return geometry.coordinates.flat(); +} + +// converts rings to individual edges (each ring is closed: last point equals the first) +function ringsToEdges(rings: Position[][]): Segment[] { + const edges: Segment[] = []; + for (const ring of rings) { + for (let i = 0; i < ring.length - 1; i++) { + edges.push({ x1: ring[i][0], y1: ring[i][1], x2: ring[i + 1][0], y2: ring[i + 1][1] }); + } + } + return edges; +} + +// Liang-Barsky line clipping as a boolean: does any part of the segment lie within the tile? +// The segment is P(t) = start + t * (dx, dy) for t in [0, 1]; the loop shrinks the window +// [t0, t1] of the part inside the tile. Touching a border counts as intersecting +// docs: https://mapcolonies.atlassian.net/wiki/spaces/MAPConflicResolution/pages/3334111234/Liang-Barsky+intersection+predicate+segmentIntersectsBounds +function segmentIntersectsBounds(segment: Segment, bounds: TileBounds): boolean { + const dx = segment.x2 - segment.x1; + const dy = segment.y2 - segment.y1; + const p = [-dx, dx, -dy, dy]; // speed toward each boundary: west, east, south, north + const q = [segment.x1 - bounds.minLon, bounds.maxLon - segment.x1, segment.y1 - bounds.minLat, bounds.maxLat - segment.y1]; // start point's signed distance inside each boundary + let t0 = 0; // window start: latest entry crossing seen so far + let t1 = 1; // window end: earliest exit crossing seen so far + + for (let i = 0; i < p.length; i++) { + // docs: #LB-2.1 + if (p[i] === 0) { + // parallel to boundary i, so it never crosses it + if (q[i] < 0) { + // and it runs on the outside of that boundary -> the whole segment is out (docs: #LB-2.1.A (path P1)) + return false; + } + } else { + const r = q[i] / p[i]; // the t at which the segment crosses boundary i + // docs: #LB-2.2 + if (p[i] < 0) { + // moving toward the inside: r is an ENTRY point + if (r > t1) { + return false; // enters only after it already exited elsewhere -> misses the tile (docs: #LB-2.2.A (path P2)) + } + if (r > t0) { + t0 = r; // docs: #LB-2.2.B + } + } else { + // moving toward the outside: r is an EXIT point (docs: #LB-2.3) + if (r < t0) { + return false; // exits before it ever entered -> misses the tile (docs: #LB-2.3.A (path P3)) + } + if (r < t1) { + t1 = r; // docs: #LB-2.3.B + } + } + } + } + // a valid window remains; t0 === t1 means a single touch point, which still counts + return true; // docs: #LB-3 (path P4) +} + +// even-odd (ray casting) point-in-polygon: cast a ray east from the point and count edge +// crossings — odd total means inside. One shared toggle across all rings is what makes +// holes and MultiPolygon parts work with no special casing (inside outer + inside hole +// toggles twice = outside). +function pointInRings(lon: number, lat: number, rings: Position[][]): boolean { + let inside = false; + for (const ring of rings) { + for (let i = 0; i < ring.length - 1; i++) { + const [x1, y1] = ring[i]; + const [x2, y2] = ring[i + 1]; + // half-open test (<= vs >): a shared vertex is counted for exactly one of its two + // edges, and horizontal edges (y1 === y2) never match — both avoid double counting + if ((y1 <= lat && y2 > lat) || (y2 <= lat && y1 > lat)) { + const xCross = x1 + ((lat - y1) * (x2 - x1)) / (y2 - y1); // longitude where the edge crosses the ray's latitude (linear interpolation) + if (xCross > lon) { + // the crossing is east of the point, i.e. actually on the ray + inside = !inside; + } + } + } + } + return inside; +} + +// bounds of a tile on the WGS84 geodetic grid: tiles span 180 / 2^z degrees, origin at the +// south-west corner (-180, -90), x grows east (two hemispheres wide), y grows north +function tileBounds(zoom: number, x: number, y: number): TileBounds { + const span = degreesPerTile(zoom); + const minLon = GRID_MIN_LON + x * span; + const minLat = GRID_MIN_LAT + y * span; + return { minLon, minLat, maxLon: minLon + span, maxLat: minLat + span }; +} + +// emits the whole subtree of a tile proven fully inside the footprint — pure arithmetic, +// no geometry tests: a coarse interior tile can stand in for millions of deep tiles +function* subtreeRanges(zoom: number, x: number, y: number, span: ZoomSpan): Generator { + // Math.max skips levels above the requested span; the block is still derived from THIS node + for (let z = Math.max(zoom, span.minZoom); z <= span.maxZoom; z++) { + const factor = Math.pow(SUBTILES_PER_AXIS, z - zoom); // each level doubles both axes + // the block this tile covers at level z: x' in [x * factor, (x + 1) * factor - 1], y' likewise + yield { zoom: z, minX: x * factor, maxX: (x + 1) * factor - 1, minY: y * factor, maxY: (y + 1) * factor - 1 }; + } +} + +// the hierarchical descent: classify a tile into one of three outcomes — +// no edges cross + center inside -> whole subtree is inside: emit it arithmetically +// no edges cross + center outside -> whole subtree is outside: prune it +// edges cross -> the tile intersects: emit it, recurse into 4 children +function* descend(zoom: number, x: number, y: number, edges: Segment[], rings: Position[][], span: ZoomSpan): Generator { + const bounds = tileBounds(zoom, x, y); + // `edges` only holds edges that crossed the PARENT tile; re-filter for this tile, so the + // list keeps shrinking as the descent narrows onto the boundary + const clipped = edges.filter((edge) => segmentIntersectsBounds(edge, bounds)); + + if (clipped.length === 0) { + // the boundary does not pass through this tile, and a tile is a connected region — so it + // lies uniformly inside or uniformly outside; any single interior point decides which + const centerLon = midpoint(bounds.minLon, bounds.maxLon); + const centerLat = midpoint(bounds.minLat, bounds.maxLat); + // the center can never lie exactly on the boundary here (no edges cross this tile) + if (pointInRings(centerLon, centerLat, rings)) { + yield* subtreeRanges(zoom, x, y, span); + } + return; // center outside -> prune: no descendant can intersect + } + + if (zoom >= span.minZoom) { + // the boundary passes through this tile, and a polygon includes its boundary -> the tile + // itself intersects; emit it as a 1x1 range (only when inside the requested span) + yield { zoom, minX: x, maxX: x, minY: y, maxY: y }; + } + if (zoom < span.maxZoom) { + const childX = SUBTILES_PER_AXIS * x; // children occupy the doubled coordinates: + const childY = SUBTILES_PER_AXIS * y; // x' in {2x, 2x+1}, y' in {2y, 2y+1} + yield* descend(zoom + 1, childX, childY, clipped, rings, span); + yield* descend(zoom + 1, childX + 1, childY, clipped, rings, span); + yield* descend(zoom + 1, childX, childY + 1, clipped, rings, span); + yield* descend(zoom + 1, childX + 1, childY + 1, clipped, rings, span); + } +} + +/** + * computes the tile ranges intersecting the given footprint for every zoom level in the + * span, using hierarchical descent. Ranges of the same zoom are disjoint; their order is + * deterministic but otherwise unspecified. A tile touching the footprint with zero overlap + * area counts as intersecting. + * @param footprint footprint to cover, holes are supported + * @param span inclusive zoom span to emit ranges for + */ +export function* footprintToTileRanges(footprint: Footprint, span: ZoomSpan): Generator { + if (span.minZoom > span.maxZoom) { + throw new RangeError(`minZoom [${span.minZoom}] must not be greater than maxZoom [${span.maxZoom}]`); + } + if (span.minZoom < MIN_SUPPORTED_ZOOM || span.maxZoom > MAX_SUPPORTED_ZOOM) { + throw new RangeError(`zoom span [${span.minZoom}, ${span.maxZoom}] exceeds the supported range [${MIN_SUPPORTED_ZOOM}, ${MAX_SUPPORTED_ZOOM}]`); + } + + const rings = ringsOf(footprint); + const edges = ringsToEdges(rings); + // the geodetic grid has two root tiles (2x1) at zoom 0 — west (x=0) and east (x=1) + // hemispheres — so the descent starts from both; the walk can only visit real grid tiles + yield* descend(0, 0, 0, edges, rings, span); + yield* descend(0, 1, 0, edges, rings, span); +} + +export type { ZoomSpan }; diff --git a/src/geo/index.ts b/src/geo/index.ts index 52d566a..e2a447d 100644 --- a/src/geo/index.ts +++ b/src/geo/index.ts @@ -1,4 +1,5 @@ export * from './bboxUtils'; +export * from './footprintToTileRanges'; export * from './geoConvertor'; export * from './geoHash'; export * from './geoIntersection'; diff --git a/tests/unit/geo/footprintToTileRanges.spec.ts b/tests/unit/geo/footprintToTileRanges.spec.ts new file mode 100644 index 0000000..239195c --- /dev/null +++ b/tests/unit/geo/footprintToTileRanges.spec.ts @@ -0,0 +1,172 @@ +import type { Feature, MultiPolygon, Polygon } from 'geojson'; +import { bboxPolygon, booleanDisjoint, bbox as turfBbox } from '@turf/turf'; +import { ITileRange } from '../../../src/models/interfaces/geo/iTile'; +import { footprintToTileRanges } from '../../../src/geo/footprintToTileRanges'; +import { Footprint } from '../../../src/geo/geoIntersection'; +import { degreesPerTile } from '../../../src/geo/tiles'; + +// exhaustive ground truth: test every candidate tile in the footprint's bbox +// (padded by one tile so boundary-touching neighbors are candidates too) with the same +// inclusive intersection predicate the production algorithm must honor +const oracleTileSet = (footprint: Footprint, minZoom: number, maxZoom: number): Set => { + const tiles = new Set(); + const [minLon, minLat, maxLon, maxLat] = turfBbox(footprint); + for (let zoom = minZoom; zoom <= maxZoom; zoom++) { + const resolution = degreesPerTile(zoom); + const cols = 2 * Math.pow(2, zoom); + const rows = Math.pow(2, zoom); + const minX = Math.max(0, Math.floor((minLon + 180) / resolution) - 1); + const maxX = Math.min(cols - 1, Math.floor((maxLon + 180) / resolution) + 1); + const minY = Math.max(0, Math.floor((minLat + 90) / resolution) - 1); + const maxY = Math.min(rows - 1, Math.floor((maxLat + 90) / resolution) + 1); + for (let x = minX; x <= maxX; x++) { + for (let y = minY; y <= maxY; y++) { + const tilePolygon = bboxPolygon([x * resolution - 180, y * resolution - 90, (x + 1) * resolution - 180, (y + 1) * resolution - 90]); + if (!booleanDisjoint(tilePolygon, footprint)) { + tiles.add(`${zoom}/${x}/${y}`); + } + } + } + } + return tiles; +}; + +const expandToTileSet = (ranges: Iterable): Set => { + const tiles = new Set(); + for (const range of ranges) { + for (let x = range.minX; x <= range.maxX; x++) { + for (let y = range.minY; y <= range.maxY; y++) { + tiles.add(`${range.zoom}/${x}/${y}`); + } + } + } + return tiles; +}; + +const squarePolygon = (minLon: number, minLat: number, maxLon: number, maxLat: number): Polygon => ({ + type: 'Polygon', + coordinates: [ + [ + [minLon, minLat], + [maxLon, minLat], + [maxLon, maxLat], + [minLon, maxLat], + [minLon, minLat], + ], + ], +}); + +describe('footprintToTileRanges', () => { + describe('zoom span validation', () => { + it('should throw RangeError when minZoom is greater than maxZoom', () => { + const footprint = squarePolygon(0, 0, 90, 90); + + const generate = () => [...footprintToTileRanges(footprint, { minZoom: 5, maxZoom: 3 })]; + + expect(generate).toThrow(RangeError); + }); + + it.each([ + [-1, 5], + [0, 23], + ])('should throw RangeError when the span [%i, %i] exceeds the supported zoom levels', (minZoom, maxZoom) => { + const footprint = squarePolygon(0, 0, 90, 90); + + const generate = () => [...footprintToTileRanges(footprint, { minZoom, maxZoom })]; + + expect(generate).toThrow(RangeError); + }); + }); + + describe('tile coverage', () => { + it('should yield exactly the two root tiles for a world-covering footprint at zoom 0', () => { + const footprint = squarePolygon(-180, -90, 180, 90); + + const tiles = expandToTileSet(footprintToTileRanges(footprint, { minZoom: 0, maxZoom: 0 })); + + expect(tiles).toEqual(new Set(['0/0/0', '0/1/0'])); + }); + + it('should include tiles that only touch the footprint boundary with zero overlap area', () => { + // footprint aligned exactly to the zoom 2 grid (tile resolution 45 degrees): it fills + // tile (4,2) and touches all 8 neighbors along edges/corners without overlapping them + const footprint = squarePolygon(0, 0, 45, 45); + + const tiles = expandToTileSet(footprintToTileRanges(footprint, { minZoom: 2, maxZoom: 2 })); + + expect(tiles).toEqual(new Set(['2/3/1', '2/4/1', '2/5/1', '2/3/2', '2/4/2', '2/5/2', '2/3/3', '2/4/3', '2/5/3'])); + }); + }); + + describe('brute-force parity', () => { + const donut: Polygon = { + type: 'Polygon', + coordinates: [ + [ + [34.7, 31.2], + [35.3, 31.2], + [35.3, 31.9], + [34.7, 31.9], + [34.7, 31.2], + ], + [ + [34.9, 31.4], + [35.1, 31.4], + [35.1, 31.7], + [34.9, 31.7], + [34.9, 31.4], + ], + ], + }; + + const disjointMultiPolygon: MultiPolygon = { + type: 'MultiPolygon', + coordinates: [squarePolygon(10, 10, 10.5, 10.5).coordinates, squarePolygon(40, -20, 40.4, -19.6).coordinates], + }; + + const thinSliver: Polygon = { + type: 'Polygon', + coordinates: [ + [ + [34.7, 31.2], + [35.3, 31.9], + [35.3005, 31.9], + [34.7005, 31.2], + [34.7, 31.2], + ], + ], + }; + + const rectangleFeature: Feature = { + type: 'Feature', + properties: {}, + geometry: squarePolygon(34.72, 31.23, 35.28, 31.87), + }; + + it.each<[string, Footprint, number, number]>([ + ['off-grid rectangle', squarePolygon(34.72, 31.23, 35.28, 31.87), 0, 8], + ['donut with a hole', donut, 0, 10], + ['disjoint MultiPolygon parts', disjointMultiPolygon, 2, 8], + ['adversarial thin sliver', thinSliver, 4, 12], + ['Feature-wrapped rectangle', rectangleFeature, 3, 8], + ])('should match the exhaustive oracle key-by-key for a %s', (_name, footprint, minZoom, maxZoom) => { + const tiles = expandToTileSet(footprintToTileRanges(footprint, { minZoom, maxZoom })); + + expect(tiles).toEqual(oracleTileSet(footprint, minZoom, maxZoom)); + }); + + it('should emit disjoint ranges per zoom level', () => { + const ranges = [...footprintToTileRanges(donut, { minZoom: 0, maxZoom: 10 })]; + + const totalTiles = ranges.reduce((sum, range) => sum + (range.maxX - range.minX + 1) * (range.maxY - range.minY + 1), 0); + expect(totalTiles).toBe(expandToTileSet(ranges).size); + }); + + it('should be deterministic across runs', () => { + const first = [...footprintToTileRanges(thinSliver, { minZoom: 0, maxZoom: 10 })]; + const second = [...footprintToTileRanges(thinSliver, { minZoom: 0, maxZoom: 10 })]; + + expect(second).toEqual(first); + }); + }); +});