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. After Lambert projection, shapes change (outer bands stretch radially, compress tangentially) but areas remain equal.

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

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).

#> Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
#> give correct results for longitude/latitude data

Key Constants

Forward Projection Steps

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

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. 8, p. 14): \[\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. 9, p. 14): \[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. 10-11, p. 14): \[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. 12-13, p. 14-15): \[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

The inverse projection cannot be solved analytically because the azimuth adjustment contains transcendental functions. A Newton-Raphson iteration finds the spherical azimuth Az from the planar azimuth Az’:

\[f(\text{Az}) = \text{agh} - \text{Az} - G + (\pi - h) = 0\]

where \(h = \arccos(\sin \text{Az} \sin G \cos E_l - \cos \text{Az} \cos G)\).

Newton-Raphson iteration for inverse projection

The iteration exhibits quadratic convergence, typically reaching machine precision in 3-5 iterations. hexify provides four precision modes:

Aperture Properties

The aperture 7 rotation angle \(\arctan(\sqrt{3/7})\) arises from the geometric constraint that 7 hexagons in a rosette pattern (1 center + 6 ring) must maintain lattice consistency (DGGRID Manual, 2023).

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 %d -> Z7 index: %s\n", cell, idx))
#> Cell 5194 -> Z7 index: 270453
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: 27 (quad 3), Digits: 0453

# 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: 27045 (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
))

# 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 %d -> Z3 index: %s\n", 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 %d -> Z-order index: %s\n", cell, idx))
#> Aperture 4: Cell 140534 -> Z-order index: 0311310300

Index Type Comparison

All three hexify encodings are bijective: each valid \((quad, i, j)\) cell maps to exactly one index string, and vice versa. hexify stores indices as character strings to support arbitrary precision and avoid integer overflow at high resolutions. As noted above, hexify’s bijective Z7 strings intentionally differ from DGGRID in the pentagon regions where DGGRID’s raw Z7 encoding collides.

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 = "gray50")
legend("topright", 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 <- 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(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, 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