The algorithm

Algorithm

Seven steps. One groundbreaking algorithm.

Step 1 of 7 · Get the buildings

Microsoft's global building map was built on an AI model trained with satellite imaging of the whole world. A scan is about 4 by 4 kilometers and we only keep the buildings whose center is inside the box.
Building outlines load one square at a time and get cut to fit the box.

The code

detectVillages.mjs

Here's the detector code every scan runs (commented) (detect-villages.js on the command line calls it). Open a step above to jump to that part, or download the whole file to use!

node detect-villages.js --south 9.02 --west 38.70 --north 9.06 --east 38.74

detectVillages.mjs2,259 lines · comments includedLines 1,952–1,987 · Get the buildings

1// Village polygons from satellite tiles, scored against OpenStreetMap.

2// No model: color and shape only.

3//

4// This is the detector itself, used two ways:

5// - from the command line: node detect-villages.js --south 9.02 --west 38.70 --north 9.06 --east 38.74

6// (detect-villages.js at the repo root is a thin wrapper around this file)

7// - inside the web app: the "Scan this view" API route imports detectVillages() directly,

8// so it also runs on Vercel, where there is no separate script to start.

9//

10// Pipeline:

11// 1. Get the buildings in the box: Microsoft's building footprints (an outline per

12// building, from an AI model run over satellite imagery). Only where Microsoft has

13// no data at all, fall back to finding roofs by color in Esri satellite tiles.

14// 2. Join buildings within ~50 m into settlements and trace each outline.

15// 3. If an OSM place sits inside or near a settlement, certainty is 100.

16// Otherwise score it from its buildings and shape, capped at 99.

17// 4. Estimate population (homes × people per home) and the money to vaccinate it.

18// 5. Price a route from the nearest clinic along OSM roads.

19// 6. Recommend the polygon with the lowest predicted cost, and that route.

20// 7. Plan one team tour: clinic → every settlement → clinic, ordered for lowest cost.

21//

22// The function returns GeoJSON. Progress lines go to stderr.

24import jpeg from "jpeg-js";

25import { existsSync, readdirSync, readFileSync } from "node:fs";

26import { mkdir, readFile, rename, stat, writeFile } from "node:fs/promises";

27import { tmpdir } from "node:os";

28import { join } from "node:path";

29import { gunzipSync } from "node:zlib";

31// Zoom 17 is about 1.2 m per pixel, small enough to see individual roofs.

32const ZOOM = 17;

33// Esri World Imagery tiles are 256×256 JPEG images.

34const TILE = 256;

35// Erode-then-dilate radius. Knocks out 1-pixel speckle before villages are joined.

36const OPEN_RADIUS = 1;

37// Dilate-then-erode radius. Bridges gaps of about 2×12 pixels (~30 m) between roofs.

38const CLOSE_RADIUS = 12;

39// Blobs smaller than this are a shed or noise, not a village. ~400 px is ~500 m².

40const MIN_PIXELS = 400;

41// --- Building footprints (Microsoft Global ML Building Footprints, ODbL) ---

42// List of every footprint file, one per zoom-9 map square (~78 km), with its download link.

43const FOOTPRINT_LINKS_URL = "https://minedbuildings.z5.web.core.windows.net/global-buildings/dataset-links.csv";

44const FOOTPRINT_QUADKEY_ZOOM = 9;

45// Downloaded files are kept here and reused; the link list is refreshed after a week.

46const FOOTPRINT_CACHE_DIR = join(tmpdir(), "hestia-footprints");

47const FOOTPRINT_LINKS_MAX_AGE_MS = 7 * 24 * 3600 * 1000;

48// Footprints are drawn on a grid this fine (meters per cell) to join and trace them.

49const FOOTPRINT_RES_M = 3;

50// Buildings closer than this belong to the same settlement.

51const FOOTPRINT_LINK_M = 50;

52// Smaller structures are granaries, sheds or latrines: part of a settlement, but nobody's home.

53const MIN_HOME_M2 = 10;

54// A building bigger than this counts as one home per this much floor area (compounds, blocks).

55const HOME_M2 = 150;

56// Average people per household in Nigeria (Nigeria DHS 2018).

57const PEOPLE_PER_HOME = 4.7;

58// The share of a settlement's outline covered by buildings that earns full marks (footprints).

59const FULL_COVERAGE = 0.15;

61// A blob with no separate roof at all is bare ground or a road, not a settlement.

62const MIN_ROOF_CLUSTERS = 1;

63// Red pixels count as roof only on a sharp edge. Red laterite soil and dirt roads are

64// as red as tile roofs, but they fade gradually; a roof has a crisp outline.

65const RED_ROOF_MIN_EDGE = 22;

66// Painted or galvanized metal roofs often read teal at zoom 17: green and blue both well

67// above red and close to each other. Houses are only 4–6 px across, so they are all edge.

68const METAL_ROOF_MIN_EDGE = 15;

69// A roof cluster bigger than this is a patch of ground, not a house. ~360 px is ~500 m².

70const MAX_ROOF_CLUSTER_PX = 360;

71// Houses are rectangles, so a cluster should fill a fair share of its bounding box.

72// Soil fragments are ragged and fill little of theirs.

73const MIN_ROOF_FILL = 0.4;

74// Average width of a cluster (pixels ÷ longest side). Crop-row lines and the bright rim

75// around a dark tree are 1–2 px wide; a house is 4–10 px across. 3 px is ~3.5 m.

76const MIN_ROOF_THICKNESS_PX = 3;

77// A house is at most a few times longer than it is wide. Longer is a field edge or track.

78const MAX_ROOF_ASPECT = 4;

79// Below this compactness a blob is a stringy web of soil and tracks, not a village.

80// An OSM place still confirms it, but then people come from roofs only, not its size.

81const MIN_COMPACTNESS = 0.06;

82// An OSM place this close to a polygon counts as the same settlement.

83const OSM_METERS = 150;

84// Douglas–Peucker tolerance in pixels. Drops stair-steps along a straight wall.

85const SIMPLIFY_PX = 2;

86// A safety cap so a huge bbox cannot download thousands of tiles.

87const MAX_TILES = 400;

88// Public satellite imagery. No API key. URL shape is .../tile/{z}/{y}/{x}.

89const ESRI =

90 "https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile";

91// Public OSM query services. The second is used when the first is busy.

92const OVERPASS_SERVERS = [

93 // The main server. Its /api/status page says when this connection's next slot frees up.

94 { url: "https://overpass-api.de/api", status: true },

95 // Public mirrors, each with its own limits. Any of them can be down; the client moves on.

96 { url: "https://overpass.kumi.systems/api", status: false },

97 { url: "https://overpass.private.coffee/api", status: false },

98 { url: "https://maps.mail.ru/osm/tools/overpass/api", status: false },

99 { url: "https://overpass.osm.jp/api", status: false },

100];

101// Overpass rejects generic clients (curl's default gets 406), so always say who we are.

102const OVERPASS_USER_AGENT = "village-polygon-detector/1.0";

103// A query is given this long before the server is treated as stuck. Road queries can be big.

104const OVERPASS_REQUEST_MS = 200 * 1000;

105// Total tries across all servers before a query gives up.

106const OVERPASS_ATTEMPTS = 10;

107// Longest a server rests after failing before it's tried again.

108const OVERPASS_MAX_COOLDOWN_MS = 5 * 60 * 1000;

109// How far outside the village box to look for a clinic, widening until one is found.

110// 0.36° is about 40 km; remote villages are often 20–50 km from the nearest facility.

111// The last step (~660 km) covers any country, so every scan gets a clinic.

112const CLINIC_SEARCH_PADS_DEG = [0.36, 0.72, 1.5, 3, 6];

113// The search width used when even the widest step finds nothing (no clinic mapped at all).

114const CLINIC_PAD_DEG = CLINIC_SEARCH_PADS_DEG[CLINIC_SEARCH_PADS_DEG.length - 1];

115// All drivable roads are downloaded this far around the box. Farther out, toward a distant

116// clinic, only main roads are, which keeps the download small and is how a team would drive.

117const LOCAL_ROADS_PAD_DEG = 0.05;

118// Clinics tried per village. The cheapest road trip of these three wins.

119const CLINICS_PER_VILLAGE = 3;

120// A road farther than this is not a usable approach. The trip is priced off-road.

121const MAX_SNAP_METERS = 8000;

122// One roof cluster stands in for one dwelling.

123const PEOPLE_PER_BUILDING = 5;

124// Used when the roof test missed houses: one dwelling per this much settlement area.

125const METERS_PER_DWELLING = 250;

126// Vaccine, syringe, and a share of the health worker's time, per person.

127const HEALTHCARE_USD_PER_PERSON = 10;

128// Cooler, records, and the stop itself, paid once even for a tiny settlement.

129const HEALTHCARE_USD_PER_STOP = 100;

130// Last mile with no road. Slow, so the team's time dominates the fuel.

131const OFFROAD_USD_PER_KM = 8;

132// The team tour skips detections below this certainty: they are often bare ground, not homes.

133const TOUR_MIN_CERTAINTY = 30;

134// 2-opt passes over the tour. Each pass tries every pair of legs; it usually settles in a few.

135const TOUR_MAX_PASSES = 50;

136

137// Check a { south, west, north, east } box. Latitudes must stay inside Web Mercator's usable range.

138function validateBbox(bbox) {

139 const keys = ["south", "west", "north", "east"];

140 if (keys.some((key) => !Number.isFinite(bbox?.[key]))) {

141 throw new Error("bbox needs numeric south, west, north and east");

142 }

143 if (!(bbox.south < bbox.north) || !(bbox.west < bbox.east)) {

144 throw new Error("bbox must have south < north and west < east");

145 }

146 if (bbox.south <= -85 || bbox.north >= 85) {

147 throw new Error("latitude must stay inside -85 and 85");

148 }

149}

150

151// Slippy-map column. Longitude maps linearly across the world.

152function lonToTileX(lon, zoom) {

153 return Math.floor(((lon + 180) / 360) * 2 ** zoom);

154}

155

156// Slippy-map row. Latitude is stretched with the Mercator formula.

157// Row 0 is the north edge of the world, so a north latitude has a smaller y.

158function latToTileY(lat, zoom) {

159 const rad = (lat * Math.PI) / 180;

160 const merc = Math.log(Math.tan(rad) + 1 / Math.cos(rad));

161 return Math.floor(((1 - merc / Math.PI) / 2) * 2 ** zoom);

162}

163

164// Turn a pixel corner in the stitched image back into longitude and latitude.

165// px/py are measured from the top-left of the mosaic, not from a single tile.

166function cornerToLonLat(px, py, grid) {

167 const world = TILE * 2 ** grid.zoom;

168 const gx = grid.x0 * TILE + px;

169 const gy = grid.y0 * TILE + py;

170 const lon = (gx / world) * 360 - 180;

171 const lat = (Math.atan(Math.sinh(Math.PI * (1 - (2 * gy) / world))) * 180) / Math.PI;

172 return [lon, lat];

173}

174

175// Inverse of cornerToLonLat: where a longitude falls inside the mosaic.

176function lonToPixelX(lon, grid) {

177 const world = TILE * 2 ** grid.zoom;

178 return ((lon + 180) / 360) * world - grid.x0 * TILE;

179}

180

181// Inverse for latitude. Used to throw away tile pixels that sit outside the bbox.

182function latToPixelY(lat, grid) {

183 const world = TILE * 2 ** grid.zoom;

184 const rad = (lat * Math.PI) / 180;

185 const merc = Math.log(Math.tan(rad) + 1 / Math.cos(rad));

186 return ((1 - merc / Math.PI) / 2) * world - grid.y0 * TILE;

187}

188

189// Tiles cover a rectangle slightly larger than the requested bbox.

190// Zero every mask pixel outside the bbox so a village just past the edge is ignored.

191function clipToBbox(mask, width, height, left, top, right, bottom) {

192 const x0 = Math.max(0, Math.floor(left));

193 const y0 = Math.max(0, Math.floor(top));

194 const x1 = Math.min(width, Math.ceil(right));

195 const y1 = Math.min(height, Math.ceil(bottom));

196 for (let y = 0; y < height; y++) {

197 for (let x = 0; x < width; x++) {

198 if (x < x0 || x >= x1 || y < y0 || y >= y1) mask[y * width + x] = 0;

199 }

200 }

201}

202

203// Download one JPEG tile and decode it to RGBA bytes. Retry a flaky response.

