Mathematical Foundations

Definition

The Lambert azimuthal equal-area projection is defined by a mathematical constraint, not a geometric construction. For a projection centered at point \(S\):

• Azimuthal property: The direction (azimuth) from \(S\) to any point \(P\) on the sphere equals the direction from the origin to \(P'\) on the plane

• Equal-area property: Any region on the sphere maps to a region of identical area on the plane

These two constraints uniquely determine the radial distance formula. For a point \(P\) at angular distance \(c\) from the center (measured along the sphere surface), the projected distance from the origin is:

\[\rho = 2R \sin\left(\frac{c}{2}\right)\]

This is not a perspective projection and has no simple geometric interpretation like “chord distance” or “ray intersection.” The formula is derived analytically from the equal-area constraint (Snyder, 1987, p. 182-185).

Forward Formulas

For the oblique aspect centered at \((\lambda_0, \phi_1)\), the full formulas are (Snyder, 1987, eq. 24-2 to 24-4, p. 185):

\[k' = \sqrt{\frac{2}{1 + \sin\phi_1\sin\phi + \cos\phi_1\cos\phi\cos(\lambda - \lambda_0)}}\]

\[x = R \cdot k' \cdot \cos\phi \cdot \sin(\lambda - \lambda_0)\]

\[y = R \cdot k' \cdot [\cos\phi_1\sin\phi - \sin\phi_1\cos\phi\cos(\lambda - \lambda_0)]\]

Equal-Area Proof

The projection preserves area because the radial and tangential scale factors satisfy \(h' \cdot k' = 1\) at every point. At angular distance \(c\) from center (Snyder, 1987, eq. 24-22, 24-23, p. 188):

\[h' = \cos\left(\frac{c}{2}\right), \quad k' = \sec\left(\frac{c}{2}\right)\]

Therefore \(h' \cdot k' = \cos(c/2) \cdot \sec(c/2) = 1\), confirming the equal-area property.

Each colored band has equal area on the sphere (the band boundaries are at \(\cos\phi_i = 1 - i/5\)). After Lambert projection the outer bands stretch radially and compress tangentially, and every annulus still has the same area: \(\rho_i^2 = 2R^2(1 - \cos\phi_i)\) grows linearly with \(i\).

Inverse Formulas

The inverse mapping \((x, y) \to (\lambda, \phi)\) recovers geographic coordinates from planar coordinates (Snyder, 1987, eq. 24-14 to 24-16, p. 187):

\[\rho = \sqrt{x^2 + y^2}\]

\[c = 2\arcsin\left(\frac{\rho}{2R}\right)\]

where \(c\) is the angular distance from the projection center. Then:

\[\phi = \arcsin\left[\cos c \cdot \sin\phi_1 + \frac{y \cdot \sin c \cdot \cos\phi_1}{\rho}\right]\]

\[\lambda = \lambda_0 + \arctan\left[\frac{x \cdot \sin c}{\rho\cos\phi_1\cos c - y\sin\phi_1\sin c}\right]\]

If \(\rho = 0\), the point is at the projection center: \(\phi = \phi_1\), \(\lambda = \lambda_0\).

Limitations

Antipode singularity: The point diametrically opposite the projection center (angular distance 180°) maps to infinity and must be excluded from the domain.

Shape distortion: While area is preserved exactly, shapes distort increasingly with distance from center. The maximum angular distortion \(\omega\) at distance \(c\) is (Snyder, 1987, eq. 24-24, p. 188):

\[\sin\left(\frac{\omega}{2}\right) = \frac{k'^2 - 1}{k'^2 + 1}\]

At \(c = 90°\), distortion reaches \(\omega \approx 70.5°\). ISEA grids limit this by using icosahedral faces subtending only ~72° from face center.

No conformality: No projection can be both equal-area and conformal (angle-preserving)—a fundamental constraint from differential geometry (Snyder, 1987, p. 16-18).

Geometry

A regular icosahedron has:

• 20 equilateral triangular faces (each covering ~1/20 of Earth’s surface)

• 12 vertices (where 5 faces meet—these become pentagon cells)

• 30 edges

Left: the icosahedron has the same surface area as the sphere, so its face plane lies at the inradius \(R_1 = 0.9104R\) from the centre O, inside the sphere, and the vertex V of the planar face lies on the ray through the sphere point at angular distance \(g = 37.38°\) from the face centre C. In this cross-section the face extends from V on one side to the midpoint of the opposite edge on the other. A point P at angular distance \(z\) from C maps to P’ on the face plane at radial distance \(\rho\) from C. Right: on the planar face the azimuth Az’ is measured from the centre-to-vertex ray, and the three-fold symmetry of the face reduces every azimuth to one 120° sector.

Vertex Latitude Derivation

The 12 vertices are located at (Coxeter, 1973, p. 52-53):

The latitude \(\arctan(1/2)\) arises from the golden ratio geometry. An icosahedron can be constructed from three mutually perpendicular golden rectangles (\(1 \times \varphi\), where \(\varphi = (1+\sqrt{5})/2\)). The non-polar vertices have \(z\)-coordinate \(1/s\) where \(s = \sqrt{1 + \varphi^2}\), yielding \(\tan\phi = 1/2\).

Standard ISEA Orientation

The default orientation places vertex 0 at longitude 11.25°, latitude \(\arctan(\varphi) \approx 58.2825°\), with azimuth 0°. This places icosahedron vertices (pentagon cells) predominantly over oceans (Sahr et al., 2003, p. 123).

Key Constants

The vertex angle has the closed form \(\tan E_l = \sqrt{14 - 6\sqrt{5}}\) (Snyder, 1992, p. 17). \(R_1\) is the radius of the icosahedron whose surface area equals that of the unit sphere: equating \(20 \cdot \frac{3\sqrt{3}}{4} R_1^2 \tan^2 E_l = 4\pi\) gives \(R_1 = \frac{2}{\tan E_l}\sqrt{\frac{\pi}{15\sqrt{3}}}\).

Forward Projection Steps

The complete algorithm comprises seven steps (Snyder, 1992, p. 12-13):

Step 1: Compute angular distance and azimuth from face center \((\lambda_0, \phi_0)\) to point \((\lambda, \phi)\):

\[z = \arccos(\sin\phi_0 \sin\phi + \cos\phi_0 \cos\phi \cos(\lambda - \lambda_0))\] \[\text{Az} = \arctan2(\cos\phi \sin(\lambda - \lambda_0), \cos\phi_0 \sin\phi - \sin\phi_0 \cos\phi \cos(\lambda - \lambda_0))\]

Step 2: Reduce azimuth to [0°, 120°) by exploiting 3-fold symmetry.

Step 3: Compute auxiliary angle \(\delta_z\) (Snyder, 1992, eq. 9, p. 13): \[\delta_z = \arctan\left(\frac{\tan E_l}{\cos \text{Az} + \cot 30° \cdot \sin \text{Az}}\right)\]

Step 4: Compute auxiliary angle \(h\) (Snyder, 1992, eq. 6, p. 12): \[h = \arccos(\sin \text{Az} \sin G \cos E_l - \cos \text{Az} \cos G)\]

Step 5: Compute adjusted azimuth \(\text{Az}'\) (Snyder, 1992, eq. 7-8, p. 12): \[A_G = \text{Az} + G + h - \pi\] \[\text{Az}' = \arctan\left(\frac{2 A_G}{R_1^2 \tan^2 E_l - 2 A_G \cot 30°}\right)\]

Step 6: Compute radial distance (Snyder, 1992, eq. 10-12, p. 13): \[f = \frac{\tan E_l}{2(\cos \text{Az}' + \cot 30° \cdot \sin \text{Az}') \sin(\delta_z / 2)}\] \[\rho = 2 R_1 f \sin(z / 2)\]

Step 7: Convert to Cartesian with sector offset restored.

The Inverse Projection

Snyder (1992, p. 13) inverts the projection by iterating on the azimuth adjustment, and states that it has no closed form. Recht (2021) restates the projection with unit vectors and inverts it in closed form; hexify uses that inverse.

Write the point’s sector as the spherical triangle \((\mathbf{v}_0, \mathbf{v}_1, \mathbf{v}_2)\) (face center, then the sector’s two vertices) and its plane image as the triangle with barycentric coordinates \((\beta_0, \beta_1, \beta_2)\). With \(c_{ij} = \mathbf{v}_i \cdot \mathbf{v}_j\), \(s_{12} = \lVert \mathbf{v}_1 \times \mathbf{v}_2 \rVert\), \(|\mathbf{V}| = \det[\mathbf{v}_0, \mathbf{v}_1, \mathbf{v}_2]\) and \(A\) the spherical area of the sector:

\[h = 1 - \beta_0, \qquad a = \frac{\beta_2}{h} A\]

\[f = \sin a \, |\mathbf{V}| + (1 - \cos a)(c_{01} c_{12} - c_{20}), \qquad g = (1 - \cos a) \, s_{12} (1 + c_{01})\]

The point \(\mathbf{p}\) where the great circle from \(\mathbf{v}_0\) through the target meets the edge \(\mathbf{v}_1 \mathbf{v}_2\) lies at arc \(2 \arctan(g / f)\) from \(\mathbf{v}_1\) along that edge. The target lies on the arc from \(\mathbf{v}_0\) to \(\mathbf{p}\), at the distance \(z\) from the face center given by

\[\sin\frac{z}{2} = h \sin\frac{z_p}{2},\]

where \(z_p\) is the arc from \(\mathbf{v}_0\) to \(\mathbf{p}\). hexify evaluates the arctangent as atan2(g, f) and \(1 - \cos a\) as \(2 \sin^2(a/2)\), which keeps the inverse at full double precision near the face center and near the sector edges. A forward projection followed by the inverse returns to within a few times \(10^{-15}\) radians.

Distortion

The projection is exactly equal-area, so its distortion falls on angles and local scale. On an icosahedron face the maximum angular deformation is \(\omega = 17.27°\), with linear scale factors between 0.860 and 1.163 (Snyder, 1992, Table 1, p. 12). The azimuth adjustment concentrates this deformation along the three radii from the face center to its vertices, where the graticule shows visible cusps (Snyder, 1992, p. 11).

Deformation is smallest, about 5°, in the ring around the face center toward the edge midpoints, and largest along the three vertex radii, where the parallels kink as they cross. The 17.27° maximum is reached at the face center approached along a vertex radius. The Jacobian determinant of the projection equals 1 everywhere on the face, which is the equal-area property in differential form.

The Truncated Icosahedron

Snyder’s construction applied to the truncated icosahedron deforms angles by at most 3.75° on its 20 hexagons and 2.65° on its 12 pentagons (Snyder, 1992, Table 1, p. 12). Cutting each corner of the icosahedron at a third of the edge places every hexagon in an icosahedron face plane, centred on the face centre (the cut keeps the face’s three-fold symmetry), and every pentagon in a plane normal to a vertex axis, centred on the vertex (five-fold symmetry). A spherical icosahedron face is therefore the spherical hexagon plus a fifth of each of the three pentagons at its corners, and the solid’s face centres are the centres of the icosahedral grid’s faces and pentagon cells. The question is whether the hexagonal cells, indexed on the icosahedron’s plane triangle as now, could be projected through the truncated icosahedron’s faces instead. Three properties of that composite settle it.

Area shares. Snyder’s right triangles have areas 0.0376062647 on the hexagon and 0.0295072264 on the pentagon (unit sphere; Snyder, 1992, eqs. 2 and 24), so the spherical hexagon covers 0.71823 of a face and each pentagon fifth 0.09392. The plane hexagon cut at a third covers 2/3 of the plane triangle and each corner 1/9. An equal-area map cannot send the solid’s vertices to the third points of the plane edges with straight sides between them. A slice map (Snyder’s construction from a polygon centre) followed by an affine map of its plane triangle is again a slice map, so the only way to carry the construction into the plane triangle with straight pieces, three-fold symmetry and exact area is to slice the 18 spherical right triangles of a face onto 18 plane triangles of proportional area. Their corners are then fixed: the solid’s vertices at 0.32044 of the edge from each corner, and the midpoint of each hexagon-pentagon edge at 0.25384 of the edge along the bisector, so the hexagon-pentagon boundary bends at the bisector. The pentagons take the hexagons’ area scale here, where Snyder’s projection gives them a 1.28 percent larger one (Snyder, 1992, p. 11).

Seams. A slice map spaces the points of the edge opposite its apex by area. On the right spherical triangle with legs \(h\) (apex to the foot \(M\) of its perpendicular to the edge) and \(x\) (from \(M\) along the edge), the area is \(2\arctan(\tan\frac{h}{2}\tan\frac{x}{2})\), so a point at \(x\) goes to the share \(\arctan(\tan\frac{h}{2}\tan\frac{x}{2}) / \arctan(\tan\frac{h}{2}\tan\frac{X}{2})\) of the half-edge \(X\). This share grows strictly with \(h\) at every interior \(x\). A hexagon-pentagon edge is perpendicular to the arc from the hexagon’s centre to the pentagon’s, so both apexes have \(M\) as their foot, at different distances: the two sides send the same point to two places. The offset reaches \(1.74 \times 10^{-5}\) of the half-edge, 21.7 m on the Earth, which is the offset Snyder reports for his own projection (“less than 0.000018 of length BC”, Snyder, 1992, p. 16). On Snyder’s icosahedron and on the vertex-oriented equal-area projection (projection = "ivea"; van Leeuwen and Strebe, 2006) every edge inside the map joins two mirror images, which space it alike. A continuous equal-area composite would have to bend the images of the slicing arcs, and no closed form for it is known.

A bound at the vertices. Five plane triangles meet at each vertex of the face frame, with 300° about it. Measure the frame in units in which its area equals the sphere’s. A circle of radius \(\rho\) about the vertex then has length \(\frac{5\pi}{3}\rho\) and encloses the area \(A = \frac{5\pi}{6}\rho^2\). Under an equal-area projection its preimage on the unit sphere is a closed curve enclosing the same area, so its length \(\ell\) satisfies the spherical isoperimetric inequality \(\ell^2 \ge A(4\pi - A)\). Where the projection has scale factors \(a \ge b\) with \(ab = 1\), the inverse stretches lengths by at most \(a = \sqrt{a/b}\). Hence some point of the circle has

\[\frac{a}{b} \ge \frac{\ell^2}{(5\pi\rho/3)^2} \ge \frac{6}{5} - \frac{\rho^2}{4}, \qquad \omega = 2\arcsin\frac{a/b - 1}{a/b + 1} \ \to\ 2\arcsin\frac{1}{11} = 10.43° \quad (\rho \to 0).\]

Every equal-area projection onto the icosahedron’s faces has points with \(\omega \ge 10.43°\) arbitrarily close to each vertex, whatever solid it passes through on the way. The truncated icosahedron’s 2.65° is measured on its own flat pentagon, which has 360° about its centre; in the icosahedron’s frame those 360° close to 300°.

Measured over a face (paper/bench/bench_trunc_icosa.R; 100,000 points on one face for the composite, 2,000,000 on the sphere for the others, and a cap of 0.5° about a vertex):

The vertex-oriented projection is continuous, has a closed-form inverse, and has smaller mean and maximum deformation than the truncated-icosahedron composite. Its maximum, reached at the vertices, lies 0.68° above the bound, which caps what any equal-area projection onto the icosahedron’s faces could gain there. hexify therefore offers no truncated-icosahedron projection.

The Earth as an Ellipsoid

Snyder’s projection is equal-area on the sphere. By default hexify, like DGGRID and H3, reads a geodetic latitude \(\varphi\) as a latitude on the sphere whose area is the WGS84 ellipsoid’s. A small region at \(\varphi\) then has \(MN/R_q^2\) times its spherical area on the ellipsoid, with \(M\) and \(N\) the meridional and prime-vertical radii of curvature and \(R_q\) the radius of that sphere. Measured geodesically on WGS84, the cells’ areas run from 0.9955 times the reported area at the equator to 1.0090 times it at the poles (paper/bench/bench_ellipsoid_area.R).

The authalic latitude \(\beta\) removes the factor. With \(e\) the eccentricity,

\[q(\varphi) = (1 - e^2)\left[\frac{\sin\varphi}{1 - e^2\sin^2\varphi} + \frac{\operatorname{artanh}(e \sin\varphi)}{e}\right], \qquad \sin\beta = \frac{q(\varphi)}{q(\pi/2)}, \qquad R_q = a\sqrt{q(\pi/2)/2},\]

the ellipsoid below \(\varphi\) holds as much area as the sphere of radius \(R_q\) below \(\beta\), so \((\lambda, \varphi) \mapsto (\lambda, \beta)\) keeps every area, and cells equal-area on that sphere are equal-area on the ellipsoid. hex_grid(ellipsoid = "WGS84") converts geodetic latitude to \(\beta\) where a point enters the projection and back where it leaves, and sizes the grid on \(R_q\) (6371.0072 km). hexify evaluates both directions by a sine series in \(2\varphi\) (or \(2\beta\)) whose coefficients it reads off the closed form above when the grid is used; on WGS84 a latitude moves by up to 0.128°. The cells’ geodesic areas then agree with the reported ones to the measurement’s own error, about \(2 \times 10^{-6}\).

The ISEA3H, ISEA7H, IVEA3H and IVEA7H definitions OGC registers are built this way, with the first vertex at 11.20°E and authalic latitude \(\arctan\phi\) (orientation = "ogc"), which is geodetic latitude 58.397145907431° on WGS84:

ogc <- hex_grid(resolution = 8, aperture = 3, ellipsoid = "WGS84",
                orientation = "ogc")
ogc
#> <HexGridInfo> aperture 3, resolution 8, WGS84 ellipsoid: 65,612 cells of 7,774 km^2

The Unfolded Projection

Applying the per-face projection to all 20 faces and unfolding them into a plane maps the whole sphere onto a strip of triangles (compare Snyder, 1992, Figure 12, p. 18). hexify_lonlat_to_plane() performs the composite transformation, in the same PLANE layout as dggridR’s dgGEO_to_PLANE(). Coastlines and graticule continue across every shared edge and are interrupted elsewhere:

draw_path <- function(lon, lat, col, lwd) {
  p <- hexify_lonlat_to_plane(lon, lat)
  jump <- c(FALSE, sqrt(diff(p$plane_x)^2 + diff(p$plane_y)^2) > 0.2)
  p$plane_x[jump] <- NA
  lines(p$plane_x, p$plane_y, col = col, lwd = lwd)
}

oldpar <- par(mar = c(0.5, 0.5, 2, 0.5), bg = "white")
plot(NULL, xlim = c(-0.1, 5.6), ylim = c(-0.05, 2.65), asp = 1,
     axes = FALSE, xlab = "", ylab = "", main = "Unfolded ISEA Projection")

# face outlines: side-1 triangles around each projected face centroid
centers <- hexify_face_centers()
pc <- hexify_lonlat_to_plane(centers$lon * 180/pi, centers$lat * 180/pi)
s3 <- sqrt(3)
for (f in 1:20) {
  base <- floor(pc$plane_y[f] / (s3/2)) * s3/2
  up <- (pc$plane_y[f] - base) < s3/4
  yy <- if (up) base + c(0, 0, s3/2) else base + c(s3/2, s3/2, 0)
  polygon(pc$plane_x[f] + c(-0.5, 0.5, 0), yy, border = "gray60", lwd = 0.8)
}

# 15-degree graticule
for (lon in seq(-180, 165, by = 15))
  draw_path(rep(lon, 181), seq(-90, 90, by = 1), "gray75", 0.6)
for (lat in seq(-75, 75, by = 15))
  draw_path(seq(-180, 180, by = 0.5), rep(lat, 721), "gray75", 0.6)

# coastlines
xy <- st_coordinates(st_cast(st_geometry(hexify_world), "MULTIPOLYGON"))
for (ring in split(seq_len(nrow(xy)),
                   interaction(xy[, "L1"], xy[, "L2"], xy[, "L3"], drop = TRUE)))
  draw_path(xy[ring, 1], xy[ring, 2], GREY, 0.8)

par(oldpar)

Aperture Properties

The aperture 7 rotation angle follows from the child lattice: the parent’s nearest neighbour sits at two steps along one child lattice axis plus one step along the next, a vector of length \(\sqrt{7}\) child spacings at angle \(\arctan\!\big(\tfrac{\sqrt{3}/2}{2 + 1/2}\big) = \arctan(\sqrt{3}/5) \approx 19.107°\) to the axis (DGGRID Manual, 2023).

How a Cell Divides Among Its Parents

Since successive resolutions do not nest, a cell can lie partly in two or three coarser cells (Carr et al., 1997). Within a face the cells of one resolution are the Voronoi cells of a hexagonal lattice and those one resolution coarser the Voronoi cells of a sublattice of index 3, 4 or 7, so the share of a cell in each coarser cell follows from the lattices alone.

After an aperture-3 step the children are centred on the parent centres and corners. The three parent edges meeting at a corner run through every second edge midpoint of the child centred there and cut it into thirds. After an aperture-4 step the children are centred on the parent centres and edge midpoints, and a parent edge runs through the centre and two opposite corners of the child on its midpoint, cutting it in half.

After an aperture-7 step every child centre lies inside a parent. Take the finer lattice as the Eisenstein integers \(\mathbb{Z}[\omega]\), \(\omega = e^{2\pi i/3}\), with unit spacing, so a cell has circumradius \(s = 1/\sqrt{3}\) and area \(\tfrac{3\sqrt{3}}{2}s^2\). The coarser lattice is \(\mu\,\mathbb{Z}[\omega]\) with \(\mu = 3 + \omega\): \(|\mu|^2 = 7\) and \(\arg \mu = \arctan(\sqrt{3}/5)\), the turn in the table above, and the other twist is its mirror image \(3 + \bar\omega\). The centre child, of circumradius \(s \approx 0.577\), lies inside its parent, of inradius \(\sqrt{7}/2 \approx 1.323\). The parent edge facing the neighbour at \(\mu\) runs between the parent corners \(\mu(1-\omega)/3\) and \(\mu(2+\omega)/3\), through its midpoint \(\mu/2\), and two of these points belong to the child tiling:

\[\frac{\mu(1-\omega)}{3} = 1 + \frac{1-\omega}{3}, \qquad \frac{\mu}{2} = \frac{1 + (2+\omega)}{2}.\]

The first is a corner \(A\) of the ring child centred at \(1\); the second is the midpoint of the edge this child shares with the child at \(2 + \omega\), the second edge from \(A\). The child’s other corners lie inside the parent, so the part beyond it is the triangle of \(A\), the next corner \(B\) and \(\mu/2\), whose sides \(s\) and \(s/2\) meet at \(120°\):

\[\frac{\tfrac{1}{2} \cdot s \cdot \tfrac{s}{2} \cdot \sin 120°}{\tfrac{3\sqrt{3}}{2}\, s^2} = \frac{\sqrt{3}\, s^2 / 8}{3\sqrt{3}\, s^2 / 2} = \frac{1}{12}.\]

The half-turn about \(\mu/2\) maps both lattices onto themselves and swaps the parents at \(0\) and \(\mu\), so the other half of the edge cuts \(1/12\) off the child at \(2 + \omega\) the other way. Every ring child lies \(11/12\) in its parent and \(1/12\) in one neighbour, and the parent’s area checks: \(1 + 6 \cdot \tfrac{11}{12} + 6 \cdot \tfrac{1}{12} = 7\) child areas.

The vertices of the solid are cell centres at every resolution, and the turn by the angle deficit about a vertex maps both lattices onto themselves, so the cells around a vertex are cut as in the plane and the cell at a vertex lies wholly in the vertex cell above it. Snyder’s projection is equal-area, so these shares hold on the sphere for every cell, across face edges and at the vertices. get_parent(overlapping = TRUE) returns the coarser cells a cell reaches into, and hex_aggregate() divides each cell’s value among them by these shares, one level at a time, each level by the shares of its own aperture:

g7 <- hex_grid(resolution = 3, aperture = 7)
ids <- grid_global(g7)$cell_id
table(lengths(get_parent(ids, g7, overlapping = TRUE)))
#> 
#>    1    2 
#>  492 2940

# A uniform density adds up to each parent's area, the twelve pentagons included
up <- hex_aggregate(ids, cell_area(ids, g7), g7)
range(up$value / cell_area(up$cell_id, hex_grid(resolution = 2, aperture = 7)))
#> [1] 1 1

Cell Count Growth

cat("Resolution  Aperture 3    Aperture 4    Aperture 7\n")
#> Resolution  Aperture 3    Aperture 4    Aperture 7
cat("---------  ----------    ----------    ----------\n")
#> ---------  ----------    ----------    ----------
for (res in 0:8) {
  cells_ap3 <- 10 * 3^res + 2
  cells_ap4 <- 10 * 4^res + 2
  cells_ap7 <- 10 * 7^res + 2
  cat(sprintf("    %d      %10s    %10s    %10s\n",
              res,
              format(cells_ap3, big.mark = ","),
              format(cells_ap4, big.mark = ","),
              format(cells_ap7, big.mark = ",")))
}
#>     0              12            12            12
#>     1              32            42            72
#>     2              92           162           492
#>     3             272           642         3,432
#>     4             812         2,562        24,012
#>     5           2,432        10,242       168,072
#>     6           7,292        40,962     1,176,492
#>     7          21,872       163,842     8,235,432
#>     8          65,612       655,362    57,648,012

Triangle to Quad

The 20 triangular faces are paired into 12 “quads” (diamond-shaped regions). Each quad contains two adjacent triangular faces sharing an edge. This pairing transforms 20 triangles into 10 diamond-shaped quads plus 2 polar quads, simplifying grid indexing because a diamond admits a rectangular \((i, j)\) lattice (DGGRID Manual, 2023).

Quad IJ: The Integer Lattice

After Snyder projection maps a point onto a triangular face, the continuous \((x, y)\) coordinates are quantized to integer \((i, j)\) indices on the hexagonal lattice. At resolution \(r\) with aperture \(a\), the lattice dimension along each axis is:

\[d = \left\lfloor a^{r/2} \right\rfloor\]

The \((i, j)\) pair uniquely identifies a cell within a quad. Combined with the quad number (0–11), this gives a globally unique cell address.

Hierarchical Index Types

hexify supports three hierarchical index encodings. Each converts the \((quad, i, j)\) triple into a compact string or integer that encodes the cell’s position in the subdivision hierarchy.

# Z7 index encoding for aperture 7
g7 <- hex_grid(resolution = 4, aperture = 7)
cell <- lonlat_to_cell(16.37, 48.21, g7)
idx <- cell_to_index(cell, g7)
cat(sprintf("Cell %s -> Z7 index: %s\n", as.character(cell), idx))
#> Cell 5194 -> Z7 index: 000324
cat(sprintf("  Leading field: %s (quad %d), Digits: %s\n",
            substr(idx, 1, 2), as.integer(substr(idx, 1, 2)) %% 12L,
            substr(idx, 3, nchar(idx))))
#>   Leading field: 00 (quad 0), Digits: 0324

# Hierarchical property: parent is obtained by dropping the last digit
parent_idx <- substr(idx, 1, nchar(idx) - 1)
parent_info <- hexify_index_to_cell(parent_idx, 7, "z7")
cat(sprintf("  Parent index: %s (face %d, i=%d, j=%d)\n",
            parent_idx, parent_info$face,
            as.integer(parent_info$i), as.integer(parent_info$j)))
#>   Parent index: 00032 (face 3, i=9, j=19)

# A valid Z7 index decodes and re-encodes without changing
decoded <- hexify_index_to_cell(idx, 7, "z7")
stopifnot(identical(
  hexify_cell_to_index(decoded$face, decoded$i, decoded$j,
                       decoded$resolution, 7, "z7"),
  idx
))

g7 <- hex_grid(resolution = 4, aperture = 7)
cell <- index_to_cell("033000", g7)
cell_to_index(cell, g7, form = "hex")
#> [1] "3600ffffffffffff"
cell_to_index(cell, g7, form = "int")
#> integer64
#> [1] 3891391553024819199
cell_to_index(cell, g7, form = "monotonic")
#> integer64
#> [1] 8232

cell <- index_to_cell("E6-317-A", ogc, form = "textZIRS")
cell_to_lonlat(cell, ogc)
#>    lon_deg  lat_deg
#> 1 34.78017 45.42938
cell_to_index(cell, ogc, form = "uint64ZIRS")
#> integer64
#> [1] 630503947831872604

# Z3 index encoding for aperture 3
g3 <- hex_grid(resolution = 8, aperture = 3)
cell <- lonlat_to_cell(16.37, 48.21, g3)
idx <- cell_to_index(cell, g3)
cat(sprintf("Cell %s -> Z3 index: %s\n", as.character(cell), idx))
#> Cell 14092 -> Z3 index: 0321202211
cat(sprintf("  Base cell: %s, Digits: %s (%d digit pairs)\n",
            substr(idx, 1, 2), substr(idx, 3, nchar(idx)),
            (nchar(idx) - 2) / 2))
#>   Base cell: 03, Digits: 21202211 (4 digit pairs)

# Z-order index for aperture 4
g4 <- hex_grid(resolution = 8, aperture = 4)
cell <- lonlat_to_cell(16.37, 48.21, g4)
idx <- cell_to_index(cell, g4)
cat(sprintf("Aperture 4: Cell %s -> Z-order index: %s\n", as.character(cell), idx))
#> Aperture 4: Cell 140534 -> Z-order index: 0311310300

Index Type Comparison

All three hexify encodings are bijective: each cell maps to exactly one index string, and each valid string to one cell. The strings hold every resolution; on aperture-7 grids of the icosahedron the Z7 index is also available as IGEO7’s 64-bit integer up to resolution 20 (see above).

SEQNUM: The Flat Cell ID

The SEQNUM provides a unique integer for each cell, independent of the hierarchical index. The total cell count at resolution \(r\) with aperture \(a\) is:

\[N(r) = 10 \times a^r + 2\]

The “+2” accounts for the two polar pentagon cells. SEQNUMs are assigned by traversing quads in order (0–11) and cells within each quad in row-major \((i, j)\) order. This numbering maintains compatibility with dggridR.

Compact uint64 Storage

For database-friendly storage, hexify can pack any hierarchical index into a single 64-bit unsigned integer:

\[\texttt{uint64} = (\texttt{face} \ll 60) \;|\; \texttt{packed\_digits}\]

The top 4 bits encode the face (0–11), and the remaining 60 bits pack the index digits. This supports resolutions up to \(\lfloor 60 / \lceil\log_2 a\rceil \rfloor\) before overflow.

Architecture

Both systems tile the sphere with hexagons on an icosahedral framework, but the projection and subdivision choices lead to different properties.

The Area Variation Trade-Off

The most significant difference is area uniformity. ISEA uses the Snyder equal-area projection, guaranteeing that all hexagonal cells at a given resolution have identical area (pentagons are 5/6 of hex area). H3 uses a gnomonic (central) projection per face, which is not equal-area. This causes cell areas to vary by latitude:

# Compare ISEA (constant area) vs H3 (variable area) at similar resolutions
lats <- seq(-85, 85, by = 5)
lons <- rep(10, length(lats))

# ISEA: aperture 7, resolution 6 (~130 km² cells)
g_isea <- hex_grid(resolution = 6, aperture = 7)
isea_cells <- lonlat_to_cell(lons, lats, g_isea)
isea_areas <- cell_area(isea_cells, g_isea)

# H3: resolution 4 (~1,770 km² cells — different scale, but shows the pattern)
g_h3 <- hex_grid(resolution = 4, type = "h3")
h3_cells <- lonlat_to_cell(lons, lats, g_h3)
h3_areas <- cell_area(h3_cells, g_h3)

# Normalize to show relative variation
isea_rel <- isea_areas / mean(isea_areas)
h3_rel <- h3_areas / mean(h3_areas)

oldpar <- par(mar = c(4, 4, 2, 1), bg = "white")
plot(lats, h3_rel, type = "l", col = "#E63946", lwd = 2.5,
     xlab = "Latitude (degrees)", ylab = "Relative cell area",
     main = "Cell Area Variation by Latitude",
     ylim = range(c(isea_rel, h3_rel)))
lines(lats, isea_rel, col = "#457B9D", lwd = 2.5)
abline(h = 1, lty = 2, col = GREY)
legend("bottomleft", legend = c("H3 (gnomonic)", "ISEA (equal-area)"),
       col = c("#E63946", "#457B9D"), lwd = 2.5, bty = "n")

par(oldpar)

For ecological and statistical applications where equal-area cells are important (species density estimation, spatial sampling), ISEA is the correct choice. For applications where hierarchical indexing speed matters more than area equality (ride-sharing, logistics), H3 may be preferable.

Resolution Mapping

Because ISEA supports multiple apertures, there is no one-to-one resolution mapping to H3. The h3_crosswalk() function finds the closest H3 resolution for a given ISEA grid:

# H3 resolution table: compare H3 and ISEA aperture-7 cell areas
cat("H3 Res  Avg Area (km²)  ISEA Ap7 Equivalent\n")
#> H3 Res  Avg Area (km²)  ISEA Ap7 Equivalent
cat("------  --------------  -------------------\n")
#> ------  --------------  -------------------
for (h3_res in 0:8) {
  g_h3 <- hex_grid(resolution = h3_res, type = "h3")
  h3_area <- summary(g_h3)$area_km2

  # Find closest ISEA ap7 resolution by brute force
  best_res <- 0
  best_diff <- Inf
  for (r in 0:15) {
    g_test <- hex_grid(resolution = r, aperture = 7)
    d <- abs(log(summary(g_test)$area_km2) - log(h3_area))
    if (d < best_diff) { best_diff <- d; best_res <- r }
  }
  g_isea <- hex_grid(resolution = best_res, aperture = 7)
  cat(sprintf("   %2d   %14.1f  res %d (%.1f km²)\n",
              h3_res, h3_area, best_res, summary(g_isea)$area_km2))
}
#>     0        4357449.4  res 1 (7084244.7 km²)
#>     1         609788.4  res 2 (1036718.7 km²)
#>     2          86801.8  res 3 (148620.5 km²)
#>     3          12393.4  res 4 (21242.1 km²)
#>     4           1770.3  res 5 (3034.8 km²)
#>     5            252.9  res 6 (433.5 km²)
#>     6             36.1  res 7 (61.9 km²)
#>     7              5.2  res 8 (8.8 km²)
#>     8              0.7  res 9 (1.3 km²)

When to Use Which