diff --git a/server/src/utils/country-borders.ts b/server/src/utils/country-borders.ts index 1332003e7..c8f2f5ed4 100644 --- a/server/src/utils/country-borders.ts +++ b/server/src/utils/country-borders.ts @@ -188,13 +188,64 @@ function pointInPackedRing( return inside; } +// Slack on the longitude sides of a part's box. The ray-cast compares the point's longitude +// with an intersection it computes in floating point, which can land a rounding step beyond +// the ring's own extremes; this is many orders of magnitude wider than that step and many +// narrower than a coordinate step of either bundle, so the box never rejects a point the +// ray-cast would have accepted. +const PART_BOX_SLACK = 1e-9; + +// [minLng, minLat, maxLng, maxLat] of every polygon's outer ring, four numbers per polygon, +// in the same divided-back doubles the ray-cast reads. Derived from the geometry on its first +// lookup and kept beside it (32 bytes per polygon), so neither bundle's packed layout changes. +const partBoxCache = new WeakMap(); + +function partBoxesOf(geom: PackedGeom, scale: number): Float64Array { + const cached = partBoxCache.get(geom); + if (cached) return cached; + const { coords, ringStarts, polyRingCounts } = geom; + const boxes = new Float64Array(polyRingCounts.length * 4); + let ri = 0; + for (const [k, rc] of polyRingCounts.entries()) { + let minLng = Infinity; + let minLat = Infinity; + let maxLng = -Infinity; + let maxLat = -Infinity; + for (let i = ringStarts[ri]; i < ringStarts[ri + 1]; i++) { + const x = coords[2 * i] / scale; + const y = coords[2 * i + 1] / scale; + if (x < minLng) minLng = x; + if (x > maxLng) maxLng = x; + if (y < minLat) minLat = y; + if (y > maxLat) maxLat = y; + } + boxes.set([minLng, minLat, maxLng, maxLat], k * 4); + ri += rc; + } + partBoxCache.set(geom, boxes); + return boxes; +} + // True when (lng,lat) falls inside a packed Polygon/MultiPolygon, honouring holes: the // first ring of each polygon is its outer boundary, the rest are holes. +// +// A polygon whose outer ring's box does not hold the point is skipped without walking its +// ring. A country is thousands of such polygons where it is mostly islands (Canada, Indonesia, +// Japan, Norway), and a point sits in the box of a handful at most, so the lookup walks a few +// rings rather than every one. The answer is the one the ray-cast alone gives: a ring has no +// edge crossing a latitude outside its own range, and none to the right of a point east of it. export function pointInPackedGeometry(lng: number, lat: number, geom: PackedGeom, scale: number): boolean { const { coords, ringStarts, polyRingCounts } = geom; + const boxes = partBoxesOf(geom, scale); let ri = 0; - for (const rc of polyRingCounts) { - if (pointInPackedRing(lng, lat, coords, scale, ringStarts[ri], ringStarts[ri + 1])) { + for (const [k, rc] of polyRingCounts.entries()) { + const b = k * 4; + const inBox = + lat >= boxes[b + 1] && + lat <= boxes[b + 3] && + lng >= boxes[b] - PART_BOX_SLACK && + lng <= boxes[b + 2] + PART_BOX_SLACK; + if (inBox && pointInPackedRing(lng, lat, coords, scale, ringStarts[ri], ringStarts[ri + 1])) { let inHole = false; for (let h = 1; h < rc; h++) { if (pointInPackedRing(lng, lat, coords, scale, ringStarts[ri + h], ringStarts[ri + h + 1])) { @@ -335,7 +386,27 @@ export function isPointInCountryBox(countryCode: string, lat: number, lng: numbe ); } +// Answers already worked out, by coordinate pair. A place keeps its coordinates between +// requests, and the places list and the trip bundle ask for every place of a trip on each +// read, so the same pairs come back all the time. The borders are fixed for the life of the +// process and the lookup reads nothing else, so an answer never goes stale. Oldest out first +// past the cap, which holds the map to a few megabytes. +const COUNTRY_ANSWERS_MAX = 20_000; +const countryAnswers = new Map(); + export function getCountryFromCoords(lat: number, lng: number): string | null { + const key = `${lat},${lng}`; + const known = countryAnswers.get(key); + if (known !== undefined) return known; + const answer = lookupCountryFromCoords(lat, lng); + if (countryAnswers.size >= COUNTRY_ANSWERS_MAX) { + countryAnswers.delete(countryAnswers.keys().next().value as string); + } + countryAnswers.set(key, answer); + return answer; +} + +function lookupCountryFromCoords(lat: number, lng: number): string | null { // Cheap prefilter: every country with a part-box containing the point. Keep the // area of the matching part so overlapping candidates can be ranked below. const candidates: { code: string; area: number }[] = []; diff --git a/server/tests/unit/utils/country-borders-lookup.test.ts b/server/tests/unit/utils/country-borders-lookup.test.ts new file mode 100644 index 000000000..e939ffcc3 --- /dev/null +++ b/server/tests/unit/utils/country-borders-lookup.test.ts @@ -0,0 +1,98 @@ +import { getCountryFromCoords, pointInPackedGeometry, type PackedGeom } from '../../../src/utils/country-borders'; + +import { describe, expect, it } from 'vitest'; + +const SCALE = 100; + +/** Packs polygons given as rings of [lng, lat] the way both bundles are packed. */ +function pack(polygons: [number, number][][][]): PackedGeom { + const rings = polygons.flat(); + const coords = new Int16Array(rings.reduce((n, ring) => n + ring.length, 0) * 2); + const ringStarts = new Uint32Array(rings.length + 1); + let p = 0; + for (const [r, ring] of rings.entries()) { + ringStarts[r] = p; + for (const [lng, lat] of ring) { + coords[2 * p] = Math.round(lng * SCALE); + coords[2 * p + 1] = Math.round(lat * SCALE); + p++; + } + } + ringStarts[rings.length] = p; + return { coords, ringStarts, polyRingCounts: Uint32Array.from(polygons.map((polygon) => polygon.length)) }; +} + +const square = (minLng: number, minLat: number, size: number): [number, number][] => [ + [minLng, minLat], + [minLng + size, minLat], + [minLng + size, minLat + size], + [minLng, minLat + size], + [minLng, minLat], +]; + +describe('pointInPackedGeometry', () => { + // An island at the origin, and a second one far east of it with a lake in the middle. + const geom = pack([[square(0, 0, 2)], [square(10, 10, 6), square(12, 12, 2)]]); + const inside = (lng: number, lat: number) => pointInPackedGeometry(lng, lat, geom, SCALE); + + it('BORDERS-001: finds a point in any polygon of the geometry, not only the first', () => { + expect(inside(1, 1)).toBe(true); + expect(inside(11, 11)).toBe(true); + expect(inside(15.5, 15.5)).toBe(true); + }); + + it('BORDERS-002: a point in a hole is outside, a point between the polygons too', () => { + expect(inside(13, 13)).toBe(false); + expect(inside(5, 5)).toBe(false); + expect(inside(1, 11)).toBe(false); + expect(inside(11, 1)).toBe(false); + }); + + it('BORDERS-003: answers a point on or beside a box edge as the ray-cast alone does', () => { + // The even-odd rule counts the west and south edges in and the east and north ones out. + expect(inside(0, 1)).toBe(true); + expect(inside(1, 0)).toBe(true); + expect(inside(2, 1)).toBe(false); + expect(inside(1, 2)).toBe(false); + expect(inside(-1e-7, 1)).toBe(false); + expect(inside(2 - 1e-7, 1)).toBe(true); + expect(inside(16 + 1e-7, 11)).toBe(false); + expect(inside(10, 16 - 1e-7)).toBe(true); + }); + + it('BORDERS-004: gives the same answers on a second pass over the same geometry', () => { + const points: [number, number][] = [ + [1, 1], + [13, 13], + [15, 11], + [9.99, 11], + [30, 30], + ]; + const first = points.map(([lng, lat]) => inside(lng, lat)); + expect(points.map(([lng, lat]) => inside(lng, lat))).toEqual(first); + expect(first).toEqual([true, false, true, false, false]); + }); +}); + +describe('getCountryFromCoords', () => { + it('BORDERS-005: answers a pair it has seen before the same way, a pair in no country included', () => { + const tokyo = getCountryFromCoords(35.68, 139.76); + const midPacific = getCountryFromCoords(-48.5, -128.5); + + expect(tokyo).toBe('JP'); + expect(midPacific).toBeNull(); + expect(getCountryFromCoords(35.68, 139.76)).toBe(tokyo); + expect(getCountryFromCoords(-48.5, -128.5)).toBeNull(); + // Latitude and longitude are not interchangeable in the remembered pair. + expect(getCountryFromCoords(139.76, 35.68)).not.toBe('JP'); + }); + + it('BORDERS-006: still answers right after more pairs than it remembers have passed through', () => { + expect(getCountryFromCoords(48.86, 2.35)).toBe('FR'); + // Open ocean south-east of Pitcairn: in no country's box, so each lookup is cheap. + for (let i = 0; i < 20_500; i++) getCountryFromCoords(-50 + (i % 100) * 0.01, -130 + Math.floor(i / 100) * 0.01); + + expect(getCountryFromCoords(48.86, 2.35)).toBe('FR'); + expect(getCountryFromCoords(-50, -130)).toBeNull(); + }); +});