204async function fetchTile(zoom, x, y) {

205 const url = `${ESRI}/${zoom}/${y}/${x}`;

206 let lastError;

207 for (let attempt = 0; attempt < 3; attempt++) {

208 try {

209 const response = await fetch(url, {

210 headers: { "User-Agent": "village-polygon-detector/1.0" },

211 });

212 if (!response.ok) throw new Error(`${response.status} for ${url}`);

213 const buffer = Buffer.from(await response.arrayBuffer());

214 // JPEG files start with the SOI marker FF D8. Anything else is an error page.

215 if (buffer[0] !== 0xff || buffer[1] !== 0xd8) {

216 throw new Error(`tile ${zoom}/${y}/${x} was not a JPEG`);

217 }

218 const image = jpeg.decode(buffer, { useTArray: true, formatAsRGBA: true });

219 if (image.width !== TILE || image.height !== TILE) {

220 throw new Error(`unexpected tile size ${image.width}x${image.height}`);

221 }

222 return image.data;

223 } catch (error) {

224 lastError = error;

225 }

226 }

227 throw lastError;

228}

229

230// Copy one decoded tile into the big RGB mosaic. Alpha is dropped.

231// dx/dy are the tile's top-left pixel inside that mosaic.

232function blit(rgb, width, pixels, dx, dy) {

233 for (let row = 0; row < TILE; row++) {

234 const src = row * TILE * 4;

235 const dst = ((dy + row) * width + dx) * 3;

236 for (let col = 0; col < TILE; col++) {

237 const from = src + col * 4;

238 const to = dst + col * 3;

239 rgb[to] = pixels[from];

240 rgb[to + 1] = pixels[from + 1];

241 rgb[to + 2] = pixels[from + 2];

242 }

243 }

244}

245

246// Run `fn` over `items` with at most `limit` downloads in flight.

247// JavaScript is single-threaded, so writing each tile into its own region is safe.

248async function mapPool(items, limit, fn) {

249 let next = 0;

250 async function worker() {

251 while (next < items.length) {

252 const index = next++;

253 await fn(items[index]);

254 }

255 }

256 const workers = Math.min(limit, items.length);

257 await Promise.all(Array.from({ length: workers }, () => worker()));

258}

259

260// Build one RGB image covering the bbox, plus the tile origin needed for lat/lng.

261async function loadSatellite(bbox) {

262 // Tile y grows southward, so the north edge is the smaller row index.

263 const x0 = lonToTileX(bbox.west, ZOOM);

264 const x1 = lonToTileX(bbox.east, ZOOM);

265 const y0 = latToTileY(bbox.north, ZOOM);

266 const y1 = latToTileY(bbox.south, ZOOM);

267 const tilesX = x1 - x0 + 1;

268 const tilesY = y1 - y0 + 1;

269 if (tilesX < 1 || tilesY < 1) throw new Error("bbox does not cover a tile");

270 if (tilesX * tilesY > MAX_TILES) {

271 throw new Error(

272 `bbox covers ${tilesX * tilesY} tiles; keep it under ${MAX_TILES}`

273 );

274 }

275

276 const width = tilesX * TILE;

277 const height = tilesY * TILE;

278 const rgb = new Uint8Array(width * height * 3);

279 const tiles = [];

280 for (let y = y0; y <= y1; y++) {

281 for (let x = x0; x <= x1; x++) tiles.push({ x, y });

282 }

283

284 let done = 0;

285 console.error(`loading ${tiles.length} satellite tiles at zoom ${ZOOM}`);

286 await mapPool(tiles, 6, async (tile) => {

287 const pixels = await fetchTile(ZOOM, tile.x, tile.y);

288 blit(rgb, width, pixels, (tile.x - x0) * TILE, (tile.y - y0) * TILE);

289 done += 1;

290 if (done === tiles.length || done % 25 === 0) {

291 console.error(`tiles ${done}/${tiles.length}`);

292 }

293 });

294

295 return { rgb, width, height, zoom: ZOOM, x0, y0 };

296}

297

298// Brightness difference between this pixel and its four neighbors.

299// Roof edges and walls score high. A flat field scores low.

300function contrastAt(rgb, width, height, x, y) {

301 const index = (y * width + x) * 3;

302 const center = (rgb[index] + rgb[index + 1] + rgb[index + 2]) / 3;

303 let maxDiff = 0;

304 const neighbors = [

305 [1, 0],

306 [-1, 0],

307 [0, 1],

308 [0, -1],

309 ];

310 for (const [dx, dy] of neighbors) {

311 const nx = x + dx;

312 const ny = y + dy;

313 if (nx < 0 || ny < 0 || nx >= width || ny >= height) continue;

314 const n = (ny * width + nx) * 3;

315 const brightness = (rgb[n] + rgb[n + 1] + rgb[n + 2]) / 3;

316 const diff = Math.abs(center - brightness);

317 if (diff > maxDiff) maxDiff = diff;

318 }

319 return maxDiff;

320}

321

322// Green dominates both other channels: trees, crops, grass.

323function isVegetation(r, g, b) {

324 return g > 40 && g > r + 12 && g > b + 12;

325}

326

327// Dark and blue: open water. Bright blue metal roofs are handled separately.

328function isWater(r, g, b) {

329 const brightness = (r + g + b) / 3;

330 return brightness < 90 && b > 40 && b > r + 12 && b >= g;

331}

332

333// Shadow, dense forest, or bare dark soil. Average brightness under 40.

334function isDark(r, g, b) {

335 return r + g + b < 120;

336}

337

338// Roof colors, checked only after vegetation, water, and darkness are rejected.

339// Red tile or flat gray that also has an edge, pale concrete, blue corrugated metal,

340// or teal-looking painted metal on an edge.

341// Gray without an edge is asphalt, so it does not count as a roof by itself.

342function isRoof(r, g, b, edge) {

343 const max = Math.max(r, g, b);

344 const min = Math.min(r, g, b);

345 const brightness = (r + g + b) / 3;

346 if (r > 90 && r > g + 25 && r > b + 25 && edge >= RED_ROOF_MIN_EDGE) return true;

347 if (min > 165 && max - min < 40) return true;

348 if (b > 90 && b > r + 18 && b > g + 8 && brightness >= 100) return true;

349 if (g > r + 25 && b > r + 25 && Math.abs(g - b) <= 15 && brightness >= 60 && edge >= METAL_ROOF_MIN_EDGE) {

350 return true;

351 }

352 return brightness > 90 && brightness < 190 && max - min < 22 && edge >= 18;

353}

354

355// Two masks from the same pass:

356// roof : pixels that look like a roof, used later to count buildings

357// built : roof pixels plus any high-contrast structure (walls, roof edges)

358function classify(rgb, width, height) {

359 const built = new Uint8Array(width * height);

360 const roof = new Uint8Array(width * height);

361 for (let y = 0; y < height; y++) {

362 for (let x = 0; x < width; x++) {

363 const pixel = y * width + x;

364 const offset = pixel * 3;

365 const r = rgb[offset];

366 const g = rgb[offset + 1];

367 const b = rgb[offset + 2];

368 if (isVegetation(r, g, b) || isWater(r, g, b) || isDark(r, g, b)) continue;

369 const edge = contrastAt(rgb, width, height, x, y);

370 const roofLike = isRoof(r, g, b, edge);

371 if (roofLike) roof[pixel] = 1;

372 if (roofLike || edge >= 28) built[pixel] = 1;

373 }

374 }

375 return { built, roof };

376}

377

378// One-row sliding-window max (dilate) or min (erode).

379// The deque keeps indexes whose values can still win the current window,

380// so a whole row costs linear time instead of window-width times the row.

381function filterH(src, width, height, radius, dst, dilate) {

382 const deque = new Int32Array(width);

383 for (let y = 0; y < height; y++) {

384 const row = y * width;

385 let head = 0;

386 let tail = 0;

387 for (let x = 0; x < width + radius; x++) {

388 if (x < width) {

389 const value = src[row + x];

390 while (tail > head) {

391 const previous = src[row + deque[tail - 1]];

392 if (dilate ? previous > value : previous < value) break;

393 tail -= 1;

394 }

395 deque[tail++] = x;

396 }

397 const outX = x - radius;

398 if (outX < 0) continue;

399 const left = outX - radius;

400 while (tail > head && deque[head] < left) head += 1;

401 dst[row + outX] = src[row + deque[head]];

402 }

403 }

404}

405

406// Same window, applied down each column. After filterH, the pair is a square.

407function filterV(src, width, height, radius, dst, dilate) {

408 const deque = new Int32Array(height);

409 for (let x = 0; x < width; x++) {

410 let head = 0;

411 let tail = 0;

412 for (let y = 0; y < height + radius; y++) {

413 if (y < height) {

414 const value = src[y * width + x];

415 while (tail > head) {

416 const previous = src[deque[tail - 1] * width + x];

417 if (dilate ? previous > value : previous < value) break;

418 tail -= 1;

419 }

420 deque[tail++] = y;

421 }

422 const outY = y - radius;

423 if (outY < 0) continue;

424 const top = outY - radius;

425 while (tail > head && deque[head] < top) head += 1;

426 dst[outY * width + x] = src[deque[head] * width + x];

427 }

428 }

429}

430

431// Square dilate (grow foreground) or erode (shrink foreground).

432// `scratch` holds the horizontal pass so src and dst can stay separate.

433function rectFilter(src, width, height, radius, dilate, scratch, dst) {

434 filterH(src, width, height, radius, scratch, dilate);

435 filterV(scratch, width, height, radius, dst, dilate);

436}

437

438// Turn scattered building pixels into settlement blobs.

439// Open (erode, then dilate) deletes speckle. Close (dilate, then erode)

440// fills the gaps between nearby roofs without growing the outer footprint.

441function settlementMask(built, width, height) {

442 const count = width * height;

443 const scratch = new Uint8Array(count);

444 const a = new Uint8Array(count);

445 const b = new Uint8Array(count);

446 rectFilter(built, width, height, OPEN_RADIUS, false, scratch, a);

447 rectFilter(a, width, height, OPEN_RADIUS, true, scratch, b);

448 rectFilter(b, width, height, CLOSE_RADIUS, true, scratch, a);

449 rectFilter(a, width, height, CLOSE_RADIUS, false, scratch, b);

450 return b;

451}

452

453// How many separate buildings sit in this blob, and how many pixels those buildings cover.

454// Roof pixels were stored before the close, so a gap between houses still splits them.

455// A cluster under 8 pixels is glare or a speck; one too big or too ragged is ground;

456// one too thin or too long is a crop-row line, a track, or the rim of a tree.

457function countRoofClusters(indices, width, height) {

458 const remaining = new Set(indices);

459 let clusters = 0;

460 let pixels = 0;

461 const stack = [];

462 for (const start of indices) {

463 if (!remaining.has(start)) continue;

464 let size = 0;

465 let minX = Infinity;

466 let minY = Infinity;

467 let maxX = -Infinity;

468 let maxY = -Infinity;

469 remaining.delete(start);

470 stack.push(start);

471 while (stack.length) {

472 const pixel = stack.pop();

473 size += 1;

474 const x = pixel % width;

475 const y = (pixel - x) / width;

476 if (x < minX) minX = x;

477 if (x > maxX) maxX = x;

478 if (y < minY) minY = y;

479 if (y > maxY) maxY = y;

480 const neighbors = [

481 [x - 1, y],

482 [x + 1, y],

483 [x, y - 1],

484 [x, y + 1],

485 ];

486 for (const [nx, ny] of neighbors) {

487 if (nx < 0 || ny < 0 || nx >= width || ny >= height) continue;

488 const next = ny * width + nx;

489 if (!remaining.has(next)) continue;

490 remaining.delete(next);

491 stack.push(next);

492 }

493 }

494 const boxW = maxX - minX + 1;

495 const boxH = maxY - minY + 1;

496 const fill = size / (boxW * boxH);

497 const thickness = size / Math.max(boxW, boxH);

498 const aspect = Math.max(boxW, boxH) / Math.min(boxW, boxH);

499 if (

500 size >= 8 &&

501 size <= MAX_ROOF_CLUSTER_PX &&

502 fill >= MIN_ROOF_FILL &&

503 thickness >= MIN_ROOF_THICKNESS_PX &&

504 aspect <= MAX_ROOF_ASPECT

505 ) {

506 clusters += 1;

507 pixels += size;

508 }

509 }

510 return { clusters, pixels };

511}

512

513// Stitch directed pixel edges into loops and keep the largest one.

514// That loop is the outer contour. Smaller loops are courtyards and holes.

515// `edges` is a flat list of x1,y1,x2,y2 segments in pixel-corner coordinates.

