30:Tb6c,md`# Centerlines from a GeoJSON Polygon with \`d3-delaunay\`

Work in progress: It works, but this notebook is somehwhat of a mess and the algorithm requires performance optimizations. Suggestions are welcome.

An attempt to convert geoJSON polygons to geoJSON LineStrings. For example for polygons resulting from [this little Mapbox tool to convert SVG to geoJSON](https://labs.mapbox.com/svg-to-geojson/) or when you have polygons of streets or river that are better off as lines. Largely based on [Noah Veltman's centerline labeling](https://observablehq.com/@veltman/centerline-labeling) notebook but using the new [\`d3-delaunay \`](https://github.com/d3/d3-delaunay) which is faster but complicates some other things in the process, and [ngraph](https://github.com/anvaka/ngraph.path) instead of node-dijkstra. [See also](https://observablehq.com/d/9ff706ed30c2aea8), the same thing but with the older d3-voronoi, which I won't continue working on. ## To Do & Ideas ### Performance 1. More efficient traversal, use distance 'inside polygon' to other nodes. * Typed arrays * Use (async?) generators * Test more ngraph pathfinders and heuristics * Put pruning operation behind a toggle * Make parts of the notebook wait on eachother? ### Nice features * Animate longest path finding process ([inspiration](https://observablehq.com/@mbostock/the-delaunays-dual)) * Select end nodes yourself by clicking them * Drag-n-drop geoJSON polygon (or FeatureCollection of polygons) into Observable notebook and have it auto-generate the centerline. * Add some polygons of meandering rivers ### Later * Rough estimate for perimeter points needed for a polygon (sinuosity?) * Make centerline all polygons seperately in a MultiPolygon instead of only the first like now * Or use convex hull to make one centerline for all polygons in MultiPolygon * Voronoi discrete skeleton ` // ### Alternatives // * A simpler(?) '[collapse](http://resources.esri.com/help/9.3/arcgisengine/java/gp_toolref/data_management_tools/collapse_dual_lines_to_centerline_data_management_.htm)' algorithm for 'simple' polygons with perpendicular lines // Why Javascript? Because it runs in the browser which means everyone with an internet connection can run it like that! // Sure, QGIS, GRASS, GDAL or ArcGIS have powerful algorithms to generate centerlines, but those are heavy, a lot harder to use (for me at least), to say nothing of the price of ArcGIS. I much prefer using [Mapshaper](https://mapshaper.org/), [turf.js](http://turfjs.org/) and [d3-geo](https://github.com/d3/d3-geo) and similiar utilities to handle geodata, generating centerlines is one that I've found missing. 31:T40d,dragtop = { const svg = d3.create("svg") .attr("viewBox", [0, 0, width, height]); const g = svg.append("g") .attr("cursor", "grab"); const origalPolygon = g.append("path") .attr("d", polygonPath) .attr("stroke", "plum") .attr("stroke-width", 0.3) .attr("fill", "none") .append("title") .text("The original polygon") g.append("path") .attr("d", path(simplifiedCenterlineGeoJSON.geometry)) .attr("stroke", lineLegend.get("Path most distant nodes").stroke) .attr("stroke-width", 0.5) .attr("fill", "none") .append("title") .text("Centerline") const bounds = g.append("path") .attr("fill", "none") .attr("stroke", "dimgray") .attr("stroke-width", 1.75) .attr("d", voronoi.renderBounds()); svg.call(d3.zoom() .extent([[0, 0], [width, height]]) .scaleExtent([1, 8]) .on("zoom", zoomed)); function zoomed(event) { g.attr("transform", event.transform); } return svg.node(); }32:T616,legendSvg = { const svg = d3.create("svg") .attr("viewBox", [0, 0, width, 80]); const points = svg.append("g") .attr("class", "legendGroup") .selectAll("g") .data(pointLegend) .join("g") .attr("transform", (d, i) => { // const width d.name.length * 10 // needs the length of the previous text // .getComputedTextLength()? return `translate(${10 + i * 175},5)` }) points.append("circle") .attr("cy", 20) .attr("r", 7) .attr("fill", d => d[1].fill) points.append("text") .attr("x", 11) .attr("y", 21) .attr("font-size", 13) .attr("alignment-baseline", "middle") .attr("text-anchor", "start") .text(d => d[0]) const lines = svg.append("g") .attr("class", "legendGroup") .attr("transform", "translate(0,60)") .selectAll("g") .data(lineLegend) .join("g") .attr("transform", (d, i) => { // const width d.name.length * 10 // needs the length of the previous text // .getComputedTextLength()? return `translate(${10 + i * 170},5)` }) lines.append("line") .attr("stroke", d => d[1].stroke) .attr("stroke-width", "5") .attr("fill", d => d[1].fill) .attr("x1", -10) .attr("y1", 0) .attr("x2", 17.5) .attr("y2", 0); lines.append("text") .attr("x", 20) .attr("y", 0) .attr("font-size", 13) .attr("alignment-baseline", "middle") .attr("text-anchor", "start") .text(d => d[0]) return svg.node() }33:T806,viewPointsAndVoronoi = { const svg = d3.create("svg") .attr("viewBox", [0, 0, width, height]); const g = svg.append("g") .attr("cursor", "grab"); const bounds = g.append("path") .attr("fill", "none") .attr("stroke", "dimgray") .attr("stroke-width", 1.75) .attr("d", voronoiBounds); // const delaunayMesh = g.append("path") // .attr("fill", "none") // .attr("stroke", lineLegend.get("Voronoi line").stroke) // .attr("opacity", 0.2) // .attr("mix-blend-mode", "lighten") // .attr("stroke-width", 0.2) // .attr("d", delaunayPath); const voronoiMesh = g.append("path") .attr("fill", "none") .attr("stroke", lineLegend.get("Voronoi line").stroke) .attr("opacity", 0.7) .attr("stroke-width", 0.3) .attr("d", voronoiPath); const origalPolygon = g.append("path") .attr("d", polygonPath) .attr("stroke", "#444") .attr("stroke-width", 0.2) .attr("fill", "none"); const allNodes = g.append("g") .selectAll("circle") .data(allCellPolygons.flat()) .join("circle") .attr("cx", d => d[0]) .attr("cy", d => d[1]) .attr("r", 0.45) .attr("fill", pointLegend.get("Voronoi node").fill); const polygonPoint = g.append("g") .selectAll("circle") .data(pointsOnPolygon) .join("circle") .attr("cx", d => d[0]) .attr("cy", d => d[1]) .attr("r", 0.55) .attr("fill", pointLegend.get("Polygon perimeter node").fill) .append("title") .text((d,i) => i.toString()); const insideNodesSvg = g.append("g") .selectAll("circle") .data(innerNodes) .join("circle") .attr("cx", d => d[0]) .attr("cy", d => d[1]) .attr("r", 0.45) .attr("fill", pointLegend.get("Inner node").fill); svg.call(d3.zoom() .extent([[0, 0], [width, height]]) .scaleExtent([1, 10]) .on("zoom", zoomed)); function zoomed(event) { g.attr("transform", event.transform); } return svg.node(); }34:T98a,viewClippedVoronoi = { const svg = d3.create("svg") .attr("viewBox", [0, 0, width, height]); const g = svg.append("g") .attr("cursor", "grab"); const bounds = g.append("path") .attr("fill", "none") .attr("stroke", "dimgray") .attr("stroke-width", 1.75) .attr("d", voronoiBounds); const voronoiMesh = g.append("path") .attr("fill", "none") .attr("stroke", "#ccc") .attr("opacity", 0.6) .attr("stroke-width", 0.5) .attr("d", voronoi.render()); const originalPolygon = g.append("path") .attr("d", polygonPath) .attr("stroke", "#444") .attr("opacity", 0.3) .attr("stroke-width", 0.2) .attr("fill", "none"); const polygonPoint = g.append("g") .selectAll("circle") .data(pointsOnPolygon) .join("circle") .attr("cx", d => d[0]) .attr("cy", d => d[1]) .attr("r", 0.25) .attr("fill", pointLegend.get("Polygon perimeter node").fill) .append("title") .text((d,i) => i.toString()); g.append("g").selectAll("line") .data(uniqueEdges) .enter().append("line") .attr("stroke", "mediumseagreen") .attr("stroke-width", 0.3) .attr("x1", e => e[0][0]) .attr("y1", e => e[0][1]) .attr("x2", e => e[1][0]) .attr("y2", e => e[1][1]); const allNodes = g.append("g") .selectAll("circle") .data(allCellPolygons.flat()) .join("circle") .attr("cx", d => d[0]) .attr("cy", d => d[1]) .attr("r", 0.35) .attr("fill", pointLegend.get("Voronoi node").fill); const insideCircles = g.append("g") .selectAll("circle") .data(innerNodes) .join("circle") .attr("cx", d => d[0]) .attr("cy", d => d[1]) .attr("r", 0.2) .attr("fill", pointLegend.get("Inner node").fill) .append("title") .text((d) => d.id); const clippedCircles = g.append("g") .selectAll("circle") .data(clippedNodes) .join("circle") .attr("cursor", "pointer") .attr("cx", d => d[0]) .attr("cy", d => d[1]) .attr("r", 0.35) .attr("fill", pointLegend.get("Clipped node").fill) .append("title") .text((d) => d.id); svg.call(d3.zoom() .extent([[0, 0], [width, height]]) .scaleExtent([1, 10]) .on("zoom", zoomed)); function zoomed(event) { g.attr("transform", event.transform); } return svg.node(); }35:T933,viewPrunedNodes = { const svg = d3.create("svg") .attr("viewBox", [0, 0, width, height]); const g = svg.append("g") .attr("cursor", "grab"); const bounds = g.append("path") .attr("fill", "none") .attr("stroke", "dimgray") .attr("stroke-width", 1.75) .attr("d", voronoiBounds); const voronoiMesh = g.append("path") .attr("fill", "none") .attr("stroke", "#ccc") .attr("opacity", 0.6) .attr("stroke-width", 0.5) .attr("d", voronoi.render()); const originalPolygon = g.append("path") .attr("d", polygonPath) .attr("stroke", "#444") .attr("opacity", 0.3) .attr("stroke-width", 0.2) .attr("fill", "none"); const polygonPoint = g.append("g") .selectAll("circle") .data(pointsOnPolygon) .join("circle") .attr("cx", d => d[0]) .attr("cy", d => d[1]) .attr("r", 0.25) .attr("fill", pointLegend.get("Polygon perimeter node").fill) .append("title") .text((d,i) => i.toString()); g.append("g").selectAll("line") .data(uniqueEdges) .enter().append("line") .attr("stroke", "mediumseagreen") .attr("stroke-width", 0.3) .attr("x1", e => e[0][0]) .attr("y1", e => e[0][1]) .attr("x2", e => e[1][0]) .attr("y2", e => e[1][1]); const allNodes = g.append("g") .selectAll("circle") .data(allCellPolygons.flat()) .join("circle") .attr("cx", d => d[0]) .attr("cy", d => d[1]) .attr("r", 0.35) .attr("fill", pointLegend.get("Voronoi node").fill); const insideCircles = g.append("g") .selectAll("circle") .data(innerNodes) .join("circle") .attr("cx", d => d[0]) .attr("cy", d => d[1]) .attr("r", 0.2) .attr("fill", pointLegend.get("Inner node").fill) .append("title") .text((d) => d.id); const prunedCircles = g.append("g") .selectAll("circle") .data(prunedNodes) .join("circle") .attr("cursor", "pointer") .attr("cx", d => d[0]) .attr("cy", d => d[1]) .attr("r", 0.9) .attr("fill", "hotpink") svg.call(d3.zoom() .extent([[0, 0], [width, height]]) .scaleExtent([1, 10]) .on("zoom", zoomed)); function zoomed(event) { g.attr("transform", event.transform); } return svg.node(); }36:T793,{ const svg = d3.create("svg") .attr("viewBox", [0, 0, width, height]); const g = svg.append("g") .attr("cursor", "grab"); const mesh = g.append("path") .attr("fill", "none") .attr("stroke", lineLegend.get("Voronoi line").stroke) .attr("opacity", 0.15) .attr("stroke-width", 0.5) .attr("d", voronoiPath); const bounds = g.append("path") .attr("fill", "none") .attr("stroke", "dimgray") .attr("stroke-width", 1.25) .attr("d", voronoiBounds); const origalPolygon = g.append("path") .attr("d", polygonPath) .attr("stroke", "#444") .attr("stroke-width", 0.2) .attr("fill", "none"); g.append("g").selectAll("line") .data(uniqueEdges) .enter().append("line") .attr("stroke", "mediumseagreen") .attr("stroke-width", 0.3) .attr("x1", e => e[0][0]) .attr("y1", e => e[0][1]) .attr("x2", e => e[1][0]) .attr("y2", e => e[1][1]); if (nodesForPathfinding === "pruned") { const prunedCircles = g.append("g") .selectAll("circle") .data(prunedNodes) .join("circle") .attr("cursor", "pointer") .attr("cx", d => d[0]) .attr("cy", d => d[1]) .attr("r", 0.6) .attr("fill", "hotpink") } const longestLine = (strategy === "longest") ? line(longestPath ) : line(mostDistantNodesPath ); const longestLineStroke = (strategy === "longest") ? lineLegend.get("Longest path").stroke : lineLegend.get("Path most distant nodes").stroke; g.append("path") .attr("d", longestLine) .attr("stroke", longestLineStroke) .attr("stroke-width", 1.25) .attr("fill", "none"); svg.call(d3.zoom() .extent([[0, 0], [width, height]]) .scaleExtent([1, 8]) .on("zoom", zoomed)); function zoomed(event) { g.attr("transform", event.transform); } return svg.node(); }37:T45c,dragfinal = { const svg = d3.create("svg") .attr("viewBox", [0, 0, width, height]); const g = svg.append("g") .attr("cursor", "grab"); const origalPolygon = g.append("path") .attr("d", polygonPath) .attr("stroke", "plum") .attr("stroke-width", 0.3) .attr("fill", "none") .append("title") .text("The original polygon") g.append("path") .attr("d", path(simplifiedCenterlineGeoJSON.geometry)) .attr("stroke", lineLegend.get("Longest path").stroke) .attr("stroke-width", 1.5) .attr("fill", "none") .append("title") .text(`The ${(strategy === "distat") ? "path between the most distant nodes" : "longest path"} on the graph`); const bounds = g.append("path") .attr("fill", "none") .attr("stroke", "dimgray") .attr("stroke-width", 1.75) .attr("d", voronoi.renderBounds()); svg.call(d3.zoom() .extent([[0, 0], [width, height]]) .scaleExtent([1, 8]) .on("zoom", zoomed)); function zoomed(event) { g.attr("transform", event.transform); } return svg.node(); }38:[{"id":"862c1ff75e2e910db7997a2199a3a5be9136682b36768877a8d44ef7d4037046ed1b885f22df43ddbfbab18a098ec15fb02ca0d7e997fec6cf996e3144b16082","href":"https://static.observableusercontent.com/files/862c1ff75e2e910db7997a2199a3a5be9136682b36768877a8d44ef7d4037046ed1b885f22df43ddbfbab18a098ec15fb02ca0d7e997fec6cf996e3144b16082","name":"amsterdamtransitmap@1.json","type":"application/json","lastModified":1593800820678.582,"size":269547},{"id":"811267b503e9174246754812702d4c5af5b04bd586967560b24db4effaa0e75212c4367f433f401297f570839c81be9eebc207345a01f8703ec4627a0bab8bf2","href":"https://static.observableusercontent.com/files/811267b503e9174246754812702d4c5af5b04bd586967560b24db4effaa0e75212c4367f433f401297f570839c81be9eebc207345a01f8703ec4627a0bab8bf2","name":"baddonut@1.geojson","type":"application/geo+json","lastModified":1599351553862.96,"size":1690}] 2e:{"title":"Centerlines from a GeoJSON Polygon with d3-delaunay","cells":[{"id":0,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":4190,"value":"// html`

