30:T434,The **pycnophylactic reallocation procedure** is a way to smooth the choropleth into a continuous surface, by allocating each of the inhabitants to one of these control points. *There are two constraints:* 1. all the control points have a non-negative value (there’s no such thing as a place with −10 people); 2. the sum of values all the control points in a given area is equal to the total number of inhabitants in that area (nobody is reallocated to a different area than where they were counted by the census). *On top of this, we need to decide on a [boundary condition](https://en.wikipedia.org/wiki/Boundary_value_problem):*

a) Is the city surrounded by “empty” spaces, in which case the boundary condition is to set every unknown point to the fixed value 0 (this case is called the Dirichlet condition);
b) or is it “similar to the closest known area”, in which case we want the *change* of value to be minimal (this is the Neumann condition).

In this case, we’ll go with Dirichlet, and ignore any inhabitants in the surroundings.31:[{"id":"8e4a66499563fbe0c1285ce239ba70f3c67c81ee51e14d56564076498851f56ca580c15c1c17dc1cac083bd2698b48b5bb8d6e48e021b7cbbc86a016a9da90dd","href":"https://static.observableusercontent.com/files/8e4a66499563fbe0c1285ce239ba70f3c67c81ee51e14d56564076498851f56ca580c15c1c17dc1cac083bd2698b48b5bb8d6e48e021b7cbbc86a016a9da90dd","name":"annarbor.geojson","type":"application/geo+json","lastModified":1637507778647.293,"size":18366}] 2e:{"title":"Ann Arbor pycnophylactic","cells":[{"id":866,"value":"# Ann Arbor pycnophylactic\n","pinned":false,"hidden":false,"mode":"md"},{"id":0,"value":"*#30DayMapChallenge — Day 26. Choropleth.*\n\nReproducing yet another of Waldo Tobler’s experiments in cartography. In this case, we begin with the population distribution in Ann Arbor (Michigan), which we display (as we should) as a choropleth map of population *density*, by dividing the absolute number of people counted in each census tract by the tract’s area.","pinned":false,"hidden":false,"mode":"md"},{"id":9,"value":"{\n await visibility();\n\n const context = DOM.context2d(width, height, 1);\n const path = d3.geoPath(projection, context);\n\n for (let i = 0; i < features.features.length; i++) {\n context.beginPath();\n path(features.features[i]);\n context.strokeStyle = context.fillStyle = color(\n total[i] / classes.get(i).length\n );\n context.fill();\n context.lineWidth = 3;\n context.stroke();\n }\n\n if (showDots) {\n path.pointRadius(1);\n for (let i = 0; i < points.length; i++) {\n if (tracts[i] === undefined) continue;\n context.beginPath();\n path({ type: \"Point\", coordinates: points[i] });\n context.fillStyle = \"#eee\";\n context.fill();\n }\n }\n context.canvas.style.maxWidth = \"400px\";\n return context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":489,"value":"On this map, we’re going to lay out an arbitrary distribution of control points (use the form below to modify and show these control points). It can be a grid, a hex grid, a poisson-distribution… whatever looks best to you!","pinned":false,"hidden":false,"mode":"md"},{"id":195,"value":"viewof gridtype = Inputs.select([\"poisson\", \"hex\", \"square\", \"random\"], {\n value: \"poisson\",\n label: \"point distribution\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":482,"value":"viewof showDots = Inputs.toggle({ label: \"Show dots\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":535,"value":"$30","pinned":false,"hidden":false,"mode":"md"},{"id":580,"value":"## Pycnophylactic allocation algorithm\n\nOur algorithm begins by linking each control point to one of the areas (and *undefined* if it’s outside the considered area).","pinned":false,"hidden":false,"mode":"md"},{"id":33,"value":"tracts = Array.from(points, (p) => {\n for (let i = 0; i < features.features.length; i++) {\n if (d3.geoContains(features.features[i], p)) return i;\n }\n\n // try and repair the few holes that were in the original geojson file\n for (let i = 0; i < features.features.length; i++) {\n if (d3.geoContains(features.features[i], [p[0] + 0.001, p[1]])) return i;\n }\n for (let i = 0; i < features.features.length; i++) {\n if (d3.geoContains(features.features[i], [p[0], p[1] + 0.001])) return i;\n }\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":657,"value":"classes = d3.group(d3.range(points.length), (i) => tracts[i])","pinned":true,"hidden":false,"mode":"ojs"},{"id":591,"value":"Then we create a graph of neighbors, with d3-geo-voronoi:","pinned":false,"hidden":false,"mode":"md"},{"id":599,"value":"voronoi = d3.geoVoronoi(points)","pinned":true,"hidden":false,"mode":"ojs"},{"id":601,"value":"The total population of each census tract is given by the field CY7001 of the GeoJSON:","pinned":false,"hidden":false,"mode":"md"},{"id":610,"value":"total = features.features.map((d) => d.properties.CY7001)","pinned":true,"hidden":false,"mode":"ojs"},{"id":613,"value":"We can now prepopulate an array of values, for each point, starting with this total divided by the number of points that belong to the same area:","pinned":false,"hidden":false,"mode":"md"},{"id":708,"value":"counts = d3.rollup(\n points,\n (v) => v.length,\n (p, i) => tracts[i]\n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":616,"value":"init = Array.from(points, (d, i) => total[tracts[i]] / counts.get(tracts[i]))","pinned":false,"hidden":false,"mode":"ojs"},{"id":629,"value":"At each step of the computation, we’ll reallocate the numbers from neighbor to neighbor, to make the distribution smooth. But we’ll also apply the two constraints. (And, for the Dirichlet boundary condition, we’ll just ignore any point that has a NaN value.)","pinned":false,"hidden":false,"mode":"md"},{"id":634,"value":"function pycnoStep(values0) {\n const values1 = Float64Array.from(values0);\n\n const neighbors = voronoi.delaunay.neighbors;\n\n for (const [tract, I] of classes) {\n if (tract === undefined) continue;\n\n // smooth by making our value closer to the mean of our neighbors\n for (const i of I) {\n values1[i] +=\n (d3.mean(neighbors[i], (j) => values0[j]) - values0[i]) * 0.7;\n }\n\n // enforce sums\n const sum = d3.sum(I, (i) => values1[i]);\n for (const i of I) {\n values1[i] += ((total[tract] - sum) / I.length) * 1;\n }\n\n // clamp each point above 0 and reallocate to others in the same class\n const neg = d3.group(I, (i) => values1[i] < 0);\n if (neg.has(true)) {\n const reallocate =\n -d3.sum(neg.get(true), (i) => values1[i]) / neg.get(false).length;\n for (const i of neg.get(true)) values1[i] = 0;\n for (const i of neg.get(false)) values1[i] += reallocate;\n }\n }\n\n return values1;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":724,"value":"viewof steps = Inputs.range([0, 25], {\n label: \"number of reallocation steps\",\n step: 1,\n value: 10\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":71,"value":"pycnoMap = {\n const context = DOM.context2d(width, height);\n const path = d3.geoPath(projection, context);\n\n const polygons = voronoi.polygons();\n\n update(init);\n\n return Object.assign(context.canvas, { update });\n\n function update(values) {\n context.clearRect(0, 0, width, height);\n\n for (let i = 0; i < points.length; i++) {\n if (tracts[i] === undefined) continue;\n context.beginPath();\n path(polygons.features[i]);\n context.fillStyle = color(values[i]);\n context.fill();\n context.strokeStyle = context.fillStyle;\n context.stroke();\n }\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":846,"value":"md`*Error:* **${error}** people are wrongly allocated — **${d3.format(\".2%\")(\n error / d3.sum(total)\n)}** of the total population (${d3.format(\",\")(d3.sum(total))}).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":736,"value":"// applies the reallocation *n* times, and updates the map\nreallocated = {\n let v = init;\n for (let i = 0; i < steps; i++) {\n v = pycnoStep(v);\n }\n pycnoMap.update(v);\n return v;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":812,"value":"error = Math.round(\n d3.sum(\n d3.rollup(\n reallocated,\n (v) => d3.sum(v),\n (d, i) => tracts[i]\n ),\n ([tract, sum]) => Math.abs(total[tract] - sum)\n )\n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":778,"value":"-------\n*boring zone*","pinned":false,"hidden":false,"mode":"md"},{"id":26,"value":"points = Array.from(pick2d(width, height, 10000, gridtype), projection.invert)","pinned":true,"hidden":false,"mode":"ojs"},{"id":673,"value":"color = d3\n .scaleSequential(d3.interpolateCividis)\n .domain(d3.extent(total.map((d, i) => d / classes.get(i).length)))","pinned":true,"hidden":false,"mode":"ojs"},{"id":8,"value":"d3 = require(\"d3@7\", \"d3-geo-voronoi@1.6\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":34,"value":"features = FileAttachment(\"annarbor.geojson\").json()","pinned":true,"hidden":false,"mode":"ojs"},{"id":15,"value":"projection = d3\n .geoMercator()\n .fitExtent([[20, 20], [width - 20, height - 20]], features)","pinned":true,"hidden":false,"mode":"ojs"},{"id":18,"value":"height = 600","pinned":true,"hidden":false,"mode":"ojs"},{"id":24,"value":"import { pick2d } from \"@fil/2d-point-distributions\"","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W31","stdlib":"1","titling":"auto"}