516function outerRing(height, edges) {

517 const stride = height + 1;

518 const nexts = new Map();

519 for (let i = 0; i < edges.length; i += 4) {

520 const x1 = edges[i];

521 const y1 = edges[i + 1];

522 const x2 = edges[i + 2];

523 const y2 = edges[i + 3];

524 const key = x1 * stride + y1;

525 let list = nexts.get(key);

526 if (!list) {

527 list = [];

528 nexts.set(key, list);

529 }

530 list.push(x2 * stride + y2);

531 }

532

533 let best = null;

534 let bestArea = 0;

535 for (const [startKey, list] of nexts) {

536 while (list.length) {

537 const ring = [];

538 let cursor = startKey;

539 const guard = edges.length;

540 for (let step = 0; step < guard; step++) {

541 const x = Math.floor(cursor / stride);

542 const y = cursor % stride;

543 ring.push([x, y]);

544 const outs = nexts.get(cursor);

545 if (!outs || outs.length === 0) break;

546 const next = outs.pop();

547 // Back at the start: the loop is closed. Do not append the start twice.

548 if (next === startKey) break;

549 cursor = next;

550 }

551 if (ring.length < 3) continue;

552 // Shoelace area. The biggest loop is the outside of the blob.

553 let area2 = 0;

554 for (let i = 0; i < ring.length; i++) {

555 const a = ring[i];

556 const b = ring[(i + 1) % ring.length];

557 area2 += a[0] * b[1] - b[0] * a[1];

558 }

559 const area = Math.abs(area2) / 2;

560 if (area > bestArea) {

561 bestArea = area;

562 best = ring;

563 }

564 }

565 }

566 if (!best || bestArea < 1) return null;

567 return best;

568}

569

570// Label each connected settlement and trace its outline.

571// Neighbors are 4-connected (up, down, left, right), matching the edge rules below.

572// A foreground pixel emits an edge on any side whose neighbor is background.

573// `labels`, when given, gets each blob's label at every one of its pixels, so callers can

574// tell which blob a point falls in. The blob's own `label` field says which value it used.

575function components(mask, roof, width, height, { minPixels = MIN_PIXELS, labels = null } = {}) {

576 const seen = new Uint8Array(mask.length);

577 const found = [];

578 const stack = [];

579 let label = 0;

580

581 for (let index = 0; index < mask.length; index++) {

582 if (!mask[index] || seen[index]) continue;

583 seen[index] = 1;

584 stack.push(index);

585 label += 1;

586 let pixels = 0;

587 const roofIndices = [];

588 const edges = [];

589

590 while (stack.length) {

591 const pixel = stack.pop();

592 pixels += 1;

593 if (labels) labels[pixel] = label;

594 const x = pixel % width;

595 const y = (pixel - x) / width;

596 if (roof[pixel]) roofIndices.push(pixel);

597

598 const left = x > 0 && mask[pixel - 1];

599 const right = x + 1 < width && mask[pixel + 1];

600 const up = y > 0 && mask[pixel - width];

601 const down = y + 1 < height && mask[pixel + width];

602 // Edges walk clockwise around the pixel: top, right, bottom, left.

603 if (!up) edges.push(x, y, x + 1, y);

604 if (!right) edges.push(x + 1, y, x + 1, y + 1);

605 if (!down) edges.push(x + 1, y + 1, x, y + 1);

606 if (!left) edges.push(x, y + 1, x, y);

607

608 if (left && !seen[pixel - 1]) {

609 seen[pixel - 1] = 1;

610 stack.push(pixel - 1);

611 }

612 if (right && !seen[pixel + 1]) {

613 seen[pixel + 1] = 1;

614 stack.push(pixel + 1);

615 }

616 if (up && !seen[pixel - width]) {

617 seen[pixel - width] = 1;

618 stack.push(pixel - width);

619 }

620 if (down && !seen[pixel + width]) {

621 seen[pixel + width] = 1;

622 stack.push(pixel + width);

623 }

624 }

625

626 if (pixels < minPixels) continue;

627 const ring = outerRing(height, edges);

628 if (!ring) continue;

629 // Only pixels in house-like clusters count toward the roof share.

630 const roofs = countRoofClusters(roofIndices, width, height);

631 found.push({

632 ring,

633 pixels,

634 label,

635 roofPixels: roofs.pixels,

636 buildings: roofs.clusters,

637 });

638 }

639 return found;

640}

641

642// Distance from a point to the infinite line through start and end.

643// Used by the simplifier to decide which corners are worth keeping.

644function perpendicularDistance(point, start, end) {

645 const dx = end[0] - start[0];

646 const dy = end[1] - start[1];

647 const length = Math.hypot(dx, dy);

648 if (length === 0) return Math.hypot(point[0] - start[0], point[1] - start[1]);

649 return Math.abs((point[0] - start[0]) * dy - (point[1] - start[1]) * dx) / length;

650}

651

652// Ramer–Douglas–Peucker. Keeps a point only when it bends the outline

653// by more than `epsilon` pixels. Straight runs collapse to their endpoints.

654function simplify(points, epsilon) {

655 const count = points.length;

656 if (count <= 2) return points.slice();

657 const keep = new Uint8Array(count);

658 keep[0] = 1;

659 keep[count - 1] = 1;

660 const stack = [[0, count - 1]];

661 while (stack.length) {

662 const [start, end] = stack.pop();

663 let maxDistance = 0;

664 let index = -1;

665 for (let i = start + 1; i < end; i++) {

666 const distance = perpendicularDistance(points[i], points[start], points[end]);

667 if (distance > maxDistance) {

668 maxDistance = distance;

669 index = i;

670 }

671 }

672 if (index !== -1 && maxDistance > epsilon) {

673 keep[index] = 1;

674 stack.push([start, index], [index, end]);

675 }

676 }

677 const simplified = [];

678 for (let i = 0; i < count; i++) if (keep[i]) simplified.push(points[i]);

679 return simplified;

680}

681

682// Local flat map in meters around latitude lat0.

683// 111320 m per degree of longitude at the equator, 110540 m per degree of latitude.

684function toMeters(lon, lat, lat0) {

685 const cos = Math.cos((lat0 * Math.PI) / 180);

686 return [lon * 111320 * cos, lat * 110540];

687}

688

689// Signed area (for winding order) and compactness, both in meters.

690// Compactness is 1 for a circle and near 0 for a long road or field edge.

691function ringStats(ring) {

692 const lat0 = ring.reduce((sum, point) => sum + point[1], 0) / ring.length;

693 const meters = ring.map(([lon, lat]) => toMeters(lon, lat, lat0));

694 let area2 = 0;

695 let perimeter = 0;

696 for (let i = 0; i < meters.length; i++) {

697 const a = meters[i];

698 const b = meters[(i + 1) % meters.length];

699 area2 += a[0] * b[1] - b[0] * a[1];

700 perimeter += Math.hypot(b[0] - a[0], b[1] - a[1]);

701 }

702 const area = Math.abs(area2) / 2;

703 const compactness = perimeter === 0 ? 0 : (4 * Math.PI * area) / (perimeter * perimeter);

704 return { area2, compactness, area };

705}

706

707// Ray cast: a point is inside when a ray to the right crosses the ring an odd number of times.

708function pointInRing(lon, lat, ring) {

709 let inside = false;

710 for (let i = 0, j = ring.length - 1; i < ring.length; j = i++) {

711 const [xi, yi] = ring[i];

712 const [xj, yj] = ring[j];

713 const crosses = yi > lat !== yj > lat;

714 if (!crosses) continue;

715 const x = ((xj - xi) * (lat - yi)) / (yj - yi) + xi;

716 if (lon < x) inside = !inside;

717 }

718 return inside;

719}

720

721// Shortest distance in meters from a lon/lat point to one edge of the polygon.

722function distanceToSegment(lon, lat, a, b) {

723 const cos = Math.cos((lat * Math.PI) / 180);

724 const px = lon * cos * 111320;

725 const py = lat * 110540;

726 const ax = a[0] * cos * 111320;

727 const ay = a[1] * 110540;

728 const bx = b[0] * cos * 111320;

729 const by = b[1] * 110540;

730 const dx = bx - ax;

731 const dy = by - ay;

732 const length2 = dx * dx + dy * dy;

733 let t = length2 === 0 ? 0 : ((px - ax) * dx + (py - ay) * dy) / length2;

734 // Clamp to the segment, not the infinite line.

735 t = Math.max(0, Math.min(1, t));

736 return Math.hypot(px - (ax + t * dx), py - (ay + t * dy));

737}

738

739// Distance from a point to the polygon boundary. Zero when the point is inside.

740function distanceToRing(lon, lat, ring) {

741 let best = Infinity;

742 for (let i = 0; i < ring.length; i++) {

743 const distance = distanceToSegment(lon, lat, ring[i], ring[(i + 1) % ring.length]);

744 if (distance < best) best = distance;

745 }

746 return best;

747}

748

749// The closest OSM place that is inside the ring or within OSM_METERS of it.

750// Places farther away are left for the satellite score.

751function nearestPlace(ring, places) {

752 let best = null;

753 let bestDistance = Infinity;

754 for (const place of places) {

755 const distance = pointInRing(place.lon, place.lat, ring)

756 ? 0

757 : distanceToRing(place.lon, place.lat, ring);

758 if (distance < bestDistance) {

759 bestDistance = distance;

760 best = place;

761 }

762 }

763 if (!best || bestDistance > OSM_METERS) return null;

764 return best;

765}

766

767// Satellite-only certainty. Weights live here so they are easy to change.

768// 45% building count : 8 or more roof clusters fills this part

769// 30% compactness : a circle is 1; 0.45 is a typical tight village

770// 25% roof share : fraction of the blob that passed the roof-color test

771// The result is capped at 99. 100 is reserved for an OSM match.

772function satelliteCertainty(buildings, compactness, roofShare) {

773 const buildingScore = Math.min(buildings / 8, 1);

774 const shapeScore = Math.min(compactness / 0.45, 1);

775 const roofScore = Math.min(roofShare / 0.35, 1);

776 const score = 100 * (0.45 * buildingScore + 0.3 * shapeScore + 0.25 * roofScore);

777 return Math.min(99, Math.round(score));

778}

779

780// Ask Overpass for JSON. Alternates servers because the public ones often answer "too busy".

781// When each server can be tried again, and how many times in a row it has failed.

782const overpassCooldown = new Map();

783const overpassFailures = new Map();

784

785const sleep = (ms) => new Promise((resolve) => setTimeout(resolve, ms));

786

787// Seconds until this connection gets a free slot on a server with a status page:

788// 0 when a slot is free now, a number when it says "in N seconds", null when it can't tell.

789async function overpassSlotWait(server) {

790 try {

791 const response = await fetch(`${server.url}/status`, {

792 headers: { "User-Agent": OVERPASS_USER_AGENT },

793 signal: AbortSignal.timeout(15_000),

794 });

795 const text = await response.text();

796 if (/\d+ slots? available now/.test(text)) return 0;

797 const waits = [...text.matchAll(/in (-?\d+) seconds/g)].map((m) => Math.max(0, Number(m[1])));

798 return waits.length ? Math.min(...waits) : null;

799 } catch {

800 return null;

801 }

802}

803

804// Rest a failing server: as long as its status page says, else longer after each failure.

805async function coolDown(server) {

806 const failures = (overpassFailures.get(server.url) ?? 0) + 1;

807 overpassFailures.set(server.url, failures);

808 let ms = Math.min(15_000 * 2 ** (failures - 1), OVERPASS_MAX_COOLDOWN_MS);

809 if (server.status) {

810 const wait = await overpassSlotWait(server);

811 if (wait !== null) ms = Math.min(Math.max(wait * 1000, 5_000), OVERPASS_MAX_COOLDOWN_MS);

812 }

813 overpassCooldown.set(server.url, Date.now() + ms);

814 console.error(`overpass: ${new URL(server.url).host} busy, resting ${Math.round(ms / 1000)}s`);

815}

816

817// The server to try next: the first one not resting, or else the one that rests the least.

818function nextServer() {

819 const now = Date.now();

820 const ready = OVERPASS_SERVERS.find((server) => (overpassCooldown.get(server.url) ?? 0) <= now);

821 if (ready) return { server: ready, waitMs: 0 };

822 const soonest = OVERPASS_SERVERS.reduce((a, b) =>

823 (overpassCooldown.get(a.url) ?? 0) <= (overpassCooldown.get(b.url) ?? 0) ? a : b,

824 );

825 return { server: soonest, waitMs: (overpassCooldown.get(soonest.url) ?? 0) - now };

826}

827

828// Ask Overpass for JSON, spreading the load over several public servers. A busy server

829// rests (for as long as its status page says, when it has one) while the others are tried.

830async function fetchOverpass(query) {

831 let lastError = new Error("Overpass is busy or returned an error");

832 for (let attempt = 0; attempt < OVERPASS_ATTEMPTS; attempt++) {

833 const { server, waitMs } = nextServer();

834 if (waitMs > 0) await sleep(waitMs);

835 // The main server tells us when a slot frees up: wait for it instead of being refused.

836 if (server.status) {

837 const wait = await overpassSlotWait(server);

838 if (wait && wait > 60) {

839 overpassCooldown.set(server.url, Date.now() + wait * 1000);

840 console.error(`overpass: ${new URL(server.url).host} has no free slot for ${wait}s, trying a mirror`);

841 continue;

842 }

843 if (wait) {

844 console.error(`overpass: waiting ${wait}s for a free slot on ${new URL(server.url).host}`);

845 await sleep(wait * 1000 + 500);

846 }

847 }

848 try {

849 const response = await fetch(`${server.url}/interpreter`, {

850 method: "POST",

851 headers: {

852 "Content-Type": "application/x-www-form-urlencoded",

853 "User-Agent": OVERPASS_USER_AGENT,

854 },

855 body: new URLSearchParams({ data: query }),

856 signal: AbortSignal.timeout(OVERPASS_REQUEST_MS),

857 });

858 const text = await response.text();

859 // A busy Overpass returns an HTML error page instead of JSON.

860 if (!response.ok || text[0] !== "{") {

861 throw new Error(`Overpass is busy or returned an error (${new URL(server.url).host} ${response.status})`);

862 }

863 overpassFailures.delete(server.url);

864 return JSON.parse(text).elements || [];

865 } catch (error) {

866

      lastError = /Overpass/.test(error.message) ? error : new Error(`Overpass is busy or returned an error (${error.message})`);

867 await coolDown(server);

868 }

869 }

870 throw lastError;

871}