`","pinned":true,"hidden":false,"mode":"ojs"},{"id":2268,"value":"md`## The Original geoJSON polygons\nPick one!\n\n0 to 8 are the polygons from Noah Veltman's [Centerline Labeling notebook](https://observablehq.com/@veltman/centerline-labeling).\n\n9 to 54 are polygons of Amsterdam's [GVB metro/tram rail transit map](https://www.gvb.nl/sites/default/files/railkaart2020.pdf) [converted to SVG and then to geoJSON](https://labs.mapbox.com/svg-to-geojson/).\n\nThe last one is a tricky horseshoe shape I drew myself. Herein nodes that are globally close are distant in the path, which trips up the usual strategy of finding the path between the most globally distant nodes. Number 31 is also a good examples of a shape in which the most-distant-node strategy fails.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1134,"value":"viewof num = radio({\n options: merged.features.map((f, i) => ({ label : f.properties.name ? f.properties.name : (i-9).toString(), value: i})),\n value: 13 //31\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7115,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":259,"value":"md`## Step 1: turn the shape into evenly spaced points\nFirst, we walk the perimeter of the shape to get a set of evenly spaced points that approximate the original shape. Adjust the slider to change how many points we're using (more points = more precision and more computation).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":273,"value":"viewof numPerimeterPoints = slider({\n min: 100,\n max: 1000,\n step: 40,\n value: 275,\n title: \"Number of perimeter points\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":278,"value":"svg``","pinned":false,"hidden":false,"mode":"ojs"},{"id":3640,"value":"md`First, apply projection to the polygon, so we can work with the SVG coordinates`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1301,"value":"outerRing = {\n const {\n geometry: { \n type,\n coordinates\n } \n } = feature \n \n if (type == \"MultiPolygon\") {\n return coordinates[0][0] // only take first feature for now\n .slice(1)\n .map(applyProjection);\n } else if (type == \"Polygon\") { // assume its a polygon\n return coordinates[0]\n .slice(1) // Ignore first point because it's in here twice to close the polygon\n .map(applyProjection); \n } else {\n throw new Error(\"Can only handle Polygon or MultiPolygon geoJSON types\");\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3648,"value":"md`Now get evenly spaced points along the perimeter of the polygon`","pinned":false,"hidden":false,"mode":"ojs"},{"id":127,"value":"pointsOnPolygon = getPointsAlongPolyline(outerRing, numPerimeterPoints)","pinned":true,"hidden":false,"mode":"ojs"},{"id":6123,"value":"lineLegend = new Map([\n [\"Voronoi line\", {stroke: \"#444\", fill: \"none\", \"stroke-width\": \"2px\" }],\n [\"Path most distant nodes\", {stroke: \"dodgerblue\", fill: \"none\", \"stroke-width\": \"2px\"}],\n [\"Original polygon\", {stroke: \"plum\", fill: \"none\", \"stroke-width\": \"2px\"}],\n [\"Longest path\", {stroke: \"gold\", fill: \"none\", \"stroke-width\": \"2px\"}]\n])","pinned":false,"hidden":false,"mode":"ojs"},{"id":3243,"value":"pointLegend = new Map([\n [\"Voronoi node\", {fill: \"dimgray\", r: 1}],\n [\"Inner node\", {fill: \"orange\", r: 1}],\n [\"Polygon perimeter node\", {fill: \"plum\", r: 1}],\n [\"Clipped node\", {fill: \"green\", r: 1}]\n]);","pinned":false,"hidden":false,"mode":"ojs"},{"id":1846,"value":"md`## Step 2. Voronoi Diagram\n\nWe use the points on the perimeter to make a Voronoi diagram with d3-delaunay\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3271,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":1671,"value":"// \nsvg``","pinned":false,"hidden":false,"mode":"ojs"},{"id":1657,"value":"voronoi = {\n const [xMin, xMax] = d3.extent(pointsOnPolygon.map(d => d[0]));\n const [yMin, yMax] = d3.extent(pointsOnPolygon.map(d => d[1]));\n \n // 1 px padding to prevent points on bounds from appearing inside polygon\n const padding = 5\n const bounds = [xMin - padding, yMin - padding, xMax + padding, yMax + padding];\n\n return d3.Delaunay\n .from(pointsOnPolygon)\n .voronoi(bounds);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2139,"value":"md`\nvoronoi.cellPolygon(i) returns the convex, closed polygon \\`[[x0, y0], [x1, y1], …, [x0, y0]]\\` representing the cell for the specified point i.\n\nvoronoi.cellPolygons() Returns an iterable over the non-empty polygons for each cell, with the cell index as property [[1](#footnote1)]. For now, let's turn this iterable into array\n\n[^](#paragraph1) *1:* [d3-delaunay docs](https://github.com/d3/d3-delaunay#voronoi_cellPolygons)\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2133,"value":"allCellPolygons = Array.from(voronoi.cellPolygons())","pinned":true,"hidden":false,"mode":"ojs"},{"id":3966,"value":"md`## Step 3. Clip Voronoi Diagram to polygon perimeter\n\n1. Loop over each Voronoi cell polygon\n2. Loop over each edge (Remember that edges are shared by adjacent polygons (half-edges))\n3. The order of the nodes in edges (i.e. the direction) matters! So check which nodes of edge are inside polygon\n3. If there is an intersection, replace the node with the intersection\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4595,"value":"addNodeIds = (edge) => {\n const [fromNode, toNode] = edge;\n fromNode.id = String(fromNode);\n toNode.id = String(toNode);\n \n return edge;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5392,"value":"isInsidePolygon = (node) => d3.polygonContains(pointsOnPolygon, node) ? true : false","pinned":true,"hidden":false,"mode":"ojs"},{"id":5389,"value":"// adds 'inside' boolean property to each node\ncheckNodesInside = (cellPolygon) => cellPolygon.map((node) => {\n node.inside = isInsidePolygon(node);\n return node\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":5384,"value":"function clipEdge(edge) {\n const [fromNode, toNode] = edge;\n\n // If only one of the two nodes is inside, we should clip \n const outward = (fromNode.inside !== toNode.inside && !fromNode.inside);\n const inward = (fromNode.inside !== toNode.inside && !toNode.inside);\n\n const clip = (outward) ? \"outward\" : (inward) ? \"inward\" : null;\n\n // find the intersection with feature polygon\n if (clip) {\n const { intersection } = findIntersectionWithFeature(fromNode, toNode);\n\n // why sometimes no intersection tho?\n if (intersection) intersection.clipped = true;\n\n // clip one edge of polygon (ie replace outside nodes) \n const newEdge = (clip === \"outward\") ? [intersection, toNode] : [fromNode, intersection];\n // newEdge.distance = distance\n return newEdge;\n }\n\n return edge;\n} ","pinned":true,"hidden":false,"mode":"ojs"},{"id":5184,"value":"function* edgeClipper(cellPolygonEdges) {\n for (const edge of cellPolygonEdges) {\n const [fromNode, toNode] = edge;\n\n // If only one of the two nodes is inside, we should clip \n const outward = (fromNode.inside !== toNode.inside && !fromNode.inside);\n const inward = (fromNode.inside !== toNode.inside && !toNode.inside);\n\n const clip = (outward) ? \"outward\" : (inward) ? \"inward\" : null;\n\n // find the intersection with feature polygon\n if (clip) {\n const { intersection } = findIntersectionWithFeature(fromNode, toNode);\n\n // why sometimes no intersection tho?\n if (intersection) intersection.clipped = true;\n\n // clip one edge of polygon (ie replace outside nodes) \n const newEdge = (clip === \"outward\") ? [intersection, toNode] : [fromNode, intersection];\n // newEdge.distance = distance\n return newEdge;\n }\n\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":4068,"value":"clipCellPolygon = (cellPolygon) => {\n // check which nodes of cellPolygon are inside feature polygon and appends boolean\n const cellPolygonWithInside = checkNodesInside(cellPolygon);\n \n // turn nodes into edges\n const cellPolygonEdges = d3.pairs(cellPolygonWithInside);\n \n const clippedCellPolygonEdges = cellPolygonEdges\n .map(clipEdge)\n .filter(([fromNode, toNode]) => { // discard outside edge \n if (fromNode && toNode) { // edge is not always defined\n const bothNodesOutside = fromNode.inside || toNode.inside;\n return bothNodesOutside;\n }\n return false;\n })\n \n // tack index back on, including index 0\n // cellPolygon.index ?? clippedCellPolygonEdges.index = cellPolygon.index // nullish coalescing operator not supported yet in Observable?\n if (typeof cellPolygon.index !== \"undefined\") clippedCellPolygonEdges.index = cellPolygon.index\n \n return clippedCellPolygonEdges;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2147,"value":"// //Array of CellPolygons (which is an array of points), not edges.\n// md`allCellPolygons consists of ${allCellPolygons.length} polygons, with polygons of ${d3.extent(allCellPolygons.map(polygon => polygon.length-1)).join(\" to \")} edges, each containing x,y points. Therefore, shared edges of polygons are in here multiple times!\n// `","pinned":false,"hidden":false,"mode":"ojs"},{"id":4290,"value":"legendSvg.cloneNode(true)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3310,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":5192,"value":"function* cellPolygonClipper(polygons) {\n for (const oneCellPolygon of polygons) {\n const clippedCellPolygon = clipCellPolygon(oneCellPolygon).map(addNodeIds)\n clippedCellPolygon.index = oneCellPolygon.index; // gets lost when mapped over\n // yield await Promises.delay(1000, clippedCellPolygon);\n yield clippedCellPolygon;\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":5194,"value":"clippedCellPolygons = Array.from(cellPolygonClipper(allCellPolygons))","pinned":true,"hidden":false,"mode":"ojs"},{"id":6140,"value":"md`#### Now filter out all duplicated edges`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5546,"value":"// dont care about cellPolygons only about unique edges\nedgesWithIdsAll = clippedCellPolygons\n .flat()\n .map(edge => edge\n .sort((a, b) => (a.id < b.id) ? 1 : -1 ) // sort nodes inside edges \n )\n .map(edge => {\n edge.id = edge[0].id + \"|\" + edge[1].id;\n return edge\n }) ","pinned":true,"hidden":false,"mode":"ojs"},{"id":5573,"value":"uniqueEdges = {\n // I love Map()!\n const edgeMap = new Map();\n edgesWithIdsAll.map(edge => edgeMap.set(edge.id, [edge[0], edge[1]] ) );\n \n return Array.from(edgeMap.values());\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":5649,"value":"uniqueNodes = {\n // Man, I really love Map()!\n const allNodes = uniqueEdges.flat()\n \n const nodesMap = new Map();\n allNodes.map(node => nodesMap.set(node.id, node )); \n \n return Array.from(nodesMap.values())\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":4256,"value":"clippedNodes = uniqueNodes.filter(e => e.clipped)","pinned":true,"hidden":false,"mode":"ojs"},{"id":4841,"value":"innerNodes = uniqueNodes.filter(e => e.inside)","pinned":true,"hidden":false,"mode":"ojs"},{"id":3429,"value":"legendSvg.cloneNode(true)","pinned":false,"hidden":false,"mode":"ojs"},{"id":5347,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":3826,"value":"md`## Step 4. Make a graph from the clipped and inner nodes\n\nConstruct a graph from the nodes using [ngraph](https://github.com/anvaka/ngraph.path)\n\n*Sometimes there are no Voronoi edges in bends?*\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4304,"value":"graph = {\n const graph = ngraph.ngraph()\n \n uniqueNodes.forEach(node => { graph.addNode(node.id, { x: node[0], y: node[1] }); })\n \n uniqueEdges.forEach(([fromNode, toNode]) => {\n const dx = fromNode[0] - toNode[0];\n const dy = fromNode[1] - toNode[1];\n \n const distance = Math.sqrt(dx * dx + dy * dy);\n \n graph.addLink(fromNode.id, toNode.id, { distance }) \n })\n \n return graph;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":5710,"value":"pathfinderAGreedy = ngraph.path.aGreedy(graph, {\n distance(fromNode, toNode) { \n // In this case we have coordinates. Lets use them as\n // distance between two nodes:\n let dx = fromNode.data.x - toNode.data.x;\n let dy = fromNode.data.y - toNode.data.y;\n\n return Math.sqrt(dx * dx + dy * dy);\n },\n heuristic(fromNode, toNode) {\n return Math.hypot(fromNode.data.x - toNode.data.x, fromNode.data.y - toNode.data.y);\n }\n });","pinned":true,"hidden":false,"mode":"ojs"},{"id":4386,"value":"pathfinderAStar = {\n let distancePathFinder = ngraph.path.aStar(graph, {\n distance(fromNode, toNode) {\n // In this case we have coordinates. Lets use them as\n // distance between two nodes:\n let dx = fromNode.data.x - toNode.data.x;\n let dy = fromNode.data.y - toNode.data.y;\n\n return Math.sqrt(dx * dx + dy * dy);\n },\n heuristic(fromNode, toNode) {\n // this is where we \"guess\" distance between two nodes.\n return Math.hypot(fromNode.data.x - toNode.data.x, fromNode.data.y - toNode.data.y);\n }\n });\n\n return distancePathFinder;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":6570,"value":"md`## Step 4a Prune perimeter nodes by using line simplification \n\nPathfinding between all nodes is expensive, so throw a bunch of nodes out. The idea comes from [Visvalingam’s algorithm line simplification algorithm](https://bost.ocks.org/mike/simplify/)\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":6704,"value":"viewof areaThreshold = slider({\n value: 0.015,\n title: \"Threshold triangle area for pruning graph nodes\",\n min: 0.001,\n max: 10\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":6902,"value":"md`Pruned ${clippedNodes.length - prunedNodes.length} out of ${clippedNodes.length} nodes`","pinned":false,"hidden":false,"mode":"ojs"},{"id":7664,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":6577,"value":"pointAreas = {\n const length = clippedNodes.length\n \n const areasArray = [];\n for (let leftIndex = 0, index = 1, rightIndex = 2; rightIndex < length; leftIndex++, index++, rightIndex++) {\n // the triangle around point b\n const area = triangleArea(clippedNodes[leftIndex], clippedNodes[index], clippedNodes[rightIndex]);\n \n areasArray.push({ area, leftIndex, index, rightIndex })\n }\n \n return areasArray;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":6870,"value":"noZeroArea = pointAreas.filter(d => d.area > areaThreshold)","pinned":true,"hidden":false,"mode":"ojs"},{"id":6878,"value":"prunedNodes = noZeroArea.map(({index}) => clippedNodes[index])","pinned":true,"hidden":false,"mode":"ojs"},{"id":6596,"value":"function triangleArea(a, b, c) {\n const [ax, ay] = a;\n const [bx, by] = b;\n const [cx, cy] = c;\n var area = Math.abs(((ay - cy) * (bx - cx) + (by - cy) * (cx - ax)) / 2);\n return area;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":4965,"value":"md`

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":6176,"value":"md`## 5. Find the path between the most distant pair of nodes`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4735,"value":"md`For each pair of perimeter nodes, find the shortest path between them. Keep the longest one we find. This strategy requires a lot more computation but it sure to give the longest path. The longest might not always be a 'centered line' though. For example for Florida or Uzbekistan the longest long will twirl up into their addendums.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":7309,"value":"viewof nodesForPathfinding = radio({\n options: [\n { label: 'All clipped nodes (expensive)', value: 'all' },\n { label: 'Pruned nodes (cheaper, but not exhaustive)', value: 'pruned' },\n ],\n value: \"pruned\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7061,"value":"viewof strategy = radio({\n options: [\n { label: 'Longest path (expensive but exhaustive)', value: 'longest' },\n { label: 'Path between most distant nodes (a good cheap guess)', value: 'distant' },\n ],\n value: \"distant\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":3205,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":7737,"value":"md`### Find the most (globally) distant nodes\n\nFor each pair of ${nodesForPathfinding} clipped nodes, calculate the global distance between them. We do this by looping over all clipped nodes and calculating the distance of that node to every other node. These results are stored in a nested map. Then find the path between the most distant pair of clipped nodes, this is often (but not always!) a good guess for longest path.\n\nIt would be even better to get node-to-node lengths inside of the polygon but I don't know how that would work.`\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":5831,"value":"globalNodeToNodeDistances = {\n const distanceToOtherClippedNodesPerNode = new Map();\n \n const nodes = (nodesForPathfinding === \"all\") ? clippedNodes : prunedNodes;\n\n for (let i = 0; i < nodes.length; i++) {\n const start = nodes[i];\n\n const distanceToOtherNodes = new Map();\n\n const nextNodeIndex = i + 1\n nodes.slice(nextNodeIndex).forEach((node, j) => {\n distanceToOtherNodes.set(nextNodeIndex + j, distanceBetween(start, node)) \n }) \n\n distanceToOtherClippedNodesPerNode.set(i, distanceToOtherNodes);\n\n }\n\n return distanceToOtherClippedNodesPerNode\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":6062,"value":"// make an array so we can sort node-to-node by distance\nsortedNodePairDistances = {\n const allNodeDistances = [];\n\n for (const [fromIndex, nodeMap] of globalNodeToNodeDistances) {\n for (const [toIndex, distance] of nodeMap) {\n allNodeDistances.push({ fromIndex, toIndex, distance })\n }\n } \n\n return allNodeDistances.sort(((a, b) => d3.descending(a.distance, b.distance)))\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":6156,"value":"mostGloballyDistantNode = sortedNodePairDistances[0]","pinned":true,"hidden":false,"mode":"ojs"},{"id":5887,"value":"mostDistantNodesGraphPath = {\n if (strategy === \"distant\") {\n const { fromIndex, toIndex } = mostGloballyDistantNode;\n const fromNodeHash = clippedNodes[fromIndex].id;\n const toNodeHash = clippedNodes[toIndex].id;\n \n console.log(fromNodeHash, toNodeHash)\n \n return pathfinderAGreedy.find(fromNodeHash, toNodeHash)\n } else {\n return []\n }\n} ","pinned":false,"hidden":false,"mode":"ojs"},{"id":6007,"value":"mostDistantNodesPath = mostDistantNodesGraphPath.map(node => ([node.data.x, node.data.y]))","pinned":true,"hidden":false,"mode":"ojs"},{"id":6180,"value":"md`### Find the longest path between pair of nodes\n\nCalculate the path between each pair of nodes and its length. Return the one that produces the longest path. \nThis requires a lot more compution but ensures the longest path is returned.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":6012,"value":"sumPathLength = line => line.map(e => distanceBetween(e[0], e[1])).reduce((acc, next) => acc + next, 0)","pinned":false,"hidden":false,"mode":"ojs"},{"id":6088,"value":"sumGraphPathLength = (graphPath) => \n sumPathLength(d3.pairs(graphPath.map(node => ([node.data.x, node.data.y]))))","pinned":false,"hidden":false,"mode":"ojs"},{"id":6077,"value":"longestGraphPath = {\n \n // heavy operation so only do it on toggle\n if (strategy === \"longest\") {\n let currentLongestPath = { path: null, length: 0 };\n\n sortedNodePairDistances.map((nodeToNode) => {\n const { fromIndex, toIndex, distance } = nodeToNode;\n\n const firstNodeHash = clippedNodes[fromIndex].id;\n const toNodeHash = clippedNodes[toIndex].id;\n\n const currentPath = pathfinderAGreedy.find(firstNodeHash, toNodeHash)\n const currentLength = sumGraphPathLength(currentPath);\n\n if (currentLength > currentLongestPath.length) {\n // console.log(currenLength)\n currentLongestPath = { path: currentPath, length: currentLength };\n }\n\n })\n\n return currentLongestPath\n } else \n return {path: []}\n} ","pinned":true,"hidden":false,"mode":"ojs"},{"id":6240,"value":"longestPath = longestGraphPath.path.map(node => ([node.data.x, node.data.y]));","pinned":true,"hidden":false,"mode":"ojs"},{"id":6226,"value":"md`## Step 6. Simplify the line & convert to geoJSON\nWe can use [simplify-js](http://mourner.github.io/simplify-js/) and a basis spline to smooth out our best candidate, to remove the jageddness.`\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":5949,"value":"viewof simplification = slider({\n min: 0.001, \n max: 10, \n value: 0.1,\n title: \"Simplification tolerance\",\n description: \"How much should we simplify the best line once we find it?\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2894,"value":"simplifiedCenterline = (strategy === \"longest\") ? simplify(longestPath) : simplify(mostDistantNodesPath)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2898,"value":"simplifiedCenterlineGeoJSON = turf.lineString(simplifiedCenterline.map(removeProjection))","pinned":true,"hidden":false,"mode":"ojs"},{"id":2960,"value":"md`#### Resulting Centerline geoJSON\n**Protip**: Download the geoJSON by clicking the three dots left of the cell and clicking \"Download JSON\"\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":6413,"value":"legendSvg.cloneNode(true)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2963,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":2,"value":"md`## Appendix`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1363,"value":"height = Math.min(width / 2.5, 500)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1360,"value":"railLinesGeoJSON = FileAttachment(\"amsterdamtransitmap@1.json\").json()","pinned":false,"hidden":false,"mode":"ojs"},{"id":6049,"value":"badDonut = FileAttachment(\"baddonut@1.geojson\").json()","pinned":false,"hidden":false,"mode":"ojs"},{"id":5025,"value":"places = d3.json(\"https://gist.githubusercontent.com/veltman/644f16a90259a20a88b036ef189d71fd/raw/d2f0a027bfc9b63ca223b509bd2cfe0cf5d138c2/places.geojson\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":5046,"value":"merged = ({type: \"FeatureCollection\", \"features\": [...places.features, ...railLinesGeoJSON.features, badDonut.features[0]]});","pinned":false,"hidden":false,"mode":"ojs"},{"id":1250,"value":"feature = merged.features[num] //railLinesGeoJSON","pinned":false,"hidden":false,"mode":"ojs"},{"id":6515,"value":"md`#### SVG Paths`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1357,"value":"polygonPath = path(feature)","pinned":false,"hidden":false,"mode":"ojs"},{"id":6519,"value":"delaunayPath = voronoi.delaunay.render()","pinned":false,"hidden":false,"mode":"ojs"},{"id":6502,"value":"voronoiPath = voronoi.render()","pinned":false,"hidden":false,"mode":"ojs"},{"id":5997,"value":"voronoiBounds = voronoi.renderBounds()","pinned":false,"hidden":false,"mode":"ojs"},{"id":2990,"value":"md`### Functions`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5937,"value":"line = d3.line().curve(d3.curveBasis)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3516,"value":"md`#### Intersection`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3853,"value":"findIntersectionWithFeature = findClosestPolygonIntersection(pointsOnPolygon)","pinned":true,"hidden":false,"mode":"ojs"},{"id":3522,"value":"findClosestPolygonIntersection = (polygon) => (startNode, endNode) => {\n \n return polygon.reduce((best, point, i) => {\n const intersection = findIntersection(startNode, endNode, point, polygon[i + 1] || polygon[0]);\n \n if (intersection) {\n const distance = distanceBetween(startNode, intersection);\n \n if (!best.distance || distance < best.distance) {\n return { intersection, distance };\n }\n \n }\n return best;\n }, {});\n \n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3519,"value":"function findIntersection(a1, a2, b1, b2) {\n // Adapted from https://github.com/Turfjs/turf-line-slice-at-intersection\n const uaT = (b2[0] - b1[0]) * (a1[1] - b1[1]) - (b2[1] - b1[1]) * (a1[0] - b1[0]);\n const ubT = (a2[0] - a1[0]) * (a1[1] - b1[1]) - (a2[1] - a1[1]) * (a1[0] - b1[0]);\n const uB = (b2[1] - b1[1]) * (a2[0] - a1[0]) - (b2[0] - b1[0]) * (a2[1] - a1[1]);\n\n if (uB !== 0) {\n const ua = uaT / uB;\n const ub = ubT / uB;\n \n if (ua > 0 && ua < 1 && ub > 0 && ub < 1) {\n return [a1[0] + ua * (a2[0] - a1[0]), a1[1] + ua * (a2[1] - a1[1])];\n }\n \n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3513,"value":"md`#### Projection`","pinned":false,"hidden":false,"mode":"ojs"},{"id":106,"value":"// Rescale to fit at the current screen width\nprojection = d3\n .geoIdentity()\n .fitExtent([[5, 5], [width - 5, height - 5]], feature)\n","pinned":true,"hidden":false,"mode":"ojs"},{"id":5966,"value":"path = d3.geoPath().projection(projection)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2753,"value":"s = projection.scale();","pinned":true,"hidden":false,"mode":"ojs"},{"id":2755,"value":"t = projection.translate();","pinned":true,"hidden":false,"mode":"ojs"},{"id":2969,"value":"applyProjection = ([x, y]) => [s * x + t[0], s * y + t[1]]","pinned":true,"hidden":false,"mode":"ojs"},{"id":2967,"value":"removeProjection = ([x, y]) => [(x - t[0])/s, (y - t[1])/s]","pinned":true,"hidden":false,"mode":"ojs"},{"id":130,"value":"// From https://observablehq.com/@veltman/centerline-labeling\nfunction getPointsAlongPolyline(polyline, count) {\n // Calculate distance between each point of polyline\n const distances = polyline.map((p, i) => {\n return distanceBetween(p, polyline[i + 1] || polyline[0]);\n });\n\n const totalLength = d3.sum(distances);\n const stepsize = totalLength / count;\n let traversed = 0;\n let next = stepsize / 2;\n\n const done = polyline.reduce((arr, currentPoint, i) => {\n const distanceBetweenCurrentPoints = distances[i];\n while (next < traversed + distanceBetweenCurrentPoints) {\n let nextPoint = polyline[i + 1] || polyline[0];\n\n\n let percent = (next - traversed) / distanceBetweenCurrentPoints;\n \n arr.push([\n currentPoint[0] + (nextPoint[0] - currentPoint[0]) * percent,\n currentPoint[1] + (nextPoint[1] - currentPoint[1]) * percent\n ]);\n \n next += stepsize;\n }\n traversed += distanceBetweenCurrentPoints;\n return arr;\n }, []);\n return done;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":981,"value":"function distanceBetween(a, b) {\n const dx = a[0] - b[0],\n dy = a[1] - b[1];\n return Math.sqrt(dx * dx + dy * dy);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2832,"value":"md`## Imports`","pinned":false,"hidden":false,"mode":"ojs"},{"id":7874,"value":"require('https://bundle.run/ngraph.path@1.4.0')","pinned":true,"hidden":false,"mode":"ojs"},{"id":4274,"value":"ngraph = ({\n ngraph: await require('https://bundle.run/ngraph.graph@20.0.0'),\n path: await require('https://bundle.run/ngraph.path@1.4.0'),\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":2596,"value":"// TODO only import turf modules used in this notebook\nturf = require(\"@turf/turf@5\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":5947,"value":"function simplify(points) {\n // Convert from [x, y] to { x, y } and back for simplify-js\n return simplifyJS(points.map(p => ({ x: p[0], y: p[1] })), simplification).map(p => [p.x, p.y]); \n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5944,"value":"simplifyJS = require(\"simplify-js\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":7,"value":"d3 = require(\"d3@v6\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":572,"value":"import {slider, radio} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1152,"value":"styles = html``","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W38","stdlib":"1","titling":"auto"}