30:T6dc,history = { const points = data.sPoints.map((pt) => pt.slice(0, 2)).flat(); const delaunay = new d3.Delaunay(points); const voronoi = delaunay.voronoi([0, 0, data.w, data.h]); const omega = 0.5; // Site to centroid projection // const clipRegion = bounds(1.17, 30); const clipRegion = bounds(input.alpha, input.dilation); const history = []; let sites = []; for (let i = 0; i < points.length; i += 2) { const cell = clip(clipRegion, [voronoi.cellPolygon(i >> 1)])[0]; if (cell !== null) { sites.push([points[i], points[i + 1]]); } } history.push(sites); for (let it = 0; it < data.numIterations; it++) { sites = []; // for (let i = 0; i < points.length; i += 2) { // const cell = clip(clipRegion, [voronoi.cellPolygon(i >> 1)])[0]; // if (cell) { // const x0 = points[i]; // const y0 = points[i + 1]; // const [x1, y1] = d3.polygonCentroid(cell); // const xNext = x0 + (x1 - x0) * omega; // const yNext = y0 + (y1 - y0) * omega; // points[i] = xNext; // points[i + 1] = yNext; // sites.push([xNext, yNext]); // } else { // sites.push([points[i], points[i + 1]]); // } // } let i = 0; for (const poly of voronoi.cellPolygons()) { clip(clipRegion, [poly]).forEach((cell) => { const x0 = points[i * 2]; const y0 = points[i * 2 + 1]; const [x1, y1] = d3.polygonCentroid(cell); const xNext = x0 + (x1 - x0) * omega; const yNext = y0 + (y1 - y0) * omega; points[i * 2] = xNext; points[i * 2 + 1] = yNext; sites.push([xNext, yNext]); }); i++; } history.push(sites); voronoi.update(); } return history; }31:T6c4,async function dataObject(dataFilename, options) { const defaults = { xCol: 0, yCol: 1, clrCol: null, geoFilename: null, featureName: null, geoId: null, paletteName: null, numIterations: 100 }; options = Object.assign({}, defaults, options); const points = await d3.csv(dataFilename, d3.autoType); const [x, y, clrId] = [ Object.keys(points[0])[options.xCol], Object.keys(points[0])[options.yCol], Object.keys(points[0])[options.clrCol] ]; const [xMin, xMax] = d3.extent(points.map((o) => o[x])); const [yMin, yMax] = d3.extent(points.map((o) => o[y])); const maxSide = 400; let w = width / 2; let h = (w * (yMax - yMin)) / (xMax - xMin); const longestSide = Math.max(w, h); w = (maxSide * w) / longestSide; h = (maxSide * h) / longestSide; const resolve = (obj, expr) => expr.split(".").reduce((o, i) => o[i], obj); let geo, getClr; if (options.geoFilename && options.topojsonFeature) { geo = await d3 .json(options.geoFilename) .then((feats) => topojson.feature(feats, resolve(feats.objects, options.topojsonFeature)) ); } if (options.paletteName && clrId) { const palette = await d3 .csv(options.paletteName, d3.autoType) .then((obs) => Object.fromEntries(obs.map((o) => [o.id, o.colour]))); getClr = (clr) => palette[clr]; } else { getClr = (clr) => `hsl(20, ${20 + rand.next() * 50}%, ${60 + rand.next() * 20}%)`; } return { w: w, h: h, geo: geo, geoId: options.geoId, points: points, sPoints: toScreen( points.map((o) => [o[x], o[y], o[clrId]]), w, h, 30 ), clr: getClr, numIterations: options.numIterations }; }32:T40f,function areaMeasures(siteHistory) { const clipRegion = bounds(116, 40); const stdev = (xs) => { let total = xs.reduce((a, b) => a + b, 0); let mean = total / xs.length; return Math.sqrt( xs.map((x) => (x - mean) * (x - mean)).reduce((a, b) => a + b, 0) / xs.length ); }; const areaStats = siteHistory.map((hState) => { const voronoi = d3.Delaunay.from(hState).voronoi([0, 0, data.w, data.h]); const areas = []; for (const poly of voronoi.cellPolygons()) { const clipped = clip(clipRegion, [poly])[0]; areas.push(Math.abs(d3.polygonArea(clipped))); } const maxArea = Math.max(...areas); const minArea = Math.min(...areas); return { stdev: stdev(areas), max: maxArea, min: minArea, range: maxArea - minArea }; }); const maxVal = Math.max(...areaStats.map((o) => o.max)); return areaStats.map((o) => ({ max: o.max / maxVal, min: o.min / maxVal, range: (o.max - o.min) / maxVal, stdev: o.stdev / maxVal })); }33:T422,function clip( clipRings, subjectRings, clipType = ClipperLib.ClipType.ctIntersection ) { const toPath = (coords) => coords.map(([x, y]) => Object.fromEntries([ ["X", x], ["Y", y] ]) ); const fromPath = (path) => path.map((obj) => Object.entries(obj).map((o) => o[1])); const clipPaths = clipRings.map((ring) => toPath(ring)); const subjectPaths = subjectRings.map((ring) => toPath(ring)); const clipper = new ClipperLib.Clipper(); const scale = 100; ClipperLib.JS.ScaleUpPaths(clipPaths, scale); ClipperLib.JS.ScaleUpPaths(subjectPaths, scale); clipper.AddPaths(clipPaths, ClipperLib.PolyType.ptClip, true); clipper.AddPaths(subjectPaths, ClipperLib.PolyType.ptSubject, true); let result = new ClipperLib.Paths(); const subjectFillType = ClipperLib.PolyFillType.pftNonZero; const clipFillType = ClipperLib.PolyFillType.pftNonZero; clipper.Execute(clipType, result, subjectFillType, clipFillType); ClipperLib.JS.ScaleDownPaths(result, scale); return result.map(fromPath); }34:[] 2e:{"title":"AVT: Equalising Regions","cells":[{"id":0,"value":"# AVT: Equalising Regions\n\nNow that we are able to transform a given point distribution into a set of bounded regions approximating their layout, we can consider how to equalise their density within those bounds.","pinned":false,"hidden":false,"mode":"md"},{"id":1433,"value":"viewof dataset = Inputs.radio(\n [\"London\", \"UK Regions\", \"Scotland LAs\", \"UK LA Districts\"],\n {\n value: \"London\",\n label: \"Dataset\"\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":135,"value":"{\n const brownFn = (att) => {\n return `hsla(20, ${20 + Math.random() * 50}%, ${\n 60 + Math.random() * 20\n }%, 0.5)`;\n };\n\n if (data.geo) {\n return new Renderer(data.w, data.h, \"svg\")\n .strokeWidth(0.2)\n .geoShape(data.geo, data.geoId, (att) => data.clr(att.split(\" \")[0]))\n .render();\n } else {\n return new Renderer(width, 50, \"svg\")\n .stroke()\n .textSize(12)\n .textBaseline(\"top\")\n .text(\"No topoJSON file provided.\", 0, 0)\n .render();\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":164,"value":"And the clipped and dilated Voronoi representation:","pinned":false,"hidden":false,"mode":"md"},{"id":1806,"value":"viewof input = Object.assign(\n Inputs.form({\n alpha: Inputs.range([0.01, 6], {\n label: \"alpha\",\n value: 5,\n step: 0.01\n }),\n dilation: Inputs.range([0, 60], {\n value: 18,\n label: \"Dilation\",\n step: 1\n }),\n preserveIslands: Inputs.toggle({ label: \"Preserve islands\", value: false })\n }),\n {\n oninput: (event) => {\n return event.isTrusted && event.stopImmediatePropagation();\n },\n onchange: (event) => event.currentTarget.dispatchEvent(new Event(\"input\"))\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":214,"value":"function bounds(a, dil) {\n const delaunay = d3.Delaunay.from(data.sPoints);\n const avR = avRadius(delaunay);\n const voronoi = delaunay.voronoi([0, 0, data.w, data.h]);\n const hull = alphaHull(delaunay, a * avR);\n // const rings = [];\n\n const k = 10; // Radius of unconnected hexagons in geo coordinates (m)\n const rings = input.preserveIslands\n ? unconnected(delaunay, a * avR).map(([x, y]) =>\n dilate(\n [\n [x + 0.866 * k, y + 0.5 * k],\n [x + 0.866 * k, y - 0.5 * k],\n [x, y - k],\n [x - 0.866 * k, y - 0.5 * k],\n [x - 0.866 * k, y + 0.5 * k],\n [x, y + k]\n ],\n input.dilation\n )\n )\n : [];\n\n hull.forEach((ring) => {\n rings.push(\n dilate(\n ring.map((v) => [delaunay.points[2 * v], delaunay.points[2 * v + 1]]),\n dil,\n ClipperLib.JoinType.jtSquare\n )\n );\n });\n return rings;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1198,"value":"function avRadius(delaunay) {\n const circumRadius = ([[x1, y1], [x2, y2], [x3, y3]]) => {\n const a2 = x1 - x2;\n const a3 = x1 - x3;\n const b2 = y1 - y2;\n const b3 = y1 - y3;\n const d1 = x1 * x1 + y1 * y1;\n const d2 = d1 - x2 * x2 - y2 * y2;\n const d3 = d1 - x3 * x3 - y3 * y3;\n const ab = (a3 * b2 - a2 * b3) * 2;\n const xa = (b2 * d3 - b3 * d2) / ab - x1;\n const ya = (a3 * d2 - a2 * d3) / ab - y1;\n if (isNaN(xa) || isNaN(ya)) {\n return 0;\n }\n return Math.sqrt(xa * xa + ya * ya);\n };\n const { points, triangles } = delaunay;\n\n const pts = [];\n triangles.forEach((i) => pts.push([points[i * 2], points[i * 2 + 1]]));\n const radii = Array.from({ length: pts.length / 3 }, (v, i) =>\n pts.slice(i * 3, i * 3 + 3)\n ).map(circumRadius);\n\n return radii.reduce((a, b) => a + b, 0) / radii.length;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":49,"value":"{\n const voronoi = d3.Delaunay.from(data.sPoints).voronoi([\n 0,\n 0,\n data.w,\n data.h\n ]);\n const clipRegion = bounds(input.alpha, input.dilation);\n const r = new Renderer(data.w, data.h, \"svg\");\n\n [...voronoi.cellPolygons()].forEach((poly, i) => {\n const clipped = clip(clipRegion, [poly])[0];\n clipped ? r.fill(data.clr(data.sPoints[i][2])).polygon(clipped) : r;\n });\n r.fill(\"black\").stroke().points(data.sPoints);\n return r.render();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":172,"value":"From this start point we can apply [Lloyd's relaxation](https://observablehq.com/@mattdzugan/lloyds-algorithm-regular-tilings) to produce more uniformly spaced regions. Unlike some applications of Lloyd's algorithm, the aim here is also to maintain as much of the original geospatial relationships observed in the initial placement of points and regions as we can.","pinned":false,"hidden":false,"mode":"md"},{"id":225,"value":"{\n resetButton; // Force recalculation when resetButton emits an input event\n const points = data.sPoints.map((pt) => pt.slice(0, 2)).flat();\n const delaunay = new d3.Delaunay(points);\n const voronoi = delaunay.voronoi([0, 0, data.w, data.h]);\n const omega = 0.5; // Site to centroid projection;\n const clipRegion = bounds(input.alpha, input.dilation);\n const r = new Renderer(data.w, data.h, \"svg\");\n\n for (let it = 0; it < data.numIterations; it++) {\n rand.reset();\n [...voronoi.cellPolygons()].forEach((poly, i) => {\n clip(clipRegion, [poly]).forEach((cell) => {\n r.fill(data.clr(data.sPoints[i][2])).polygon(cell);\n const x0 = points[i * 2];\n const y0 = points[i * 2 + 1];\n const [x1, y1] = d3.polygonCentroid(cell);\n points[i * 2] = x0 + (x1 - x0) * omega;\n points[i * 2 + 1] = y0 + (y1 - y0) * omega;\n r.fill(\"black\").points([[x1, y1]]);\n });\n });\n yield r.render();\n voronoi.update();\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":384,"value":"viewof resetButton = Inputs.button(\"Reset\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":399,"value":"We can examine the transformation of the region sites more carefully to see what happens on each iteration of the relaxation. Below we calculate and store the positions of all sites during each Lloyd's iteration.","pinned":false,"hidden":false,"mode":"md"},{"id":403,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":684,"value":"Once stored, we can then play back through the iterations.","pinned":false,"hidden":false,"mode":"md"},{"id":423,"value":"{\n const voronoi = d3.Delaunay.from(history[iteration]).voronoi([\n 0,\n 0,\n data.w,\n data.h\n ]);\n // const clipRegion = bounds(1.17, 30);\n const clipRegion = bounds(input.alpha, input.dilation);\n\n const r = new Renderer(data.w, data.h, \"svg\");\n rand.reset();\n [...voronoi.cellPolygons()].forEach((poly, i) => {\n clip(clipRegion, [poly]).forEach((cell) => {\n if (data.sPoints[i]) {\n r.fill(data.clr(data.sPoints[i][2])).polygon(cell);\n }\n });\n });\n r.fill(\"black\").stroke().points(history[iteration]);\n return r.render();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":376,"value":"viewof iteration = Scrubber(\n history.map((_, i) => i),\n { autoplay: false, alternate: true, delay: 50 }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":802,"value":"With this pre-computed set of locations, it becomes easy to show how various measures of geospatial/topological character change with each iteration of the relaxation. Below we show maximum region size, minimum region size, standard deviation of region size and average displacement at each iterations. Most of the change in properties occur with the first few iterations, so we use a square root scale on the x-axis.","pinned":false,"hidden":false,"mode":"md"},{"id":832,"value":"{\n const padL = 10;\n const padR = 320;\n const padB = 20;\n const h = 200;\n const xSc = d3\n .scaleSqrt()\n .domain([0, data.numIterations])\n .range([padL, width - padR]);\n const ySc = d3\n .scaleLinear()\n .domain([0, 1])\n .range([h - padB, 0]);\n\n return new Renderer(width, h, \"svg\")\n .setTransform(([x, y]) => [xSc(x), ySc(y)])\n .fill(\"rgb(0, 0, 0, 0.05)\")\n .stroke()\n .polygon(\n stats\n .map((o, i) => [i, o.areaMax])\n .concat(stats.map((o, i) => [i, o.areaMin]).reverse())\n )\n .stroke(\"#2ca02c\")\n .line(stats.map((o, i) => [i, o.areaMax]))\n .stroke(\"#1f77b4\")\n .line(stats.map((o, i) => [i, o.areaMin]))\n .stroke(\"#ff7f0e\")\n .line(stats.map((o, i) => [i, o.areaStdev]))\n .stroke(\"#d62728\")\n .line(stats.map((o, i) => [i, o.avDisp]))\n .stroke(\"grey\")\n .line([\n [iteration, 0],\n [iteration, 1]\n ])\n .fill(\"black\")\n .axisBottom(xSc, ySc)\n .render();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":8,"value":"---\n\n## Appendix\n\n### Data","pinned":false,"hidden":false,"mode":"md"},{"id":1415,"value":"Function that takes at least a point CSV file containing positions and at least one attribute of each point to allocate and creates a data object from it containing the scaled point coordinates. By default columns 0 and 1 are assumed to store x and y coordinates respectively, but if not, these can be optionally specified with `xCol` and `yCol` parameters. Can optionally provide topoJSON file and colour palette representing the same region. If a colour palette is provided, the `clrCol` parameter should indicate which column in the point CSV contains the attribute to map to colour.","pinned":false,"hidden":false,"mode":"md"},{"id":1221,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":1421,"value":"From this we can create our actual data objects.","pinned":false,"hidden":false,"mode":"md"},{"id":695,"value":"london = dataObject(\"https://jwolondon.github.io/web/data/londonBoroughs.csv\", {\n geoFilename:\n \"https://gicentre.github.io/data/geoTutorials/londonBoroughs.json\",\n topojsonFeature: \"boroughs\",\n geoId: \"id\",\n paletteName: \"https://jwolondon.github.io/web/data/londonBoroughsPalette.csv\",\n clrCol: 4\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1430,"value":"ukRegions = dataObject(\"https://jwolondon.github.io/web/data/ukRegions.csv\", {\n geoFilename: \"https://jwolondon.github.io/web/data/ukRegions.json\",\n topojsonFeature: \"regions\",\n geoId: \"id\",\n paletteName: \"https://jwolondon.github.io/web/data/ukRegionsPalette.csv\",\n clrCol: 2\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1562,"value":"scotlandLAs = dataObject(\n \"https://jwolondon.github.io/web/data/scotlandLAs.csv\",\n {\n paletteName: \"https://jwolondon.github.io/web/data/scotlandLAsPalette.csv\",\n clrCol: 2,\n numIterations: 200\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1689,"value":"ukLADistricts = dataObject(\n \"https://gicentre.github.io/data/ramp/ukLADistricts2019Centroids.csv\",\n {\n xCol: 1,\n yCol: 2,\n geoFilename: \"https://gicentre.github.io/data/ramp/UKLADistricts2019.json\",\n topojsonFeature: \"UKLADistricts2019\",\n geoId: \"properties.LAD19CD\",\n numIterations: 10\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1436,"value":"data = {\n switch (dataset) {\n case \"London\":\n return london;\n case \"UK Regions\":\n return ukRegions;\n case \"Scotland LAs\":\n return scotlandLAs;\n case \"UK LA Districts\":\n return ukLADistricts;\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1386,"value":"### Seeded Random Numbers\n\nUsed when no palette is provided so that instead we generate a persistent, but randomly allocate set of region colours.","pinned":false,"hidden":false,"mode":"md"},{"id":1320,"value":"function RandGen(seed = \"apples\") {\n this.next = randGenerator(seed);\n this.reset = function () {\n this.next = randGenerator(seed);\n };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1381,"value":"rand = new RandGen()","pinned":false,"hidden":false,"mode":"ojs"},{"id":587,"value":"### Summary Statistics","pinned":false,"hidden":false,"mode":"md"},{"id":595,"value":"Area statistics are calculated together as they involve the comparatively expensive construction of Voronoi regions for each set of locations. All stats use a common scale in range [0,1] to aid visual comparison.","pinned":false,"hidden":false,"mode":"md"},{"id":437,"value":"stats = {\n const areaStats = areaMeasures(history);\n const avDisps = avDisplacement(history);\n const stats = areaStats.map((areaStat, i) => ({\n iteration: i,\n areaStdev: areaStat.stdev,\n areaMax: areaStat.max,\n areaMin: areaStat.min,\n areaRange: areaStat.range,\n avDisp: avDisps[i]\n }));\n\n const maxDisp = Math.max(...stats.map((o) => o.avDisp));\n return stats.map((o) => ({ ...o, avDisp: o.avDisp / maxDisp }));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":432,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":591,"value":"function avDisplacement(siteHistory) {\n const distSq = ([x0, y0], [x1, y1]) =>\n (x0 - x1) * (x0 - x1) + (y0 - y1) * (y0 - y1);\n return siteHistory.map((hState) => {\n let totalDistSq = 0;\n for (let i = 0; i < hState.length; i++) {\n totalDistSq += distSq(hState[0], hState[i]);\n }\n return totalDistSq / hState.length;\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":50,"value":"### Point Transformations\n\nScale input points so they fit to screen space (defined by `w` and `h` with optional `border`, all in pixel values).","pinned":false,"hidden":false,"mode":"md"},{"id":56,"value":"function toScreen(pts, w, h, border = 20) {\n const xs = pts.map((pt) => pt[0]);\n const ys = pts.map((pt) => pt[1]);\n const xMin = Math.min(...xs);\n const yMin = Math.min(...ys);\n const xMax = Math.max(...xs);\n const yMax = Math.max(...ys);\n\n const sf = Math.min(\n (w - 2 * border) / (xMax - xMin),\n (h - 2 * border) / (yMax - yMin)\n );\n return pts.map((pt) =>\n [border + sf * (pt[0] - xMin), h - border - sf * (pt[1] - yMin)].concat(\n pt.slice(2)\n )\n );\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":42,"value":"### Clipping and Dilation\n\nFor constraining region by the alpha shape of the point distribution with some dilation of the outer boundary.","pinned":false,"hidden":false,"mode":"md"},{"id":46,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":76,"value":"function dilate(poly, offset = 0, jt = ClipperLib.JoinType.jtMiter) {\n const toPath = (coords) =>\n coords.map(([x, y]) =>\n Object.fromEntries([\n [\"X\", x],\n [\"Y\", y]\n ])\n );\n const fromPath = (path) =>\n path.map((obj) => Object.entries(obj).map((o) => o[1]));\n\n let path = toPath(poly);\n\n // Ensure we have no very close vertices before dilating\n const cleanDelta = 0.1;\n path = ClipperLib.Clipper.CleanPolygon(path, cleanDelta);\n\n // Create the offset object.\n const miterLimit = 2;\n const arcTolerance = 0.25;\n const co = new ClipperLib.ClipperOffset(miterLimit, arcTolerance);\n co.AddPath(path, jt, ClipperLib.EndType.etClosedPolygon);\n\n const result = new ClipperLib.Path();\n co.Execute(result, offset);\n return fromPath(result[0]);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":26,"value":"### Imports","pinned":false,"hidden":false,"mode":"md"},{"id":32,"value":"import { Renderer } from \"@jwolondon/renderer\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":35,"value":"import { alphaHull, unconnected } from \"@jwolondon/alpha-shapes\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1297,"value":"import { randGenerator } from \"@jwolondon/seeded-random-numbers\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":373,"value":"import { Scrubber } from \"@mbostock/scrubber\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":132,"value":"topojson = require(\"topojson-client@3\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":39,"value":"ClipperLib = require(\"https://makio135.gitlab.io/js/clipper.min.js\").catch(\n () => window.ClipperLib\n)","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}