872

873// Ask Overpass for settlement nodes in the bbox.

874async function fetchOsmPlaces(bbox) {

875 const local = localRegionFor(bbox);

876 if (local) {

877 return localJson(local, "places.json")

878 .filter((p) => p.lon >= bbox.west && p.lon <= bbox.east && p.lat >= bbox.south && p.lat <= bbox.north)

879 .map((p) => ({ id: p.id, name: p.name, lon: p.lon, lat: p.lat }));

880 }

881 const query = `[out:json][timeout:60];

882node["place"~"^(village|hamlet|isolated_dwelling|town|suburb)$"](${bbox.south},${bbox.west},${bbox.north},${bbox.east});

883out body;`;

884 const elements = await fetchOverpass(query);

885 return elements

886 .filter((element) => element.type === "node")

887 .map((element) => ({

888 id: element.id,

889 name: element.tags && element.tags.name ? element.tags.name : null,

890 lon: element.lon,

891 lat: element.lat,

892 }));

893}

894

895// --- Local OpenStreetMap data ---------------------------------------------

896// prepare-osm.js turns a Geofabrik extract into osm-data/<region>/: places, clinics and

897// road tiles. When a box is inside a prepared region, those answer instead of Overpass:

898// no rate limits, no "busy", and seconds instead of minutes. Elsewhere Overpass is used.

899

900// Run from the repo root (scan-batch.js) or from web/ (the map server).

901

const OSM_DIRS = [process.env.HESTIA_OSM_DIR, join(process.cwd(), "osm-data"), join(process.cwd(), "..", "osm-data")].filter(Boolean);

902let localRegions = null;

903const localCache = new Map();

904

905function loadLocalRegions() {

906 if (localRegions) return localRegions;

907 localRegions = [];

908 for (const dir of OSM_DIRS) {

909 if (!existsSync(dir)) continue;

910 for (const name of readdirSync(dir)) {

911 const metaFile = join(dir, name, "meta.json");

912 if (!existsSync(metaFile)) continue;

913 const meta = JSON.parse(readFileSync(metaFile, "utf8"));

914 localRegions.push({ dir: join(dir, name), meta });

915 console.error(`osm data: local ${meta.region} (built ${meta.builtAt.slice(0, 10)})`);

916 }

917 if (localRegions.length) break;

918 }

919 return localRegions;

920}

921

922// The prepared region whose extract contains the center of this box, if any.

923function localRegionFor(box) {

924 const lon = (box.west + box.east) / 2;

925 const lat = (box.south + box.north) / 2;

926 return (

927 loadLocalRegions().find(

928 ({ meta: { bbox } }) => lon >= bbox.west && lon <= bbox.east && lat >= bbox.south && lat <= bbox.north,

929 ) ?? null

930 );

931}

932

933function localJson(region, file) {

934 const path = join(region.dir, file);

935 if (!localCache.has(path)) localCache.set(path, JSON.parse(readFileSync(path, "utf8")));

936 return localCache.get(path);

937}

938

939// One road tile: [[id, highway, surface, oneway, [[lon, lat], …]], …]. Missing tile: no roads.

940function localTile(path) {

941 if (!localCache.has(path)) {

942 let ways = [];

943 if (existsSync(path)) {

944 ways = gunzipSync(readFileSync(path))

945 .toString()

946 .split("\n")

947 .filter(Boolean)

948 .map((line) => JSON.parse(line));

949 }

950 localCache.set(path, ways);

951 }

952 return localCache.get(path);

953}

954

955// Roads touching the boxes, shaped like Overpass ways. Roads that meet share the exact

956// same coordinate, so the coordinate doubles as the node id that joins them in the graph.

957function localRoads(region, boxes, mainOnly) {

958 const size = mainOnly ? region.meta.mainTile : region.meta.roadTile;

959 const folder = mainOnly ? "main" : "roads";

960 const ways = new Map();

961 for (const box of boxes) {

962 for (let x = Math.floor(box.west / size); x <= Math.floor(box.east / size); x++) {

963 for (let y = Math.floor(box.south / size); y <= Math.floor(box.north / size); y++) {

964

        for (const [id, highway, surface, oneway, coords] of localTile(join(region.dir, folder, `${x}_${y}.ndjson.gz`))) {

965 if (ways.has(id)) continue;

966

          if (!coords.some(([lon, lat]) => lon >= box.west && lon <= box.east && lat >= box.south && lat <= box.north)) {

967 continue;

968 }

969 const tags = { highway };

970 if (surface) tags.surface = surface;

971 if (oneway) tags.oneway = oneway;

972 ways.set(id, {

973 type: "way",

974 id,

975 nodes: coords.map(([lon, lat]) => `${lon},${lat}`),

976 geometry: coords.map(([lon, lat]) => ({ lon, lat })),

977 tags,

978 });

979 }

980 }

981 }

982 }

983 return [...ways.values()];

984}

985

986// --- Cost and route --------------------------------------------------------

987// Predicted cost is the dollars to vaccinate a village plus the dollars to

988// get a team there. The recommendation is the polygon where that sum is lowest.

989

990// Local distance in meters. Accurate enough inside one clinic's reach.

991function metersBetween(lon1, lat1, lon2, lat2) {

992 const cos = Math.cos((((lat1 + lat2) / 2) * Math.PI) / 180);

993 const dx = (lon2 - lon1) * 111320 * cos;

994 const dy = (lat2 - lat1) * 110540;

995 return Math.hypot(dx, dy);

996}

997

998// Grow a bbox by a fixed number of degrees, staying inside Web Mercator.

999function padBbox(bbox, degrees) {

1000 return {

1001 south: Math.max(-85, bbox.south - degrees),

1002 west: Math.max(-180, bbox.west - degrees),

1003 north: Math.min(85, bbox.north + degrees),

1004 east: Math.min(180, bbox.east + degrees),

1005 };

1006}

1007

1008// Five people per roof cluster, or five per 250 m², whichever is larger.

1009// The area figure covers adjoining roofs that fused into one cluster.

1010function estimatePopulation(buildings, areaM2) {

1011 const fromRoofs = buildings * PEOPLE_PER_BUILDING;

1012 const fromArea = Math.round(areaM2 / METERS_PER_DWELLING) * PEOPLE_PER_BUILDING;

1013 return Math.max(fromRoofs, fromArea, PEOPLE_PER_BUILDING);

1014}

1015

1016// Dollars to vaccinate everyone in the settlement, including the stop itself.

1017function healthcareCost(population) {

1018 return population * HEALTHCARE_USD_PER_PERSON + HEALTHCARE_USD_PER_STOP;

1019}

1020

1021// Average of the ring corners. The repeated closing point is left out.

1022function ringCentroid(ring) {

1023 const count = ring.length > 1 && ring[0][0] === ring[ring.length - 1][0] &&

1024 ring[0][1] === ring[ring.length - 1][1]

1025 ? ring.length - 1

1026 : ring.length;

1027 let lon = 0;

1028 let lat = 0;

1029 for (let i = 0; i < count; i++) {

1030 lon += ring[i][0];

1031 lat += ring[i][1];

1032 }

1033 return [lon / count, lat / count];

1034}

1035

1036// Hospitals, clinics, doctors, and health posts. Ways use their center point.

1037// The nearest facilities, searching wider and wider until at least one turns up.

1038async function fetchClinics(bbox) {

1039 for (const pad of CLINIC_SEARCH_PADS_DEG) {

1040 const clinics = await fetchClinicsWithin(padBbox(bbox, pad));

1041 if (clinics.length) {

1042 console.error(`clinic search: ${clinics.length} within ${Math.round(pad * 111)} km`);

1043 return clinics;

1044 }

1045 }

1046 return [];

1047}

1048

1049// Hospitals, clinics, doctors, and health posts inside one box.

1050async function fetchClinicsWithin(box) {

1051 const local = localRegionFor(box);

1052 if (local) {

1053 const found = localJson(local, "clinics.json")

1054 .filter((c) => c.lon >= box.west && c.lon <= box.east && c.lat >= box.south && c.lat <= box.north)

1055 .map((c) => ({ osmId: c.id, name: c.name, kind: c.kind, lon: c.lon, lat: c.lat }));

1056 return mergeSameFacility(found);

1057 }

1058 const query = `[out:json][timeout:90];

1059(

1060 node["amenity"~"^(hospital|clinic|doctors|health_post)$"](${box.south},${box.west},${box.north},${box.east});

1061 way["amenity"~"^(hospital|clinic|doctors|health_post)$"](${box.south},${box.west},${box.north},${box.east});

1062 node["healthcare"~"^(hospital|clinic|doctor|centre|center)$"](${box.south},${box.west},${box.north},${box.east});

1063);

1064out center;`;

1065 const elements = await fetchOverpass(query);

1066 const clinics = [];

1067 for (const element of elements) {

1068 const lat = Number.isFinite(element.lat) ? element.lat : element.center && element.center.lat;

1069 const lon = Number.isFinite(element.lon) ? element.lon : element.center && element.center.lon;

1070 if (!Number.isFinite(lat) || !Number.isFinite(lon)) continue;

1071 const tags = element.tags || {};

1072 clinics.push({

1073 osmId: element.id,

1074 name: tags.name || null,

1075 kind: tags.amenity || tags.healthcare || "clinic",

1076 lon,

1077 lat,

1078 });

1079 }

1080 return mergeSameFacility(clinics);

1081}

1082

1083// A building outline and its entrance node are the same facility.

1084function mergeSameFacility(clinics) {

1085 const kept = [];

1086 for (const clinic of clinics) {

1087 const match = kept.find(

1088 (other) => metersBetween(other.lon, other.lat, clinic.lon, clinic.lat) < 40

1089 );

1090 if (!match) {

1091 kept.push(clinic);

1092 continue;

1093 }

1094 if (!match.name && clinic.name) {

1095 match.name = clinic.name;

1096 match.osmId = clinic.osmId;

1097 match.kind = clinic.kind;

1098 }

1099 }

1100 return kept;

1101}

1102

1103// The few clinics actually worth routing to, so the road download stays local.

1104function clinicsNearVillages(features, clinics) {

1105 const chosen = new Map();

1106 for (const feature of features) {

1107 const [lon, lat] = ringCentroid(feature.geometry.coordinates[0]);

1108 const nearest = clinics

1109 .map((clinic) => ({

1110 clinic,

1111 meters: metersBetween(lon, lat, clinic.lon, clinic.lat),

1112 }))

1113 .sort((a, b) => a.meters - b.meters)

1114 .slice(0, CLINICS_PER_VILLAGE);

1115 for (const item of nearest) chosen.set(item.clinic.osmId, item.clinic);

1116 }

1117 return [...chosen.values()];

1118}

1119

1120// Box that contains the villages and the clinics the routes will start from.

1121function bboxAround(bbox, clinics) {

1122 let south = bbox.south;

1123 let west = bbox.west;

1124 let north = bbox.north;

1125 let east = bbox.east;

1126 for (const clinic of clinics) {

1127 south = Math.min(south, clinic.lat);

1128 north = Math.max(north, clinic.lat);

1129 west = Math.min(west, clinic.lon);

1130 east = Math.max(east, clinic.lon);

1131 }

1132 // A little extra so a road just past a clinic is still in the graph.

1133 return padBbox({ south, west, north, east }, 0.01);

1134}

1135

1136// Drivable roads, including tracks. Footpaths are not a vaccination route.

1137const ALL_ROADS =

1138

  "motorway|motorway_link|trunk|trunk_link|primary|primary_link|secondary|secondary_link|tertiary|tertiary_link|unclassified|residential|living_street|service|track";

1139// Roads a team uses to cover long distances between a far clinic and the settlements.

1140

const MAIN_ROADS = "motorway|motorway_link|trunk|trunk_link|primary|primary_link|secondary|secondary_link|tertiary|tertiary_link";

1141

1142// Roads of the given types inside any of the boxes, in a single Overpass request.

