30:Tf36,dggs = md`[H3](https://eng.uber.com/h3/) is an open-source [discrete global grid system](https://en.wikipedia.org/wiki/Discrete_global_grid) released by Uber in 2018, with applications such as the quick and precise indexation of points on the sphere, useful for all kinds of analysis of geographic data, and for beautiful “hexgrid” maps. Creating a DGGS is definitely [not an easy task](http://webpages.sou.edu/~sahrk/sqspc/pubs/gdggs03.pdf): ideally we would want a regular grid, something made of triangles, squares or hexagons ⬡ that would all have the same shape. This works great on the plane — but on the sphere it is a mathematical impossibility: any such system will have to include special points to “take in” the sphere’s curvature. In the case of an hexagonal mesh, [as Leonhard Euler demonstrated](https://quantixed.org/2018/06/27/pentagrammarspin-why-twelve-pentagons/), the grid must include exactly 12 pentagons ⬠. And darn! these pentagons have to be put somewhere… even worse, the whole grid has to be built around them. H3 does not escape this problems, and follows the usual “soccer ball” construction. It starts with a regular icosahedron: a shape of 20 equilateral triangles joining the 12 Euler points. The shape is rotated, as per Buckminster Fuller’s original design, so that the Euler points all fall in the water at the corners of the [Airocean projection](https://observablehq.com/@fil/airocean-projection). Then, on each of these triangular faces, H3 defines, for each “resolution” (a number between 1 and 16) a certain number of polygons. Seen as a graph, each of these polygons has 6 neighbours (in that sense, they are topologically hexagons) — except those centered on the 12 Euler points where the triangular faces meet, which have five neighbors (pentagons). In H3, the spherical polygons are not directly defined on the sphere, but are drawn on a [gnomonic projection](https://uber.github.io/h3/#/documentation/core-library/coordinate-systems) of each face (centered on the face’s centroid), as a regular grid of hexagons, which is then projected back to the sphere. The gnomonic projection’s role is to map all geodesic (aka “great circle”) lines on the sphere to straight lines on the plane, and vice versa. As a consequence, the spherical polygons created this way have a neat spherical definition. However, this technique has two problematic implications: First, the gnomonic projection creates important distortions between the centers and the edges, leading to [quite a lot of areal variation](https://github.com/uber/h3/issues/168) between the polygons: in this case, about 2 decibels (dB) — in other words the largest hexagons (in red, near the faces’ centers) are about 1.6 times bigger than the smallest (in green, near the control points). This means that geostatistics (done by aggregating values on these cells) will not consider the same areas if you work on, say, Perth (near the Euler point) and Sydney (near the face’s center) — two cities in the same country [or, for other faces: Boston and Tucson, Tripoli and Lagos…]. **Extra vertices.** Second, at odd resolution values, extra “[distortion vertices](https://github.com/uber/h3/issues/209)” must be added the spherical polygons that reside on two or three planar/gnomonic faces, to make sure that they don’t look “bent”. So, technically, those hexagons will be heptagons and octogons, and pentagons will be spherical decagons. This does not happen at even resolutions, because the hexagons positions are symmetrical with respect to the triangles’ edges. **Gray-Fuller grid.** Compare this with the [Gray-Fuller spatial grid](https://observablehq.com/@fil/gray-fuller-grid): [visually](https://observablehq.com/@fil/gray-fuller-grid-oddities) and with [metrics](https://observablehq.com/@fil/gray-fuller-grid-metrics). `31:T75d,function chart() { const context = DOM.context2d(width, height), path = d3.geoPath(projection).context(context); context.beginPath(); path({ type: "Sphere" }); context.fillStyle = color(0); context.fill(); context.lineWidth = 3; context.stroke(); context.clip(); context.lineWidth = 0.5; for (const c of contours) { context.beginPath(); path(c); context.strokeStyle = context.fillStyle = color(c.value); context.fill(); context.stroke(); } context.beginPath(); context.fillStyle = "lightblue"; for (const penta of hexagons.features.filter(d => d.properties.pentagon)) path(penta); context.fill(); context.beginPath(); for (const hepta of hexagons.features.filter( d => d.geometry.coordinates[0].length == 8 )) { path({ type: "Point", coordinates: d3.geoCentroid(hepta) }); } context.fillStyle = "red"; context.fill(); context.beginPath(); context.fillStyle = "black"; for (const octo of hexagons.features.filter( d => d.geometry.coordinates[0].length == 9 )) path({ type: "Point", coordinates: d3.geoCentroid(octo) }); context.fill(); context.beginPath(); context.fillStyle = "lightblue"; for (const nona of hexagons.features.filter( d => d.geometry.coordinates[0].length == 10 )) path({ type: "Point", coordinates: d3.geoCentroid(nona) }); context.fill(); context.beginPath(); context.fillStyle = "darkorange"; for (const deca of hexagons.features.filter( d => d.geometry.coordinates[0].length == 11 )) path({ type: "Point", coordinates: d3.geoCentroid(deca) }); context.fill(); context.strokeStyle = "#fff"; context.beginPath(); path(hexagons); context.lineWidth = 0.5; context.stroke(); context.strokeStyle = "black"; context.beginPath(); path(land); context.lineWidth = 1.5; context.stroke(); return context.canvas; }32:T5b4,projection = ({ Airocean: d3 .geoAirocean() .angle(-150) .fitExtent([[2, 2], [width - 2, height - 2]], { type: "Sphere" }), "Gnomonic icosahedral": d3 .geoIcosahedral() .faceProjection(face => { const c = face.site.map(d => -d); c[2] = Math.abs(c[1] - 52.62) < 1 || Math.abs(c[1] + 10.81) < 1 ? 60 : 0; return d3 .geoProjection(d3.geoGnomonicRaw) .rotate(c) .translate([0, 0]); }) .rotate([-83.65929, 25.44458, -87.45184]) .angle(60) .fitExtent([[2, 2], [width - 2, height - 2]], { type: "Sphere" }), "Gray-Fuller icosahedral": d3 .geoIcosahedral() .faceProjection(face => { var c = face.site.map(d => -d); c[2] = Math.abs(c[1] - 52.62) < 1 || Math.abs(c[1] + 10.81) < 1 ? 60 : 0; return d3 .geoProjection(d3.geoGrayFullerRaw()) .rotate(c) .translate([0, 0]); }) .rotate([-83.65929, 25.44458, -87.45184]) .angle(60) .fitExtent([[2, 2], [width - 2, height - 2]], { type: "Sphere" }), Orthographic: d3 .geoOrthographic() .rotate([-110, -35]) .fitExtent([[2, 2], [width - 2, height - 2]], { type: "Sphere" }), "Azimuthal equal-area": d3 .geoAzimuthalEqualArea() .rotate([-110, -35]) .clipAngle(54) .fitExtent([[2, 2], [width - 2, height - 2]], { type: "Sphere" }), Mercator: d3 .geoMercator() .fitExtent([[2, 2], [width - 2, height - 2]], { type: "Sphere" }) }[projname])33:[] 2e:{"title":"H3 oddities","cells":[{"id":400,"value":"md`# H3 oddities\n\n**Summary:** just like there cannot be a perfect map projection, creating a regular discrete global grid on a sphere is mathematically impossible. As a consequence, any such system has to make choices with various trade-offs. In this notebook we examine Uber’s H3 grid system, which paves the globe with lots of hexagons ⬡ and the unavoidable 12 pentagons ⬠.\n\nThe map below shows the H3 grid at resolution 1, with area distortion indicated by the color scale for the common hexagons ⬡ ; the pentagons ⬠ are painted in light blue. Special polygons are marked with a dot: red • hexagons with 7 sides 🌶, black • hexagons with 8 sides 🎩, and orange • pentagons with… 10 sides 🤯.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":441,"value":"viewof projname = select({\n options: [\n \"Airocean\",\n \"Gnomonic icosahedral\",\n \"Gray-Fuller icosahedral\",\n \"Orthographic\",\n \"Azimuthal equal-area\",\n \"Mercator\"\n ],\n description: \"Projection for the map (doesn’t affect the analysis)\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":316,"value":"map = chart()","pinned":false,"hidden":false,"mode":"ojs"},{"id":407,"value":"legend({ color, title: \"Areal distortion (dB)\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":380,"value":"viewof resolution = select({\n options: [0, 1, 2, 3],\n value: 1,\n description: \"H3 resolution\"\n})\n// you can up it to 4 if you're very patient\n// at resolution 5 your browser will likely hang for a few dozen minutes…","pinned":false,"hidden":false,"mode":"ojs"},{"id":0,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":594,"value":"md`---\n\n_A note about the map’s projection:_\n\nOnly the visual rendering is affected when you modify the map’s projection. I chose **Airocean** because it's “almost” equal-area and fits well, by H3’s design, with the position of the 12 pentagons. However it can be hard to read because it itself has oddities. So I added an **azimuthal equal-area** projection to facilitate the reading of the hexagons’s varying areas, and, as [suggested by Max Friedrich Hartmann](https://observablehq.com/@mxfh/h3-oddities) (thanks!), a **gnomonic icosahedral** projection that shows from which geometry H3’s grid is coming. In that geometry, we can see that the hexagonal grid is perfectly regular. Max Friedrich also suggested a **Gray-Fuller icosahedral** projection with the same arrangement of faces, for easier comparison (the formula for the faces being the same as that of the Airocean projection). The **orthographic** and **Mercator** projections are there because they’re the most common.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":429,"value":"---\n\nThanks to [DJ Collier](https://observablehq.com/@dcollier) for updating this notebook to the new H3 API.\n\n_code_","pinned":false,"hidden":false,"mode":"md"},{"id":2,"value":"h3 = import(\"h3-js@4\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":6,"value":"hexworld = level => ({\n type: \"FeatureCollection\",\n features: h3.getRes0Cells()\n .map(i => h3.cellToChildren(i, level))\n .flat()\n .map(d => ({\n type: \"Feature\",\n properties: { id: d, pentagon: h3.isPentagon(d) },\n geometry: {\n type: \"Polygon\",\n coordinates: [h3.cellToBoundary(d, true).reverse()]\n }\n }))\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":194,"value":"hexagons = hexworld(resolution)","pinned":true,"hidden":false,"mode":"ojs"},{"id":9,"value":"contours = d3\n .geoContour()\n .x(d => d3.geoCentroid(d)[0])\n .y(d => d3.geoCentroid(d)[1])\n .value((d, i) => 10 * Math.log10(areas[i] / median))\n .thresholds(100)(hexagons.features)","pinned":true,"hidden":false,"mode":"ojs"},{"id":341,"value":"areas = hexagons.features.map(d => d3.geoArea(d))","pinned":true,"hidden":false,"mode":"ojs"},{"id":338,"value":"median = d3.median(areas)","pinned":true,"hidden":false,"mode":"ojs"},{"id":251,"value":"color = d3.scaleSequential(t => d3.interpolateRdYlGn(1 - t)).domain([-1, 1])","pinned":true,"hidden":false,"mode":"ojs"},{"id":22,"value":"$31","pinned":true,"hidden":false,"mode":"ojs"},{"id":12,"value":"d3 = require(\"d3@7\", \"d3-geo-voronoi@2\", \"d3-geo-polygon@1.8\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":280,"value":"land = {\n const world = await fetch(\"https://cdn.jsdelivr.net/npm/visionscarto-world-atlas@0.1.0/world/110m.json\").then(d => d.json());\n return topojson.feature(world, world.objects.countries);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":449,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":26,"value":"height = (width * (projname === \"Airocean\" ? 2.1 : 1)) | 0","pinned":true,"hidden":false,"mode":"ojs"},{"id":405,"value":"import { legend } from \"@d3/color-legend\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":443,"value":"import { select } from \"@jashkenas/inputs\"","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}