2d:null 2e:["$","title",null,{"children":"degenerate ring analysis | jleedev | Observable"}] 33:T497,render_globe = (countries, land) => { const sphere = { type: "Sphere" }; let projection = d3.geoNaturalEarth1(); projection.fitSize([Math.min(width, 400), 300], sphere); let bounds = d3.geoPath(projection).bounds(sphere); let cx = DOM.context2d( bounds[1][0] - bounds[0][0], bounds[1][1] - bounds[0][1] ); cx.translate(-bounds[0][0], -bounds[0][1]); let path = d3.geoPath(projection, cx); cx.fillStyle = "cornsilk"; cx.fillRect(0, 0, 400, 300); // Ocean cx.fillStyle = "#9cf"; cx.beginPath(), path({ type: "Sphere" }), cx.fill(); cx.save(), cx.clip(); // Land cx.strokeStyle = "#66c"; cx.fillStyle = "darkolivegreen"; cx.lineWidth = 3; cx.lineJoin = "round"; cx.beginPath(), path(land), cx.stroke(), cx.fill(); // Countries cx.strokeStyle = "#fff"; cx.lineWidth = 1; cx.fillStyle = "#000"; cx.beginPath(), path(countries), cx.fill(), cx.stroke(); cx.restore(); // Outline cx.strokeStyle = "#000"; cx.lineWidth = 0.5; cx.beginPath(), path({ type: "Sphere" }), cx.stroke(); cx.strokeStyle = "#0004"; cx.lineWidth = 0.5; cx.beginPath(), path(d3.geoGraticule10()), cx.stroke(); return cx.canvas; }34:T740,showRing = (ring) => { const size = 128; let margin = 2; let cx = DOM.context2d(size + margin * 2, size + margin * 2); let proj = d3.geoMercator(); let center = turf.center(ring); let area = d3.geoArea(ring); let area_m2 = (area > 2 * Math.PI ? 4 * Math.PI - area : area) * 6371008.8 ** 2; let length_m = d3.geoLength(ring) * 6371008.8; let turfArea = turf.area(ring); // Buffer a radius around the point to avoid d3 attempting a zero bounding box. // Explode to points to avoid turf being unable to rewind properly. // Fit to points and buf to avoid turf being unable to buffer properly. let buf = turf.rewind(turf.buffer(ring, 0.05), { reverse: true }); proj.fitSize( [size, size], turf.featureCollection([ ...turf.explode(buf).features, ...turf.explode(ring).features ]) ); let path = d3.geoPath(proj, cx); cx.translate(margin, margin); cx.fillStyle = "cornsilk"; cx.fillRect(0, 0, size, size); cx.strokeStyle = "cornflowerblue"; cx.beginPath(), path(buf), cx.stroke(); cx.strokeStyle = "#000"; cx.fillStyle = "#fff"; cx.beginPath(), path(ring), cx.stroke(); cx.beginPath(), path(turf.explode(ring)), cx.stroke(), cx.fill(); let kinks = turf.kinks(ring); cx.fillStyle = "red"; cx.beginPath(), path(kinks), cx.fill(), cx.stroke(); let captions = [ html`${ring.geometry.coordinates[0].length} points`, html`${+(area / Math.PI).toExponential(3)}π sr enclosed`, html`${+turfArea.toExponential(3)} m²`, html`length ${+length_m.toExponential(2)}m`, html`${kinks.features.length} kinks`, html`${dms(turf.getCoord(center))}` ]; return html`

${cx.canvas}

${ring.id ?? ""} ${ring.properties.name ?? ""}

${cap}

`; }35:T4d5,viewArc = (line) => { let { id, coords } = line; let path = d3.geoPath(d3.geoIdentity()).pointRadius(0.2); let bounds = path.bounds(line).flat(); let [x, y] = bounds.slice(0, 2); let w = bounds[2] - bounds[0]; let h = bounds[3] - bounds[1]; let scale = Math.max(w, h); path.pointRadius(scale / 10); let pad = 0.2 * scale; x -= pad; y -= pad; w += 2 * pad; h += 2 * pad; let viewBox = [x, y, w, h].join(); let points = turf.explode(line); let numbers = d3 .create("svg:g") .attr("font-size", 0.2 * scale) .attr("text-anchor", "middle") .attr("transform", `translate(0, ${scale * 0.07})`); numbers .selectAll("text") .data(points.features) .join("text") .attr("x", (p) => turf.getCoord(p)[0]) .attr("y", (p) => turf.getCoord(p)[1]) .text((_, i) => i); return html`

${id}

`; }36:[{"id":"eee5236563c7d3c84f8b5b6c28d7437ea62f4240f68439ae2c54b8cb8400fd5c76908f7fe90630cb1c0ed7c4060697e8237228c6165053af570ec9be417175f4","href":"https://static.observableusercontent.com/files/eee5236563c7d3c84f8b5b6c28d7437ea62f4240f68439ae2c54b8cb8400fd5c76908f7fe90630cb1c0ed7c4060697e8237228c6165053af570ec9be417175f4","name":"image@1.png","type":"image/png","lastModified":1661538795430.68,"size":68645},{"id":"a68d741a7430f2fb6ae3017cad0456e7dff06db99bb28e09c64af0aecc0d77fb2775aacdde872c14ff4aafec79783e49d6756ffbb6584656fb27493198d7c48f","href":"https://static.observableusercontent.com/files/a68d741a7430f2fb6ae3017cad0456e7dff06db99bb28e09c64af0aecc0d77fb2775aacdde872c14ff4aafec79783e49d6756ffbb6584656fb27493198d7c48f","name":"image.png","type":"image/png","lastModified":1661538792881.259,"size":31147}] 31:{"id":"b4e1b6c97d82b5cb","slug":"degenerate-ring-analysis","title":"degenerate ring analysis","thumbnail":"67dcbd54348957d4da4666bfcec4ed563a35e897879248bbc144da851106c472","access_level":"public","create_time":"$D2022-08-26T11:32:59.368Z","update_time":"$D2022-08-26T19:19:46.808Z","view_time":null,"first_public_time":"$D2022-08-26T15:25:36.693Z","has_unpublished_edits":null,"fork":false,"stars":0,"starred":false,"owner":{"type":"team","id":"93f88a392b09c0b2","login":"jleedev","name":"","avatar_url":"https://avatars-next.observableusercontent.com/68aa27cdda263ac7113092ff8aacbd9058bedc56790cabb734e9c32a2246dff1","tier":"starter_2024"},"version":1088,"publish_version":1088,"role":null,"authored":false,"authors":[{"type":"individual","id":"b5e0663ad99246fa","login":"jleedev","name":"","avatar_url":"https://avatars-next.observableusercontent.com/68aa27cdda263ac7113092ff8aacbd9058bedc56790cabb734e9c32a2246dff1","tier":null}],"members":[{"type":"individual","id":"b5e0663ad99246fa","login":"jleedev","name":"","avatar_url":"https://avatars-next.observableusercontent.com/68aa27cdda263ac7113092ff8aacbd9058bedc56790cabb734e9c32a2246dff1","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":0,"collection":null,"type":"notebook","body":{"title":"degenerate ring analysis","cells":[{"id":0,"value":"# degenerate ring analysis","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":603,"value":"The issue here is that `countries-10m` from [world-atlas] contains some inside-out rings, so d3 dutifully draws these polygons as encompassing the entire earth.\n\n[world-atlas]: https://www.npmjs.com/package/world-atlas","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":965,"value":"

![image.png](${await FileAttachment(\"image.png\").url()})\n

Broken data

\n\n

![image@1.png](${await FileAttachment(\"image@1.png\").url()}) \n

Fixed data

","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":565,"value":"// I’m kind of expensive to render, so that’s a screenshot!\n// md`${render_globe(countries, land)} ${render_globe(countries_fixed, land_fixed)}`","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":561,"value":"$33","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":599,"value":"This is visible in two ways:\n\n- `fill()` fills the entire planet.\n- `stroke()` draws the outline of the feature _and_ an outline around the entire planet.\n\nd3 is very sophisticated and can correctly handle geometry which crosses the poles or the antimeridian, while other tools break down — Leaflet and Maplibre, for instance, will draw antimeridian crossings as entire lines of latitude, and a polygon that contains the north pole will be inverted in a pretty hard-to-understand way.\n\nThe GeoJSON spec requires correct winding:\n\n> ```\n> o A linear ring MUST follow the right-hand rule with respect to the\n> area it bounds, i.e., exterior rings are counterclockwise, and\n> holes are clockwise.\n> ```\n\nBut tools which treat the earth as a rectangle are free to ignore this. d3 requires correct winding.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":79,"value":"areaTable = (countries) =>\n Inputs.table(\n countries.features.map((feat, i) => ({\n i,\n name: feat.properties.name,\n area: d3.geoArea(feat) * 6371008.8 ** 2\n })),\n { sort: \"area\", reverse: true, height: 500 }\n )","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":5,"value":"areaTable(countries)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":439,"value":"reversePolygon = ({ geometry: { coordinates, ...geometry }, ...feature }) => ({\n ...feature,\n geometry: {\n ...geometry,\n coordinates: coordinates.map(reversed)\n }\n})","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":23,"value":"reversed = (coordinates) => [...coordinates].reverse()","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":20,"value":"flattenRings = ({\n geometry: {\n coordinates: [outer, ...inners]\n },\n ...polygonFeature\n}) =>\n turf.featureCollection(\n [outer, ...inners.map(reversed)].map((ring) => ({\n ...polygonFeature,\n ...turf.polygon([ring]),\n properties: { ...polygonFeature.properties }\n }))\n )","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":83,"value":"multiFlatten = (collection) =>\n turf.featureCollection(\n turf\n .flatten(collection)\n .features.map(flattenRings)\n .flatMap((c) => c.features)\n )","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":14,"value":"degenerateRings = turf.featureCollection(\n multiFlatten(countries)\n .features.filter(\n (feat) => d3.geoArea(feat) > 2 * Math.PI || d3.geoArea(feat) == 0\n )\n .sort((a, b) => d3.geoArea(a) - d3.geoArea(b))\n)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":114,"value":"areaTable(degenerateRings)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":119,"value":"$34","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":339,"value":"Blue outline is a 50m buffer.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":122,"value":"md`## rings as given\n\n

${degenerateRings.features.map(showRing)}

`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":381,"value":"maldives_inverted = degenerateRings.features.slice(17)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":429,"value":"Inputs.table(\n degenerateRings.features.map((feature, i) => ({\n i,\n name: feature.properties.name,\n area: d3.geoArea(feature),\n reverseArea: d3.geoArea(reversePolygon(feature))\n })),\n { height: 500 }\n)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":452,"value":"It looks like d3 will decide that some of these polygons are sufficiently flat that they have no inside or outside.\n\nBut one of the Maldives polygons is very much not flat, and topojson wound it differently from d3; two of them sure look pretty flat to me, and the libraries disagree on which side of the triangle that third vertex falls on.\n\nAlso, Italy and Vatican are the same ring, they should be absolutely identical. Why is one 188m and one 376m? How does `d3.geoLength` depend on winding order?\n\nTopoJSON doesn’t actually rewind anything, it encodes the coordinates in the same order as in the input. But, it applies a quantization step, and that’s what produces inside-out rings.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":474,"value":"md`## all rings reversed\n\n

${degenerateRings.features\n .map(reversePolygon)\n .map(showRing)}

`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":482,"value":"Spain’s buffer is now inside out, but others are not! And the output of `d3.geoLength` has changed.\n\nApparently `turf.buffer` is actually capable of thinking that a polygon’s interior is its exterior. Probably not on purpose.\n\nLet’s back up to the encoded, quantized TopoJSON data:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":655,"value":"getOneArc = (arc) =>\n topojson.feature(world, { type: \"MultiPolygon\", arcs: [arc] })","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":660,"value":"badRings = world.objects.countries.geometries\n .flatMap((geom) => {\n if (geom.type == \"MultiPolygon\") {\n return geom.arcs.map((arcs) => ({ ...geom, type: \"Polygon\", arcs }));\n } else {\n return [geom];\n }\n })\n .flatMap((geom) => {\n let geo = getOneArc(geom.arcs);\n if (d3.geoArea(geo) > 2 * Math.PI) {\n return [geom];\n } else {\n return [];\n }\n })","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":669,"value":"arcIndices = badRings.map(({ arcs }) => arcs).flat(2)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":928,"value":"lookupArc = (id) => ({ id, arc: world.arcs[id < 0 ? ~id : id] })","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":679,"value":"arcs = arcIndices.map(lookupArc)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1052,"value":"getArc = ({ id, arc }) => {\n let coords = arc.map(topojson.transform(world.transform));\n return { id, coords, ...turf.lineString(coords) };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":686,"value":"transformArcRaw = ({ id, arc }) => {\n let coords = [[0, 0], ...arc.slice(1)].map(\n topojson.transform({ scale: [1, 1], translate: [0, 0] })\n );\n return {\n id,\n coords,\n ...turf.lineString(coords)\n };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":716,"value":"$35","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":756,"value":"html`

${arcs.map(transformArcRaw).map(viewArc)}

`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1033,"value":"kinky_arcs = world.objects.countries.geometries\n .flatMap((x) => (x.properties.name == \"Maldives\" ? [x] : []))\n .flatMap((obj) => (obj.type == \"MultiPolygon\" ? obj.arcs.flat() : obj.arcs))\n .map(lookupArc)\n .map(getArc)\n .flatMap((line) => (turf.kinks(line).features.length ? [line] : []))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":922,"value":"md`## kinky arcs\n\n

${kinky_arcs.map(viewArc)}

`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1039,"value":"md`## kinky rings\n\n

${kinky_arcs\n .map((line) => ({ ...line, ...turf.lineStringToPolygon(line) }))\n .map(showRing)}

`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":810,"value":"So at this point it’s obvious that TopoJSON is encoding integer Cartesian coordinates, and three points that fall on a straight line form a triangle whose area and two angles are precisely zero. The error is introduced when TopoJSON transforms the encoded coordinates to decimal latitude and longitude.\n\nExterior rings must be counterclockwise — visually, TopoJSON sure seems wrong about the first one.\n\nTurf produces some really weird buffer results when polygons self-intersect.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":130,"value":"","pinned":true,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":173,"value":"formatcoords = (await import(\"https://cdn.skypack.dev/formatcoords\")).default","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":184,"value":"dms = ([x, y]) => formatcoords(y, x).format(\"DDMMssX\", { decimalPlaces: 0 })","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":204,"value":"html`${dms([-80, 40.44])}`","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":494,"value":"world = fetch(\"https://cdn.jsdelivr.net/npm/world-atlas@2/countries-10m.json\").then((r) => r.json())","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":3,"value":"countries = topojson.feature(world, world.objects.countries)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":539,"value":"land = topojson.feature(world, world.objects.land)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":551,"value":"import { countries_fixed, land_fixed } from \"15347f40054e706a\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":12,"value":"turf = require(\"@turf/turf\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W36","stdlib":"1","titling":"auto"}} 32:"$31:body"