1143async function fetchRoads(boxes, types = ALL_ROADS) {

1144 const local = boxes.length ? localRegionFor(boxes[0]) : null;

1145 if (local && boxes.every((box) => localRegionFor(box) === local)) {

1146 return localRoads(local, boxes, types === MAIN_ROADS);

1147 }

1148 const parts = boxes

1149 .map((b) => ` way["highway"~"^(${types})$"](${b.south},${b.west},${b.north},${b.east});`)

1150 .join("\n");

1151 const query = `[out:json][timeout:180];

1152(

1153${parts}

1154);

1155out geom;`;

1156 const elements = await fetchOverpass(query);

1157 return elements.filter((element) => element.type === "way");

1158}

1159

1160// Every road near the boxes, plus main roads across the whole area out to the clinics.

1161// `boxes` is one or more scan boxes (a joined clinic tour spans several).

1162export async function fetchRoadNetwork(boxes, clinics) {

1163 const ways = new Map();

1164 const add = (list) => list.forEach((way) => ways.set(way.id, way));

1165 add(await fetchRoads(boxes.map((box) => padBbox(box, LOCAL_ROADS_PAD_DEG))));

1166 const all = boxes.reduce((acc, box) => ({

1167 south: Math.min(acc.south, box.south),

1168 west: Math.min(acc.west, box.west),

1169 north: Math.max(acc.north, box.north),

1170 east: Math.max(acc.east, box.east),

1171 }));

1172 const extent = bboxAround(all, clinics);

1173 const local = boxes.length === 1 ? padBbox(boxes[0], LOCAL_ROADS_PAD_DEG) : null;

1174 const reachesOut =

1175

    !local || extent.south < local.south || extent.west < local.west || extent.north > local.north || extent.east > local.east;

1176 if (reachesOut) {

1177 add(await fetchRoads([extent], MAIN_ROADS));

1178 // The clinic itself may sit on a small road, so add every road right around each one.

1179 if (clinics.length) {

1180 add(

1181 await fetchRoads(

1182 clinics.map((c) => padBbox({ south: c.lat, west: c.lon, north: c.lat, east: c.lon }, LOCAL_ROADS_PAD_DEG)),

1183 ),

1184 );

1185 }

1186 }

1187 return [...ways.values()];

1188}

1189

1190// Dollars per kilometer. A track costs more because it is slower and harder on the vehicle.

1191function roadUsdPerKm(tags) {

1192 const highway = String(tags.highway || "").replace(/_link$/, "");

1193 const rates = {

1194 motorway: 1.2,

1195 trunk: 1.3,

1196 primary: 1.4,

1197 secondary: 1.6,

1198 tertiary: 1.8,

1199 unclassified: 2.2,

1200 residential: 2,

1201 living_street: 2.2,

1202 service: 2.4,

1203 track: 4.5,

1204 };

1205 let perKm = rates[highway] ?? 2.5;

1206 const surface = tags.surface || "";

1208 if (unpaved && highway !== "track") perKm *= 1.5;

1209 const slow = highway === "track" || unpaved;

1210 return { perKm, slow };

1211}

1212

1213function addDirected(adj, from, to, meters, usd, slow) {

1214 let list = adj.get(from);

1215 if (!list) {

1216 list = [];

1217 adj.set(from, list);

1218 }

1219 list.push({ to, meters, usd, slow });

1220}

1221

1222// OSM ways become a graph. Intersections share node ids, so ways join there.

1223function buildRoadGraph(ways) {

1224 const nodes = new Map();

1225 const adj = new Map();

1226 const segments = [];

1227 for (const way of ways) {

1228 if (!way.nodes || !way.geometry || way.nodes.length < 2) continue;

1229 if (way.nodes.length !== way.geometry.length) continue;

1230 const tags = way.tags || {};

1231 const { perKm, slow } = roadUsdPerKm(tags);

1232 const oneway = tags.oneway;

1233 const forward = oneway !== "-1" && oneway !== "reverse";

1234 const backward = oneway !== "yes" && oneway !== "true" && oneway !== "1";

1235 for (let i = 0; i < way.nodes.length; i++) {

1236 const point = way.geometry[i];

1237 if (!nodes.has(way.nodes[i])) {

1238 nodes.set(way.nodes[i], { lon: point.lon, lat: point.lat });

1239 }

1240 }

1241 for (let i = 1; i < way.nodes.length; i++) {

1242 const a = way.nodes[i - 1];

1243 const b = way.nodes[i];

1244 const pa = nodes.get(a);

1245 const pb = nodes.get(b);

1246 const meters = metersBetween(pa.lon, pa.lat, pb.lon, pb.lat);

1247 if (meters === 0) continue;

1248 const usd = (meters / 1000) * perKm;

1249 if (forward) addDirected(adj, a, b, meters, usd, slow);

1250 if (backward) addDirected(adj, b, a, meters, usd, slow);

1251 segments.push({ a, b, meters, usd, slow, forward, backward });

1252 }

1253 }

1254 return { nodes, adj, segments, nextTemp: -1 };

1255}

1256

1257// Closest point on a road segment, in meters, plus where it falls (0 at a, 1 at b).

1258function projectToSegment(lon, lat, a, b) {

1259 const cos = Math.cos((lat * Math.PI) / 180);

1260 const px = lon * cos * 111320;

1261 const py = lat * 110540;

1262 const ax = a.lon * cos * 111320;

1263 const ay = a.lat * 110540;

1264 const bx = b.lon * cos * 111320;

1265 const by = b.lat * 110540;

1266 const dx = bx - ax;

1267 const dy = by - ay;

1268 const length2 = dx * dx + dy * dy;

1269 let t = length2 === 0 ? 0 : ((px - ax) * dx + (py - ay) * dy) / length2;

1270 t = Math.max(0, Math.min(1, t));

1271 const sx = ax + t * dx;

1272 const sy = ay + t * dy;

1273 return {

1274 t,

1275 meters: Math.hypot(px - sx, py - sy),

1276 lon: sx / (cos * 111320),

1277 lat: sy / 110540,

1278 };

1279}

1280

1281// Join a clinic or village to the road graph at the nearest point on a segment.

1282// Returns null when every road is farther than MAX_SNAP_METERS.

1283function snapToRoad(graph, lon, lat) {

1284 if (graph.segments.length === 0) return null;

1285 let best = null;

1286 let bestSegment = null;

1287 for (const segment of graph.segments) {

1288 const snap = projectToSegment(lon, lat, graph.nodes.get(segment.a), graph.nodes.get(segment.b));

1289 if (!best || snap.meters < best.meters) {

1290 best = snap;

1291 bestSegment = segment;

1292 }

1293 }

1294 if (!best || best.meters > MAX_SNAP_METERS) return null;

1295 // Already on a vertex. Reuse that node instead of splitting the edge.

1296 if (best.t <= 0.02) {

1297 const node = graph.nodes.get(bestSegment.a);

1298 return { id: bestSegment.a, lon: node.lon, lat: node.lat, meters: best.meters };

1299 }

1300 if (best.t >= 0.98) {

1301 const node = graph.nodes.get(bestSegment.b);

1302 return { id: bestSegment.b, lon: node.lon, lat: node.lat, meters: best.meters };

1303 }

1304 const id = graph.nextTemp;

1305 graph.nextTemp -= 1;

1306 graph.nodes.set(id, { lon: best.lon, lat: best.lat });

1307 const metersA = bestSegment.meters * best.t;

1308 const metersB = bestSegment.meters * (1 - best.t);

1309 const usdA = bestSegment.usd * best.t;

1310 const usdB = bestSegment.usd * (1 - best.t);

1311 // Moving toward A is against the way direction, so it follows `backward`.

1312 if (bestSegment.backward && metersA > 0) {

1313 addDirected(graph.adj, id, bestSegment.a, metersA, usdA, bestSegment.slow);

1314 }

1315 if (bestSegment.forward && metersA > 0) {

1316 addDirected(graph.adj, bestSegment.a, id, metersA, usdA, bestSegment.slow);

1317 }

1318 if (bestSegment.forward && metersB > 0) {

1319 addDirected(graph.adj, id, bestSegment.b, metersB, usdB, bestSegment.slow);

1320 }

1321 if (bestSegment.backward && metersB > 0) {

1322 addDirected(graph.adj, bestSegment.b, id, metersB, usdB, bestSegment.slow);

1323 }

1324 return { id, lon: best.lon, lat: best.lat, meters: best.meters };

1325}

1326

1327function heapPush(heap, item) {

1328 heap.push(item);

1329 let i = heap.length - 1;

1330 while (i > 0) {

1331 const parent = (i - 1) >> 1;

1332 if (heap[parent][0] <= heap[i][0]) break;

1333 const swap = heap[parent];

1334 heap[parent] = heap[i];

1335 heap[i] = swap;

1336 i = parent;

1337 }

1338}

1339

1340function heapPop(heap) {

1341 const top = heap[0];

1342 const last = heap.pop();

1343 if (heap.length === 0) return top;

1344 heap[0] = last;

1345 let i = 0;

1346 for (;;) {

1347 const left = i * 2 + 1;

1348 const right = left + 1;

1349 let smallest = i;

1350 if (left < heap.length && heap[left][0] < heap[smallest][0]) smallest = left;

1351 if (right < heap.length && heap[right][0] < heap[smallest][0]) smallest = right;

1352 if (smallest === i) break;

1353 const swap = heap[i];

1354 heap[i] = heap[smallest];

1355 heap[smallest] = swap;

1356 i = smallest;

1357 }

1358 return top;

1359}

1360

1361// One search from every clinic. The first time a village node is settled,

1362// that path is the cheapest way to reach it.

1363function cheapestRoadTrips(graph, sources, targets) {

1364 const dist = new Map();

1365 const prev = new Map();

1366 const origin = new Map();

1367 const heap = [];

1368 sources.forEach((source, index) => {

1369 const cost = source.offroadUsd;

1370 if (dist.has(source.nodeId) && dist.get(source.nodeId) <= cost) return;

1371 dist.set(source.nodeId, cost);

1372 origin.set(source.nodeId, index);

1373 prev.set(source.nodeId, null);

1374 heapPush(heap, [cost, source.nodeId, index]);

1375 });

1376

1377 const settled = new Set();

1378 const found = new Map();

1379 let remaining = 0;

1380 for (const list of targets.values()) remaining += list.length;

1381

1382 while (heap.length && remaining > 0) {

1383 const [cost, node, sourceIndex] = heapPop(heap);

1384 if (settled.has(node)) continue;

1385 if (dist.get(node) !== cost || origin.get(node) !== sourceIndex) continue;

1386 settled.add(node);

1387 const villages = targets.get(node);

1388 if (villages) {

1389 for (const villageIndex of villages) {

1390 if (found.has(villageIndex)) continue;

1391 found.set(villageIndex, { node, sourceIndex });

1392 remaining -= 1;

1393 }

1394 }

1395 const edges = graph.adj.get(node) || [];

1396 for (const edge of edges) {

1397 if (settled.has(edge.to)) continue;

1398 const next = cost + edge.usd;

1399 if (dist.has(edge.to) && dist.get(edge.to) <= next) continue;

1400 dist.set(edge.to, next);

1401 origin.set(edge.to, sourceIndex);

1402 prev.set(edge.to, { node, meters: edge.meters, usd: edge.usd, slow: edge.slow });

1403 heapPush(heap, [next, edge.to, sourceIndex]);

1404 }

1405 }

1406 return { found, prev };

1407}

1408

1409// Walk predecessors back to the clinic's road node.

1410function traceRoad(graph, prev, node) {

1411 const points = [];

1412 let roadMeters = 0;

1413 let slowMeters = 0;

1414 let roadUsd = 0;

1415 let cursor = node;

1416 const seen = new Set();

1417 const guard = prev.size + 2;

1418 for (let step = 0; step < guard; step++) {

1419 if (seen.has(cursor)) break;

1420 seen.add(cursor);

1421 const point = graph.nodes.get(cursor);

1422 points.push([point.lon, point.lat]);

1423 const from = prev.get(cursor);

1424 if (!from) break;

1425 roadMeters += from.meters;

1426 roadUsd += from.usd;

1427 if (from.slow) slowMeters += from.meters;

1428 cursor = from.node;

1429 }

1430 points.reverse();

1431 return { points, roadMeters, slowMeters, roadUsd };

1432}

1433

1434// 0 is a short paved hop. 100 is a long way off the road network.

1435// The off-road gap at the village itself is half the score.

1436function reachDifficulty(villageOffroadMeters, travelMeters, slowMeters) {

1437 const offroad = Math.min(villageOffroadMeters / 1000, 1) * 50;

1438 const distance = Math.min(travelMeters / 30000, 1) * 30;

1439 const share = travelMeters === 0 ? 1 : slowMeters / travelMeters;

1440 const surface = Math.min(share, 1) * 20;

1441 const score = Math.max(0, Math.min(100, Math.round(offroad + distance + surface)));

1442 const label = score <= 33 ? "easy" : score <= 66 ? "moderate" : "hard";

1443 return { score, label };

1444}

