30:T4f6,# PolyLabel _Modified version of [polylabel](https://github.com/mapbox/polylabel) based on [Frank Brill's enhancements](https://observablehq.com/@fzbrill/copy-of-polylabel). This version tidies up some of the JavaScript code, increases flexibility of input polygons and adds some examples for working with Topo/GeoJSON objects._ Given a polygon comprising in the form `[[x0,y0], [x1,y1], ...]` or a multi-polygon in the form `[[[r0x0,r0y0], [r0x1,r0y1], ...],[[r1x0,r1y0], [r1x1,r1y1], ...], ...]`, find the largest incircle that can be placed without crossing any boundaries. The result is a single array in the form `[cx,cy,r]` in the same units as the polygon coordinates. The precision setting (second parameter, default 1) determines the resolution of candidate regions that might contain the incircle centre. Larger values less precise, but less computationally demanding. Measured in whatever units are used to define vertex coordinates. The centre of this circle, also known as the _pole of inaccessibility_, can be used for positioning labels. The radius can be used to summarise how elongated the shape is, for example circle area / polygon area. Can be easily used in other documents with ```javascript import {polylabel} from "@jwolondon/polylabel" ```31:Tecf,/* ISC License Copyright (c) 2016 Mapbox Permission to use, copy, modify, and/or distribute this software for any purpose with or without fee is hereby granted, provided that the above copyright notice and this permission notice appear in all copies. THE SOFTWARE IS PROVIDED "AS IS" AND ISC DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL ISC BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ function polylabel(poly, precision = 1, debug = false) { let polygon = poly; const depth = arrayDepth(poly); if (depth < 2) { // No polygon provided. return null; } else if (depth == 2) { // Single polygon. polygon = [poly]; } else { // List of polygons. let bestpoint = polylabel(polygon[0], precision, debug); for (let p = 1; p < polygon.length; p++) { const thispoint = polylabel(polygon[p], precision, debug); if (thispoint[2] > bestpoint[2]) { bestpoint = thispoint; } } return bestpoint; } // find the bounding box of the outer ring let minX, minY, maxX, maxY; for (let i = 0; i < polygon[0].length; i++) { const p = polygon[0][i]; if (!i || p[0] < minX) { minX = p[0]; } if (!i || p[1] < minY) { minY = p[1]; } if (!i || p[0] > maxX) { maxX = p[0]; } if (!i || p[1] > maxY) { maxY = p[1]; } } const [width, height] = [maxX - minX, maxY - minY]; const cellSize = Math.min(width, height); let h = cellSize / 2; if (cellSize === 0) { const degeneratePoleOfInaccessibility = [minX, minY]; degeneratePoleOfInaccessibility.distance = 0; return degeneratePoleOfInaccessibility; } // Priority queue of cells in order of their "potential" (max distance to polygon) const cellQueue = new TinyQueue(undefined, (a, b) => b.max - a.max); // Cover polygon with initial cells for (let x = minX; x < maxX; x += cellSize) { for (let y = minY; y < maxY; y += cellSize) { cellQueue.push(new Cell(x + h, y + h, h, polygon)); } } // Take centroid as the first best guess let bestCell = getCentroidCell(polygon); // Special case for rectangular polygons const bboxCell = new Cell(minX + width / 2, minY + height / 2, 0, polygon); if (bboxCell.d > bestCell.d) { bestCell = bboxCell; } let numProbes = cellQueue.length; while (cellQueue.length) { // Pick the most promising cell from the queue const cell = cellQueue.pop(); // Update the best cell if we found a better one if (cell.d > bestCell.d) { bestCell = cell; if (debug) { console.log( "found best %d after %d probes", Math.round(1e4 * cell.d) / 1e4, numProbes ); } } // Do not drill down further if there's no chance of a better solution if (cell.max - bestCell.d <= precision) { continue; } // Split the cell into four cells h = cell.h / 2; cellQueue.push(new Cell(cell.x - h, cell.y - h, h, polygon)); cellQueue.push(new Cell(cell.x + h, cell.y - h, h, polygon)); cellQueue.push(new Cell(cell.x - h, cell.y + h, h, polygon)); cellQueue.push(new Cell(cell.x + h, cell.y + h, h, polygon)); numProbes += 4; } if (debug) { console.log("num probes: " + numProbes); console.log("best distance: " + bestCell.d); console.log("best point: [" + bestCell.x + ", " + bestCell.y + "]"); console.log(polygon); } return [bestCell.x, bestCell.y, bestCell.d]; }32:T8e1,/*ISC License Copyright (c) 2017, Vladimir Agafonkin Permission to use, copy, modify, and/or distribute this software for any purpose with or without fee is hereby granted, provided that the above copyright notice and this permission notice appear in all copies. THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ class TinyQueue { constructor(data = [], compare = (a, b) => (a < b ? -1 : a > b ? 1 : 0)) { this.data = data; this.length = this.data.length; this.compare = compare; if (this.length > 0) { for (let i = (this.length >> 1) - 1; i >= 0; i--) { this._down(i); } } } push(item) { this.data.push(item); this._up(this.length++); } pop() { if (this.length === 0) { return undefined; } const top = this.data[0]; const bottom = this.data.pop(); if (--this.length > 0) { this.data[0] = bottom; this._down(0); } return top; } peek() { return this.data[0]; } _up(pos) { const { data, compare } = this; const item = data[pos]; while (pos > 0) { const parent = (pos - 1) >> 1; const current = data[parent]; if (compare(item, current) >= 0) { break; } data[pos] = current; pos = parent; } data[pos] = item; } _down(pos) { const { data, compare } = this; const halfLength = this.length >> 1; const item = data[pos]; while (pos < halfLength) { let bestChild = (pos << 1) + 1; // initially it is the left child const right = bestChild + 1; if (right < this.length && compare(data[right], data[bestChild]) < 0) { bestChild = right; } if (compare(data[bestChild], item) >= 0) { break; } data[pos] = data[bestChild]; pos = bestChild; } data[pos] = item; } }33:[] 2e:{"title":"PolyLabel","cells":[{"id":0,"value":"$30","pinned":false,"hidden":false,"mode":"md"},{"id":629,"value":"* **polylabel(_coords_, _[precision]_, _[debug]_)** where _coords_ is a list of polygon coordinate pairs or a list of polygon rings making a multipolygon; _precision_ in map units (default `1`), determines how precisely to locate the incircle; _debug_ (default `false`) reports details to console if true.","pinned":false,"hidden":false,"mode":"md"},{"id":6,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":144,"value":"## Tests","pinned":false,"hidden":false,"mode":"md"},{"id":368,"value":"_Precision can be set interactively, but here only applied to the non geographic shapes._","pinned":false,"hidden":false,"mode":"md"},{"id":269,"value":"viewof precision = Inputs.range([1, 100000], {\n label: \"Precision\",\n value: 1,\n step: 1,\n transform: Math.log\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":250,"value":"showIncircle(\n [\n [0, 0],\n [10, 0],\n [10, 10],\n [0, 10]\n ],\n \"Simple convex polygon\"\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":247,"value":"showIncircle(\n [\n [105000, 14000],\n [116000, 9000],\n [121000, 0],\n [110000, 5000],\n [101000, 3000]\n ],\n \"Simple concave polygon\"\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":372,"value":"showIncircle(\n [\n [\n [10, 0],\n [20, 0],\n [15, 30]\n ],\n [\n [21, 0],\n [60, 0],\n [55, 30]\n ],\n [\n [19, 30],\n [50, 30],\n [23, 7]\n ]\n ],\n \"Multi-polygon\"\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":289,"value":"showIncircle(\n [\n [\n [100, 400],\n [350, 350],\n [800, 450],\n [800, 10],\n [100, 100]\n ],\n [\n [200, 300],\n [250, 150],\n [500, 150],\n [600, 350]\n ],\n [\n [350, 270],\n [440, 300],\n [490, 200],\n [300, 180]\n ]\n ],\n \"Multi-polygon with holes and islands\"\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":405,"value":"showGeoIncircle(london, d3.geoTransverseMercator().rotate([2, 0]))","pinned":false,"hidden":false,"mode":"ojs"},{"id":609,"value":"showGeoIncircle(projLondon)","pinned":false,"hidden":false,"mode":"ojs"},{"id":83,"value":"## Appendices","pinned":false,"hidden":false,"mode":"md"},{"id":89,"value":"_[TinyQueue](https://github.com/mourner/tinyqueue) is a small footprint binary heap priority queue:_","pinned":false,"hidden":false,"mode":"md"},{"id":14,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":42,"value":"function Cell(x, y, h, polygon) {\n this.x = x; // cell center x\n this.y = y; // cell center y\n this.h = h; // half the cell size\n this.d = pointToPolygonDist(x, y, polygon); // distance from cell center to polygon\n this.max = this.d + this.h * Math.SQRT2; // max distance to polygon within a cell\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":111,"value":"_Signed distance from a point to a polygon boundary (-ve if outside):_","pinned":false,"hidden":false,"mode":"md"},{"id":43,"value":"function pointToPolygonDist(x, y, polygon) {\n let inside = false;\n let minDistSq = Infinity;\n\n for (let k = 0; k < polygon.length; k++) {\n const ring = polygon[k];\n\n for (let i = 0, len = ring.length, j = len - 1; i < len; j = i++) {\n const [a, b] = [ring[i], ring[j]];\n if (\n a[1] > y !== b[1] > y &&\n x < ((b[0] - a[0]) * (y - a[1])) / (b[1] - a[1]) + a[0]\n ) {\n inside = !inside;\n }\n minDistSq = Math.min(minDistSq, getSegDistSq(x, y, a, b));\n }\n }\n return minDistSq === 0 ? 0 : (inside ? 1 : -1) * Math.sqrt(minDistSq);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":120,"value":"_Cell containing the polygon centroid:_","pinned":false,"hidden":false,"mode":"md"},{"id":45,"value":"function getCentroidCell(polygon) {\n let [area, x, y] = [0, 0, 0];\n const points = polygon[0];\n\n for (let i = 0, len = points.length, j = len - 1; i < len; j = i++) {\n const [a, b] = [points[i], points[j]];\n const f = a[0] * b[1] - b[0] * a[1];\n x += (a[0] + b[0]) * f;\n y += (a[1] + b[1]) * f;\n area += f * 3;\n }\n if (area === 0) {\n return new Cell(points[0][0], points[0][1], 0, polygon);\n }\n return new Cell(x / area, y / area, 0, polygon);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":127,"value":"_Squared distance from a point to a segment:_","pinned":false,"hidden":false,"mode":"md"},{"id":47,"value":"function getSegDistSq(px, py, a, b) {\n let [x, y] = [a[0], a[1]];\n let [dx, dy] = [b[0] - x, b[1] - y];\n\n if (dx !== 0 || dy !== 0) {\n const t = ((px - x) * dx + (py - y) * dy) / (dx * dx + dy * dy);\n if (t > 1) {\n x = b[0];\n y = b[1];\n } else if (t > 0) {\n x += dx * t;\n y += dy * t;\n }\n }\n dx = px - x;\n dy = py - y;\n return dx * dx + dy * dy;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":322,"value":"_Polgyons might be arrays of coordinate pairs, or arrays of arrays of coordinate pairs. Find the depth of the polygon array and convert to a single array of coordinates for some operations:_","pinned":false,"hidden":false,"mode":"md"},{"id":181,"value":"function arrayDepth(v) {\n return Array.isArray(v) ? 1 + Math.max(0, ...v.map(arrayDepth)) : 0;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":317,"value":"function flatCoords(poly) {\n return arrayDepth(poly) <= 2 ? poly : flatCoords(poly.flat());\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":575,"value":"_For testing, use the renderer to show incircles:_","pinned":false,"hidden":false,"mode":"md"},{"id":222,"value":"import { Renderer } from \"@jwolondon/renderer\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":225,"value":"function showIncircle(poly, title = \"\") {\n const incircle = polylabel(poly, precision);\n return new Renderer(300, 200)\n .stroke()\n .text(title, 0, 20)\n .fit(flatCoords(poly), 25)\n .fill(\"rgba(180,90,45,0.2)\")\n .stroke(\"black\")\n .polygon(poly)\n .fill()\n .stroke(\"red\")\n .ellipse(...incircle)\n .stroke()\n .fill(\"red\")\n .points([[incircle[0], incircle[1]]])\n .render();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":596,"value":"function showGeoIncircle(geo, proj) {\n const brownFn = (att) => {\n return `hsla(20, ${20 + Math.random() * 50}%, ${\n 60 + Math.random() * 20\n }%, 0.3)`;\n };\n\n // Draw map polygons\n const r = new Renderer(width / 2, width / 3)\n .stroke()\n .text(`Incircles ${proj ? \"before\" : \"after\"} projection`, 5, 15)\n .stroke(\"black\")\n .setProjection(proj ? proj : d3.geoIdentity())\n .strokeWidth(0.2)\n .geoShape(geo, \"properties.GSSCode\", brownFn);\n\n // Draw incircles.\n const coords = geo.features.forEach((feat) => {\n const incircle = polylabel(feat.geometry.coordinates, 0.001);\n r.fill()\n .stroke(\"red\")\n .strokeWidth(1)\n .ellipse(...incircle)\n .stroke()\n .fill(\"red\")\n .points([[incircle[0], incircle[1]]]);\n });\n\n return r.render();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":552,"value":"_Create a cloned copy of a given geo object with all coordinates projected with the given map projection:_","pinned":false,"hidden":false,"mode":"md"},{"id":530,"value":"function projectGeo(geo, proj) {\n const projectCoords = (geom, proj) => {\n if (arrayDepth(geom) > 2) {\n geom.forEach((g) => projectCoords(g, proj));\n } else {\n geom.forEach((coords, i) => (geom[i] = proj(coords)));\n }\n };\n const newGeo = structuredClone(geo);\n newGeo.features.forEach((feat) =>\n projectCoords(feat.geometry.coordinates, proj)\n );\n return newGeo;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":472,"value":"_London features (in geographic and projected coordinates) used for testing:_","pinned":false,"hidden":false,"mode":"md"},{"id":477,"value":"london = await d3\n .json(\"https://gicentre.github.io/data/geoTutorials/londonBoroughs.json\")\n .then((feats) => topojson.feature(feats, feats.objects.boroughs))","pinned":false,"hidden":false,"mode":"ojs"},{"id":556,"value":"projLondon = projectGeo(london, d3.geoTransverseMercator().rotate([2, 0]))","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}