1445

1446function pushCoord(coordinates, lon, lat) {

1447 const point = [round6(lon), round6(lat)];

1448 const last = coordinates[coordinates.length - 1];

1449 if (last && last[0] === point[0] && last[1] === point[1]) return;

1450 coordinates.push(point);

1451}

1452

1453// Direct line, priced entirely as off-road. Used when no road connects the two.

1454function straightTrip(clinic, lon, lat) {

1455 const meters = metersBetween(clinic.lon, clinic.lat, lon, lat);

1456 const reachUsd = Math.round((meters / 1000) * OFFROAD_USD_PER_KM);

1457 const difficulty = reachDifficulty(meters, meters, meters);

1458 return {

1459 reachUsd,

1460 reachDifficulty: difficulty.score,

1461 reachLabel: difficulty.label,

1462 travelMeters: meters,

1463 offroadMeters: meters,

1464 slowMeters: meters,

1465 route: [

1466 [round6(clinic.lon), round6(clinic.lat)],

1467 [round6(lon), round6(lat)],

1468 ],

1469 clinic,

1470 };

1471}

1472

1473function tripFromRoad(clinic, source, village, traced, destLon, destLat) {

1474 const travelMeters = traced.roadMeters + source.meters + village.meters;

1475 const slowMeters = traced.slowMeters + source.meters + village.meters;

1476 const reachUsd = Math.round(traced.roadUsd + source.offroadUsd + (village.meters / 1000) * OFFROAD_USD_PER_KM);

1477 const difficulty = reachDifficulty(village.meters, travelMeters, slowMeters);

1478 const route = [];

1479 pushCoord(route, clinic.lon, clinic.lat);

1480 for (const [lon, lat] of traced.points) pushCoord(route, lon, lat);

1481 pushCoord(route, destLon, destLat);

1482 return {

1483 reachUsd,

1484 reachDifficulty: difficulty.score,

1485 reachLabel: difficulty.label,

1486 travelMeters,

1487 offroadMeters: village.meters,

1488 slowMeters,

1489 route,

1490 clinic,

1491 };

1492}

1493

1494function cheaperTrip(current, candidate) {

1495 if (!current || candidate.reachUsd < current.reachUsd) return candidate;

1496 return current;

1497}

1498

1499// Fill reach fields on every polygon and return the one with the lowest predicted cost.

1500function scoreAccess(features, clinics, graph) {

1501 const arrivals = features.map((feature) => ringCentroid(feature.geometry.coordinates[0]));

1502 const clinicChoices = arrivals.map(([lon, lat]) =>

1503 clinics

1504 .map((clinic) => ({ clinic, meters: metersBetween(lon, lat, clinic.lon, clinic.lat) }))

1505 .sort((a, b) => a.meters - b.meters)

1506 .slice(0, CLINICS_PER_VILLAGE)

1507 );

1508

1509 const sources = [];

1510 const sourceKey = new Map();

1511 for (const choices of clinicChoices) {

1512 for (const choice of choices) {

1513 if (sourceKey.has(choice.clinic.osmId)) continue;

1514 const snap = snapToRoad(graph, choice.clinic.lon, choice.clinic.lat);

1515 if (!snap) continue;

1516 sourceKey.set(choice.clinic.osmId, sources.length);

1517 sources.push({

1518 clinic: choice.clinic,

1519 nodeId: snap.id,

1520 meters: snap.meters,

1521 offroadUsd: (snap.meters / 1000) * OFFROAD_USD_PER_KM,

1522 });

1523 }

1524 }

1525

1526 const villageSnaps = arrivals.map(([lon, lat]) => snapToRoad(graph, lon, lat));

1527 const targets = new Map();

1528 villageSnaps.forEach((snap, index) => {

1529 if (!snap) return;

1530 let list = targets.get(snap.id);

1531 if (!list) {

1532 list = [];

1533 targets.set(snap.id, list);

1534 }

1535 list.push(index);

1536 });

1537

1538 const road = sources.length && targets.size

1539 ? cheapestRoadTrips(graph, sources, targets)

1540 : { found: new Map(), prev: new Map() };

1541

1542 return features.map((feature, index) => {

1543 const [lon, lat] = arrivals[index];

1544 let best = null;

1545 const hit = road.found.get(index);

1546 if (hit && villageSnaps[index]) {

1547 const source = sources[hit.sourceIndex];

1548 const traced = traceRoad(graph, road.prev, hit.node);

1549 best = tripFromRoad(source.clinic, source, villageSnaps[index], traced, lon, lat);

1550 }

1551 for (const choice of clinicChoices[index]) {

1552 best = cheaperTrip(best, straightTrip(choice.clinic, lon, lat));

1553 }

1554 if (!best) {

1555 // No clinic mapped anywhere in the widest search. Price that whole distance as off-road.

1556 const meters = CLINIC_PAD_DEG * 110540;

1557 const difficulty = reachDifficulty(meters, meters, meters);

1558 best = {

1559 reachUsd: Math.round((meters / 1000) * OFFROAD_USD_PER_KM),

1560 reachDifficulty: difficulty.score,

1561 reachLabel: difficulty.label,

1562 travelMeters: meters,

1563 offroadMeters: meters,

1564 slowMeters: meters,

1565 route: [[round6(lon), round6(lat)]],

1566 clinic: null,

1567 };

1568 }

1569 const population = feature.properties.population;

1570 const healthcareUsd = feature.properties.healthcareUsd;

1571 return {

1572 population,

1573 healthcareUsd,

1574 reachDifficulty: best.reachDifficulty,

1575 reachLabel: best.reachLabel,

1576 reachUsd: best.reachUsd,

1577 travelMeters: Math.round(best.travelMeters),

1578 offroadMeters: Math.round(best.offroadMeters),

1579 slowMeters: Math.round(best.slowMeters),

1580 predictedCostUsd: healthcareUsd + best.reachUsd,

1581 route: best.route,

1582 clinic: best.clinic

1583 ? {

1584 name: best.clinic.name,

1585 lon: round6(best.clinic.lon),

1586 lat: round6(best.clinic.lat),

1587 osmId: best.clinic.osmId,

1588 }

1589 : null,

1590 };

1591 });

1592}

1593

1594function isLowerCost(score, other) {

1595 if (score.predictedCostUsd !== other.predictedCostUsd) {

1596 return score.predictedCostUsd < other.predictedCostUsd;

1597 }

1598 if (score.reachDifficulty !== other.reachDifficulty) {

1599 return score.reachDifficulty < other.reachDifficulty;

1600 }

1601 return score.population > other.population;

1602}

1603

1604// Route every village, then keep the polygon whose predicted cost is lowest.

1605async function recommendAccess(features, bbox) {

1606 if (features.length === 0) return { recommendation: null, tour: null };

1607

1608 let clinics = [];

1609 try {

1610 clinics = await fetchClinics(bbox);

1611 } catch (error) {

1612 console.error(`clinics unavailable: ${error.message || error}`);

1613 }

1614 console.error(`clinics ${clinics.length}`);

1615

1616 const routed = clinicsNearVillages(features, clinics);

1617 let graph = buildRoadGraph([]);

1618 // Roads are needed even with no clinic: the team tour drives between settlements.

1619 try {

1620 const ways = await fetchRoadNetwork([bbox], routed);

1621 graph = buildRoadGraph(ways);

1622 console.error(`road segments ${graph.segments.length}`);

1623 } catch (error) {

1624 console.error(`roads unavailable: ${error.message || error}`);

1625 }

1626

1627 const scores = scoreAccess(features, routed, graph);

1628 let bestIndex = 0;

1629 for (let i = 0; i < scores.length; i++) {

1630 const score = scores[i];

1631 features[i].properties.reachDifficulty = score.reachDifficulty;

1632 features[i].properties.reachLabel = score.reachLabel;

1633 features[i].properties.reachUsd = score.reachUsd;

1634 // The three inputs to reachDifficulty, so a map can explain the score.

1635 features[i].properties.travelMeters = score.travelMeters;

1636 features[i].properties.offroadMeters = score.offroadMeters;

1637 features[i].properties.slowMeters = score.slowMeters;

1638 features[i].properties.predictedCostUsd = score.predictedCostUsd;

1639 // Each polygon keeps its own trip so a map can draw every route, not just the winner.

1640 features[i].properties.route = { type: "LineString", coordinates: score.route };

1641 features[i].properties.clinic = score.clinic;

1642 if (isLowerCost(score, scores[bestIndex])) bestIndex = i;

1643 }

1644

1645 const best = scores[bestIndex];

1646 features[bestIndex].properties.recommended = true;

1647 const clinicName = best.clinic ? best.clinic.name || "unnamed clinic" : "no clinic in range";

1648 console.error(

1649 `lowest cost polygon: $${best.predictedCostUsd} ` +

1650 `(healthcare $${best.healthcareUsd}, reach $${best.reachUsd}, ` +

1651 `${best.reachLabel} ${best.reachDifficulty}/100, population ${best.population}, ${clinicName})`

1652 );

1653 const recommendation = {

1654 polygon: features[bestIndex].geometry.coordinates[0],

1655 population: best.population,

1656 healthcareUsd: best.healthcareUsd,

1657 reachDifficulty: best.reachDifficulty,

1658 reachLabel: best.reachLabel,

1659 reachUsd: best.reachUsd,

1660 predictedCostUsd: best.predictedCostUsd,

1661 route: { type: "LineString", coordinates: best.route },

1662 clinic: best.clinic,

1663 };

1664 const tour = planTour(features, scores, routed, graph);

1665 return { recommendation, tour };

1666}

1667

1668// --- Team tour -------------------------------------------------------------

1669// One vaccination team leaves a clinic, visits every settlement, and returns.

1670// Each leg costs the cheaper of driving (roads + off-road gaps) or a straight

1671// off-road line. The visiting order is nearest-neighbor, then improved by 2-opt.

1672

1673// Cheapest cost from one road node to every node, stopping once `targets` are settled.

1674function roadCosts(graph, source, targets) {

1675 const dist = new Map([[source, 0]]);

1676 const prev = new Map([[source, null]]);

1677 const heap = [[0, source]];

1678 const settled = new Set();

1679 let remaining = targets.size;

1680 while (heap.length && remaining > 0) {

1681 const [cost, node] = heapPop(heap);

1682 if (settled.has(node) || dist.get(node) !== cost) continue;

1683 settled.add(node);

1684 if (targets.has(node)) remaining -= 1;

1685 for (const edge of graph.adj.get(node) || []) {

1686 const next = cost + edge.usd;

1687 if (settled.has(edge.to) || (dist.has(edge.to) && dist.get(edge.to) <= next)) continue;

1688 dist.set(edge.to, next);

1689 prev.set(edge.to, { node, meters: edge.meters, usd: edge.usd, slow: edge.slow });

1690 heapPush(heap, [next, edge.to]);

1691 }

1692 }

1693 return { dist, prev };

1694}

1695

1696// The clinic most settlements already route from, or else the nearest one.

1697function pickDepot(stops, scores, clinics) {

1698 const votes = new Map();

1699 for (const stop of stops) {

1700 const clinic = scores[stop.index].clinic;

1701 if (clinic) votes.set(clinic.osmId, (votes.get(clinic.osmId) || 0) + 1);

1702 }

1703 let best = null;

1704 for (const clinic of clinics) {

1705 const count = votes.get(clinic.osmId) || 0;

1706 if (count && (!best || count > best.count)) best = { clinic, count };

1707 }

1708 return best ? best.clinic : null;

1709}

1710

1711// Total of a closed tour under the symmetric cost matrix.

1712function tourLength(order, cost) {

1713 let total = 0;

1714 for (let i = 0; i < order.length - 1; i++) total += cost[order[i]][order[i + 1]];

1715 return total;

1716}

1717

1718// Reverse order[i..j] whenever that makes the loop cheaper. order[0] and the end stay put.

1719function twoOpt(order, cost) {

1720 for (let pass = 0; pass < TOUR_MAX_PASSES; pass++) {

1721 let improved = false;

1722 for (let i = 1; i < order.length - 2; i++) {

1723 for (let j = i + 1; j < order.length - 1; j++) {

1724 const a = order[i - 1];

1725 const b = order[i];

1726 const c = order[j];

1727 const d = order[j + 1];

1728 const delta = cost[a][c] + cost[b][d] - cost[a][b] - cost[c][d];

1729 if (delta < -1e-9) {

1730 for (let lo = i, hi = j; lo < hi; lo++, hi--) [order[lo], order[hi]] = [order[hi], order[lo]];

1731 improved = true;

1732 }

1733 }

1734 }

1735 if (!improved) break;

1736 }

1737 return order;

1738}

1739

1740// `depot`, when given, is the clinic the tour must start from (a joined clinic tour).

1741function planTour(features, scores, clinics, graph, { depot: fixedDepot = null } = {}) {

1742 const stops = [];

1743 let skipped = 0;

1744 features.forEach((feature, index) => {

1745 if (feature.properties.certainty < TOUR_MIN_CERTAINTY) {

1746 skipped += 1;

1747 return;

1748 }

1749 const [lon, lat] = ringCentroid(feature.geometry.coordinates[0]);

1750 stops.push({ index, lon, lat, snap: snapToRoad(graph, lon, lat) });

1751 });

1752 if (stops.length === 0) return null;

1753

1754 // Node 0 is the depot: a clinic, or the settlement closest to a road when no clinic is in range.

1755 const clinic = fixedDepot ?? pickDepot(stops, scores, clinics);

1756 let depot;

1757 if (clinic) {

1758 depot = {

1759 lon: clinic.lon,

1760 lat: clinic.lat,

1761 snap: snapToRoad(graph, clinic.lon, clinic.lat),

1762

      info: { kind: "clinic", name: clinic.name, osmId: clinic.osmId, lon: round6(clinic.lon), lat: round6(clinic.lat) },

1763 };

1764 } else {

1765 const start = stops.reduce((a, b) => ((a.snap?.meters ?? Infinity) <= (b.snap?.meters ?? Infinity) ? a : b));

1766 depot = {

1767 lon: start.lon,

1768 lat: start.lat,

1769 snap: start.snap,

1770

      info: { kind: "settlement", name: features[start.index].properties.osmName, lon: round6(start.lon), lat: round6(start.lat) },

1771 };

1772 }

1773 const nodes = [depot, ...stops];

1774 const n = nodes.length;

1775 const offroadUsd = (snap) => (snap.meters / 1000) * OFFROAD_USD_PER_KM;

1776

1777 // Directed leg costs. `road[i][j]` is false when the straight off-road line is cheaper.

1778 const targets = new Set(nodes.filter((node) => node.snap).map((node) => node.snap.id));

1779 const direct = Array.from({ length: n }, () => new Float64Array(n));

1780 const road = Array.from({ length: n }, () => new Uint8Array(n));

1781 for (let i = 0; i < n; i++) {

1782 const from = nodes[i];

1783 const reach = from.snap ? roadCosts(graph, from.snap.id, targets).dist : null;

1784 for (let j = 0; j < n; j++) {

1785 if (i === j) continue;

1786 const to = nodes[j];

1787 const straight = (metersBetween(from.lon, from.lat, to.lon, to.lat) / 1000) * OFFROAD_USD_PER_KM;

1788 const viaRoad =

1789 reach && to.snap && reach.has(to.snap.id)

1790 ? offroadUsd(from.snap) + reach.get(to.snap.id) + offroadUsd(to.snap)

1791 : Infinity;

1792 direct[i][j] = Math.min(straight, viaRoad);

1793 road[i][j] = viaRoad < straight ? 1 : 0;

1794 }

1795 }

1796 // Roads are almost all two-way out here, so plan on the cheaper direction of each pair.

1797 const cost = Array.from({ length: n }, (_, i) =>

1798 Array.from({ length: n }, (_, j) => Math.min(direct[i][j], direct[j][i]))

1799 );

1800

1801 // Nearest neighbor from the depot, back to the depot, then 2-opt.

1802 const order = [0];

1803 const left = new Set(Array.from({ length: n - 1 }, (_, k) => k + 1));

1804 while (left.size) {

1805 const last = order[order.length - 1];

1806 let next = null;

1807 for (const k of left) if (next === null || cost[last][k] < cost[last][next]) next = k;

1808 order.push(next);

1809 left.delete(next);

1810 }

1811 order.push(0);

1812 const greedyUsd = tourLength(order, cost);

1813 twoOpt(order, cost);

1814

1815 // Draw each leg and add up road and off-road distance.

1816 const coordinates = [];

1817 let roadMeters = 0;

1818 let roadUsd = 0;

1819 let offroadMeters = 0;

1820 let offroadUsdTotal = 0;

1821 for (let k = 0; k < order.length - 1; k++) {

1822 // Walk the leg in whichever direction the plan priced, then flip it back if needed.

1823 let a = order[k];

1824 let b = order[k + 1];

1825 const flipped = direct[b][a] < direct[a][b];

1826 if (flipped) [a, b] = [b, a];

1827 const from = nodes[a];

1828 const to = nodes[b];

1829 let points;

1830 if (road[a][b]) {

1831 const { prev } = roadCosts(graph, from.snap.id, new Set([to.snap.id]));

1832 const traced = traceRoad(graph, prev, to.snap.id);

1833 points = [[from.lon, from.lat], ...traced.points, [to.lon, to.lat]];

1834 roadMeters += traced.roadMeters;

1835 roadUsd += traced.roadUsd;

1836 offroadMeters += from.snap.meters + to.snap.meters;

1837 offroadUsdTotal += offroadUsd(from.snap) + offroadUsd(to.snap);

1838 } else {

1839 const meters = metersBetween(from.lon, from.lat, to.lon, to.lat);

1840 points = [[from.lon, from.lat], [to.lon, to.lat]];

1841 offroadMeters += meters;

1842 offroadUsdTotal += (meters / 1000) * OFFROAD_USD_PER_KM;

1843 }

1844 if (flipped) points.reverse();

1845 for (const [lon, lat] of points) pushCoord(coordinates, lon, lat);

1846 }

1847

1848 // Visit position on each settlement, 1 = first stop after the depot.

1849 const visit = order.slice(1, -1).map((node) => nodes[node].index);

1850 visit.forEach((featureIndex, position) => {

1851 features[featureIndex].properties.tourStop = position + 1;

1852 });

1853

1854 const costUsd = Math.round(roadUsd + offroadUsdTotal);

1855 console.error(

1856 `team tour: ${stops.length} stops, $${costUsd} ` +

1857 `(${(roadMeters / 1000).toFixed(1)} km road, ${(offroadMeters / 1000).toFixed(1)} km off-road), ` +

1858 `nearest-neighbor $${Math.round(greedyUsd)} improved by 2-opt`

1859 );

1860 return {

1861 depot: depot.info,

1862 minCertainty: TOUR_MIN_CERTAINTY,

1863 stops: stops.length,

1864 skipped,

1865 order: visit,

1866 costUsd,

1867 roadUsd: Math.round(roadUsd),

1868 offroadUsd: Math.round(offroadUsdTotal),

1869 roadMeters: Math.round(roadMeters),

1870 offroadMeters: Math.round(offroadMeters),

1871 greedyUsd: Math.round(greedyUsd),

1872 route: { type: "LineString", coordinates },

1873 };

1874}

1875

1876// Six decimal places is about 0.1 m. Enough for a roof, small enough to read.

1877function round6(value) {

1878 return Math.round(value * 1e6) / 1e6;

1879}

1880

1881// How many pixels in a 0/1 mask are set. Logged so a run can be sanity-checked.

1882function countOnes(mask) {

1883 let count = 0;

1884 for (let i = 0; i < mask.length; i++) count += mask[i];

1885 return count;

1886}

1887

1888// --- Building footprints -----------------------------------------------------

1889

1890// Zoom-9 quadkey of a point: the name Microsoft uses for each footprint file's map square.

1891function quadkey(lon, lat, zoom) {

1892 const n = 2 ** zoom;

1893 const x = Math.min(n - 1, Math.floor(((lon + 180) / 360) * n));

1894 const rad = (lat * Math.PI) / 180;

1895 const y = Math.min(n - 1, Math.floor(((1 - Math.log(Math.tan(rad) + 1 / Math.cos(rad)) / Math.PI) / 2) * n));

1896 let key = "";

1897 for (let i = zoom; i > 0; i--) {

1898 const bit = 1 << (i - 1);

1899 key += String((x & bit ? 1 : 0) + (y & bit ? 2 : 0));

1900 }

1901 return key;

1902}

1903

1904// Every quadkey the box touches. A 4 km box usually sits in one, at most four.

1905function quadkeysForBbox(bbox) {

1906 const keys = new Set();

1907 for (const lon of [bbox.west, bbox.east]) {

1908 for (const lat of [bbox.south, bbox.north]) keys.add(quadkey(lon, lat, FOOTPRINT_QUADKEY_ZOOM));

1909 }

1910 return [...keys];

1911}

1912

1913// Download once, then reuse the copy in the cache folder.

1914async function cachedDownload(url, file, maxAgeMs = Infinity) {

1915 const path = join(FOOTPRINT_CACHE_DIR, file);

1916 try {

1917 const info = await stat(path);

1918 if (Date.now() - info.mtimeMs < maxAgeMs) return readFile(path);

1919 } catch {

1920 // Not cached yet.

1921 }

1922 const response = await fetch(url, { headers: { "User-Agent": "village-polygon-detector/1.0" } });

1923 if (!response.ok) throw new Error(`footprints: ${response.status} for ${url}`);

1924 const buffer = Buffer.from(await response.arrayBuffer());

1925 try {

1926 await mkdir(FOOTPRINT_CACHE_DIR, { recursive: true });

1927 await writeFile(`${path}.tmp`, buffer);

1928 await rename(`${path}.tmp`, path);

1929 } catch {

1930 // No writable cache (unusual). Still fine for this run.

1931 }

1932 return buffer;

1933}

1934

1935let footprintLinks = null;

1936

1937// quadkey → download links. One square can have several files (parts, or country borders).

1938async function loadFootprintLinks() {

1939 if (footprintLinks) return footprintLinks;

1940 const csv = (await cachedDownload(FOOTPRINT_LINKS_URL, "dataset-links.csv", FOOTPRINT_LINKS_MAX_AGE_MS)).toString();

1941 const links = new Map();

1942 for (const line of csv.split("\n").slice(1)) {

1943 const [, key, url] = line.split(",");

1944 if (!key || !url) continue;

1945 if (!links.has(key)) links.set(key, []);

1946 links.get(key).push(url.trim());

1947 }

1948 footprintLinks = links;

1949 return links;

1950}

1951

1952// Buildings whose center is inside the box, or null when Microsoft has no data for this

1953// part of the world (then the caller falls back to finding roofs by color).

1954// A missing square inside a covered region just means no buildings were found there.

1955async function loadFootprints(bbox) {

1956 const links = await loadFootprintLinks();

1957 const keys = quadkeysForBbox(bbox);

1958 const covered = keys.some((key) => {

1959 const region = key.slice(0, 6);

1960 for (const known of links.keys()) if (known.startsWith(region)) return true;

1961 return false;

1962 });

1963 if (!covered) return null;

1964

1965 const buildings = [];

1966 const margin = 0.01;

1967 for (const key of keys) {

1968 for (const [i, url] of (links.get(key) ?? []).entries()) {

1969 const gz = await cachedDownload(url, `${key}-${i}.csv.gz`);

1970 for (const line of gunzipSync(gz).toString().split("\n")) {

1971 // Cheap check on the first corner before parsing the whole line.

1972 const at = line.indexOf("[[[");

1973 if (at < 0) continue;

1974 const comma = line.indexOf(",", at);

1975 const lon = Number(line.slice(at + 3, comma));

1976 const lat = Number(line.slice(comma + 1, line.indexOf("]", comma)));

1977

        if (lon < bbox.west - margin || lon > bbox.east + margin || lat < bbox.south - margin || lat > bbox.north + margin) {

1978 continue;

1979 }

1980 const ring = JSON.parse(line).geometry.coordinates[0];

1981 const [cx, cy] = ringCentroid(ring);

1982 if (cx < bbox.west || cx > bbox.east || cy < bbox.south || cy > bbox.north) continue;

1983 buildings.push({ ring, lon: cx, lat: cy, areaM2: ringStats(ring).area });

1984 }

1985 }

1986 }

1987 return buildings;

1988}

1989

1990// How many homes a building stands for: none for a granary or shed, one for a house,

1991// several for a big compound or block.

1992function homesIn(areaM2) {

1993 if (areaM2 < MIN_HOME_M2) return 0;

1994 return Math.max(1, Math.round(areaM2 / HOME_M2));

1995}

1996

1997// Fill a polygon (pixel coordinates) into a mask with an even-odd scanline.

1998function fillPolygon(mask, width, height, points) {

1999 let minY = Infinity;

2000 let maxY = -Infinity;

2001 for (const [, y] of points) {

2002 minY = Math.min(minY, y);

2003 maxY = Math.max(maxY, y);

2004 }

2005 for (let y = Math.max(0, Math.floor(minY)); y <= Math.min(height - 1, Math.ceil(maxY)); y++) {

2006 const yc = y + 0.5;

2007 const xs = [];

2008 for (let i = 0, j = points.length - 1; i < points.length; j = i++) {

2009 const [xi, yi] = points[i];

2010 const [xj, yj] = points[j];

2011 if (yi > yc !== yj > yc) xs.push(xi + ((yc - yi) * (xj - xi)) / (yj - yi));

2012 }

2013 xs.sort((a, b) => a - b);

2014 for (let k = 0; k + 1 < xs.length; k += 2) {

2015 for (let x = Math.max(0, Math.round(xs[k])); x < Math.min(width, Math.round(xs[k + 1])); x++) {

2016 mask[y * width + x] = 1;

2017 }

2018 }

2019 }

2020}

2021

2022// Settlements from footprints: draw every building on a fine grid, join buildings within

2023// FOOTPRINT_LINK_M by closing the gaps, trace each group, and count its homes.

2024function footprintSettlements(bbox, buildings) {

2025 const lat0 = (bbox.south + bbox.north) / 2;

2026 const mPerLon = 111320 * Math.cos((lat0 * Math.PI) / 180);

2027 const mPerLat = 110540;

2028 const pad = FOOTPRINT_LINK_M;

2029 const width = Math.ceil(((bbox.east - bbox.west) * mPerLon + 2 * pad) / FOOTPRINT_RES_M);

2030 const height = Math.ceil(((bbox.north - bbox.south) * mPerLat + 2 * pad) / FOOTPRINT_RES_M);

2031 const toPx = (lon, lat) => [

2032 ((lon - bbox.west) * mPerLon + pad) / FOOTPRINT_RES_M,

2033 ((bbox.north - lat) * mPerLat + pad) / FOOTPRINT_RES_M,

2034 ];

2035 const toLonLat = (x, y) => [

2036 bbox.west + (x * FOOTPRINT_RES_M - pad) / mPerLon,

2037 bbox.north - (y * FOOTPRINT_RES_M - pad) / mPerLat,

2038 ];

2039

2040 const mask = new Uint8Array(width * height);

2041 const centerPixel = [];

2042 for (const b of buildings) {

2043 fillPolygon(mask, width, height, b.ring.map(([lon, lat]) => toPx(lon, lat)));

2044 // A building smaller than one cell still marks its center.

2045 const [x, y] = toPx(b.lon, b.lat).map(Math.floor);

2046 const pixel = y * width + x;

2047 mask[pixel] = 1;

2048 centerPixel.push(pixel);

2049 }

2050

2051 // Close (grow, then shrink) by half the link distance: gaps narrower than the link

2052 // distance fill in, and the outline still hugs the outer buildings.

2053 const radius = Math.ceil(FOOTPRINT_LINK_M / 2 / FOOTPRINT_RES_M);

2054 const scratch = new Uint8Array(mask.length);

2055 const grown = new Uint8Array(mask.length);

2056 const closed = new Uint8Array(mask.length);

2057 rectFilter(mask, width, height, radius, true, scratch, grown);

2058 rectFilter(grown, width, height, radius, false, scratch, closed);

2059

2060 const labels = new Int32Array(mask.length);

2061 const blobs = components(closed, new Uint8Array(mask.length), width, height, { minPixels: 1, labels });

2062 const byLabel = new Map(blobs.map((blob) => [blob.label, { blob, buildings: [] }]));

2063 buildings.forEach((b, i) => byLabel.get(labels[centerPixel[i]])?.buildings.push(b));

2064

2065 const settlements = [];

2066 for (const { blob, buildings: members } of byLabel.values()) {

2067 const homes = members.reduce((sum, b) => sum + homesIn(b.areaM2), 0);

2068 // Only sheds and granaries: nobody lives here.

2069 if (homes === 0) continue;

2070 const simplified = simplify(blob.ring, 1);

2071 if (simplified.length < 3) continue;

2072 const ring = simplified.map(([x, y]) => toLonLat(x, y));

2073 const builtM2 = members.reduce((sum, b) => sum + b.areaM2, 0);

2074 settlements.push({ ring, buildings: members.length, homes, builtM2 });

2075 }

2076 return settlements;

2077}

2078

2079// Footprint certainty. Same three parts and weights as satelliteCertainty, but the last

2080// part is the share of the outline covered by building footprints, which is real data.

2081function footprintCertainty(buildings, compactness, coverage) {

2082 const buildingScore = Math.min(buildings / 8, 1);

2083 const shapeScore = Math.min(compactness / 0.45, 1);

2084 const coverageScore = Math.min(coverage / FULL_COVERAGE, 1);

2085 const score = 100 * (0.45 * buildingScore + 0.3 * shapeScore + 0.25 * coverageScore);

2086 return Math.min(99, Math.round(score));

2087}

2088

2089// --- Color fallback --------------------------------------------------------

2090

2091// Settlements by finding roofs in Esri satellite pixels. Used only where there are no footprints.

2092async function pixelSettlements(bbox) {

2093 // Download and stitch the satellite tiles that cover that box.

2094 const grid = await loadSatellite(bbox);

2095

2096 // Split the image into two masks: any structure, and roof-colored pixels.

2097 // Roof pixels are kept separate so buildings can be counted after roofs are joined.

2098 const { built, roof } = classify(grid.rgb, grid.width, grid.height);

2099 // The raw image is large and unused after classification, so release it.

2100 grid.rgb = null;

2101

2102 // Tiles are bigger than the box. Convert each edge into mosaic pixels

2103 // so the next step can erase the margin outside the request.

2104 const left = lonToPixelX(bbox.west, grid);

2105 const top = latToPixelY(bbox.north, grid);

2106 const right = lonToPixelX(bbox.east, grid);

2107 const bottom = latToPixelY(bbox.south, grid);

2108 // Drop built pixels outside the box, or a village just past the edge would be included.

2109 clipToBbox(built, grid.width, grid.height, left, top, right, bottom);

2110 // Same clip on the roof mask, so building counts stay inside the box too.

2111 clipToBbox(roof, grid.width, grid.height, left, top, right, bottom);

2112 // Stderr only. Confirms the classifier saw something before polygons are built.

2113 console.error(`built pixels ${countOnes(built)}, roof pixels ${countOnes(roof)}`);

2114

2115 // Join nearby roofs into one blob per settlement and remove speckle.

2116 const mask = settlementMask(built, grid.width, grid.height);

2117 const settlements = [];

2118 for (const blob of components(mask, roof, grid.width, grid.height)) {

2119 if (blob.buildings < MIN_ROOF_CLUSTERS) continue;

2120 // Pixel outlines follow every stair-step. Keep only corners that change direction.

2121 const simplified = simplify(blob.ring, SIMPLIFY_PX);

2122 // A ring needs three corners. Fewer than that is a degenerate sliver.

2123 if (simplified.length < 3) continue;

2124 // Pixel corners are in image space. GeoJSON needs longitude, latitude.

2125 const ring = simplified.map(([x, y]) => cornerToLonLat(x, y, grid));

2126 // Share of the blob that looked like a roof. Roads and fields score low here.

2127 const roofShare = blob.pixels === 0 ? 0 : blob.roofPixels / blob.pixels;

2128 settlements.push({ ring, buildings: blob.buildings, roofShare });

2129 }

2130 return settlements;

2131}

2132

2133// Run the whole pipeline on one box and return the GeoJSON FeatureCollection.

2134// Progress lines go to console.error. Throws when tiles or OpenStreetMap can't be reached.

2135export async function detectVillages(bbox) {

2136 // The only input is the bounding box. Everything else is derived from it.

2137 validateBbox(bbox);

2138

2139 // Real building outlines where Microsoft has them; roof colors only as a fallback.

2140 let buildingSource = "footprints";

2141 let found;

2142 let footprints = null;

2143 try {

2144 footprints = await loadFootprints(bbox);

2145 } catch (error) {

2146 console.error(`footprints unavailable: ${error.message || error}`);

2147 }

2148 if (footprints) {

2149 console.error(`footprint buildings ${footprints.length}`);

2150 found = footprintSettlements(bbox, footprints);

2151 } else {

2152 buildingSource = "pixels";

2153 found = await pixelSettlements(bbox);

2154 }

2155 console.error(`settlements ${found.length} (${buildingSource})`);

2156

2157 // OSM places are the only source that can produce certainty 100.

2158 const places = await fetchOsmPlaces(bbox);

2159 console.error(`osm places ${places.length}`);

2160

2161 // One GeoJSON feature per village polygon.

2162 const features = [];

2163 for (const settlement of found) {

2164 const geographic = settlement.ring;

2165 // Signed area tells the winding. GeoJSON exteriors must be counter-clockwise.

2166 const { area2 } = ringStats(geographic);

2167 // Negative area is clockwise, so flip the ring. Shape and area stay the same.

2168 if (area2 < 0) geographic.reverse();

2169 // Compactness is measured after the flip, on the ring that will be written out.

2170 const stats = ringStats(geographic);

2171

2172 // Inside the polygon, or within 150 m, this place confirms the village.

2173 const place = nearestPlace(geographic, places);

2174

2175 let certaintyScore;

2176 let population;

2177 let roofShare;

2178 if (buildingSource === "footprints") {

2179 // Share of the outline covered by building footprints.

2180 roofShare = stats.area === 0 ? 0 : Math.min(1, settlement.builtM2 / stats.area);

2181 certaintyScore = footprintCertainty(settlement.buildings, stats.compactness, roofShare);

2182 population = Math.max(1, Math.round(settlement.homes * PEOPLE_PER_HOME));

2183 } else {

2184 // A stringy web of soil and tracks. Keep it only if OSM says people live there.

2185 const stringy = stats.compactness < MIN_COMPACTNESS;

2186 if (stringy && !place) continue;

2187 roofShare = settlement.roofShare;

2188 certaintyScore = satelliteCertainty(settlement.buildings, stats.compactness, roofShare);

2189 // People from roofs and from settlement area. The larger figure is kept.

2190 // A stringy outline's area is mostly empty ground, so it counts roofs only.

2191 population = estimatePopulation(settlement.buildings, stringy ? 0 : stats.area);

2192 }

2193 // A matched OSM place is 100. Anything else is a satellite guess, never 100.

2194 const certainty = place ? 100 : certaintyScore;

2195 const healthcareUsd = healthcareCost(population);

2196

2197 // Round so coordinates stay readable. Six decimals is about 0.1 m.

2198 const closed = geographic.map(([lon, lat]) => [round6(lon), round6(lat)]);

2199 // GeoJSON polygons are closed rings: the first point is repeated at the end.

2200 closed.push(closed[0]);

2201 features.push({

2202 type: "Feature",

2203 properties: {

2204 // 100 when OSM confirmed the blob, otherwise the satellite score.

2205 certainty,

2206 // Which evidence produced that number.

2207 source: place ? "osm" : "satellite",

2208 // Name and id are present only when a place matched. Otherwise null.

2209 osmName: place ? place.name : null,

2210 osmId: place ? place.id : null,

2211 // People and vaccination cost. Reach cost is filled in after routing.

2212 population,

2213 healthcareUsd,

2214 // Where the buildings came from: "footprints" (Microsoft) or "pixels" (roof colors).

2215 buildingSource,

2216 // Inputs to the population estimate and the certainty, kept for explanations.

2217 buildings: settlement.buildings,

2218 homes: settlement.homes ?? null,

2219 areaM2: Math.round(stats.area),

2220 compactness: Math.round(stats.compactness * 1000) / 1000,

2221 roofShare: Math.round(roofShare * 1000) / 1000,

2222 },

2223 // One outer ring. Holes were dropped when the contour was traced.

2224 geometry: { type: "Polygon", coordinates: [closed] },

2225 });

2226 }

2227

2228 // Price a clinic route for every polygon, keep the cheapest one, and plan the team tour.

2229 const { recommendation, tour } = await recommendAccess(features, bbox);

2230

2231 // Wrap the polygons in the standard GeoJSON collection type.

2232 // `recommendation` is the route section: points, population, healthcare

2233 // money, reach difficulty, and the cost to get there.

2234 // `tour` is one round trip from a clinic through every settlement worth visiting.

2235 return { type: "FeatureCollection", recommendation, tour, features };

2236}

2237

2238// --- Joined clinic tours ---------------------------------------------------

2239

2240// One tour from `clinic` through the settlements of several scans, on one road network that

2241// covers all of them. `scans` is [{ id, bbox, features }]. Stops come back as

2242// { scanId, index } (index into that scan's features), in visiting order.

2243export async function planClinicTour(clinic, scans) {

2244 const features = [];

2245 const owners = [];

2246 for (const scan of scans) {

2247 scan.features.forEach((feature, index) => {

2248 // Copies, so planning never writes into the saved scans.

2249 features.push({ ...feature, properties: { ...feature.properties } });

2250 owners.push({ scanId: scan.id, index });

2251 });

2252 }

2253 if (!features.length) return null;

2254 const graph = buildRoadGraph(await fetchRoadNetwork(scans.map((scan) => scan.bbox), [clinic]));

2255

  console.error(`joined tour: ${scans.length} scans, ${features.length} settlements, ${graph.segments.length} road segments`);

2256 const tour = planTour(features, null, [clinic], graph, { depot: clinic });

2257 if (!tour) return null;

2258 return { ...tour, order: tour.order.map((i) => owners[i]) };

2259}

Every scan runs this file. Help us scan a region!

Scan a region