H3 Grid System

H3 Grid System

Code 
using GlobalGrids, GeoMakie, CairoMakie

grid = H3Grid()

fig = Figure()
ax = GeoAxis(fig[1,1], dest = "+proj=ortho +lon_0=-50 +lat_0=30")
lines!(ax, [grid[i] for i in 1:121])
poly!(ax, collect(GlobalGrids.pentagons(grid, 0)), color=(:black, .4))

ep = surface!(ax,
    -180..180, -90..90,
    zeros(axes(rotr90(GeoMakie.earth())));
    shading = NoShading, color = rotr90(GeoMakie.earth())
)
translate!(ep, 0, 0, -1)

fig
Precompiling packages...
   1005.2 ms  ✓ QuartoNotebookWorkerTablesExt (serial)
  1 dependency successfully precompiled in 1 seconds
Precompiling packages...
   1863.2 ms  ✓ QuartoNotebookWorkerLaTeXStringsExt (serial)
  1 dependency successfully precompiled in 2 seconds
Precompiling packages...
   1175.1 ms  ✓ QuartoNotebookWorkerJSONExt (serial)
  1 dependency successfully precompiled in 1 seconds
Precompiling packages...
   6530.1 ms  ✓ QuartoNotebookWorkerMakieExt (serial)
  1 dependency successfully precompiled in 7 seconds
Precompiling packages...
   1163.2 ms  ✓ QuartoNotebookWorkerJSON3Ext (serial)
  1 dependency successfully precompiled in 1 seconds
Precompiling packages...
   5683.4 ms  ✓ QuartoNotebookWorkerCairoMakieExt (serial)
  1 dependency successfully precompiled in 6 seconds

H3 is a DGGS developed by Uber. Here’s a brief summary of the H3 grid system:

  • There are 16 resolutions (0-15 coarse-to-fine).
  • At each resolution, there are exactly 12 pentagon cells in the grid. Every other cell is a hexagon.
  • Pentagon locations are determined such that their centers are not on land.
  • Refer to https://h3geo.org/docs for more information.

Creating an H3Cell

  • Every cell in the grid is indexed by a single UInt64. An H3Cell is simply a wrapper around a UInt64 index.
  • The index encodes the base cell (121 possible values, shown in the figure above), resolution, and digits (0-6).
  • Each hexagon gets split into 7 hexagons by the next highest resolution. An additional digit corresponds to one of the smaller hexagons derived from its parent.

From H3Grid

grid = H3Grid()

# Select one of the 121 base cells:
grid[121]
  ⬡ H3Cell 0 121-
# Nesting getindex to dive down the hierarchy:
grid[121][1][5][0]
  ⬡ H3Cell 3 121-150

From Base Cell and Digits

H3Cell(121, [1, 5, 0])
  ⬡ H3Cell 3 121-150

From Coordinates and Resolution:

coords = LonLat(-75.0, 54.0)

H3Cell(coords, 10)
  ⬡ H3Cell 10 7-1203346514

From Hexadecimal String

c = H3Cell(coords)

str = GlobalGrids.h3_string(c)
@info str

H3Cell(str)
[ Info: 8a0e506e6a67fff

  ⬡ H3Cell 10 7-1203346514

Inspecting H3Cells

Code 
@show resolution(c)
@show base_cell(c)
@show is_pentagon(c)
@show GlobalGrids.digits(c)
@show encode(c)
@show decode(c)
@show GlobalGrids.h3_face_numbers(c)
nothing
resolution(c) = 10
base_cell(c) = 7
is_pentagon(c) = false
GlobalGrids.digits(c) = [1, 2, 0, 3, 3, 4, 6, 5, 1, 4]
encode(c) = "8a0e506e6a67fff"
decode(c) = "7-1203346514"
GlobalGrids.h3_face_numbers(c) = Int32[2]

GeoInterface Integration:

Code 
import GeoInterface as GI

@show GI.coordinates(c)
@show GI.centroid(c)
@show GI.area(c)  # meters ^ 2
nothing
GI.coordinates(c) = GlobalGrids.LonLat{Float64}[LonLat{Float64} (-74.99962251383738, 54.0001926668228), LonLat{Float64} (-75.0008585102879, 54.0000194052262), LonLat{Float64} (-75.00122192991167, 53.999307797036366), LonLat{Float64} (-75.00034938834023, 53.99876945809693), LonLat{Float64} (-74.99911343000446, 53.99894271390214), LonLat{Float64} (-74.99874997512619, 53.99965431443802), LonLat{Float64} (-74.99962251383738, 54.0001926668228)]
GI.centroid(c) = LonLat{Float64} (-74.99998595778595, 53.99948106205936)
GI.area(c) = 17803.420253871413

Grid Hierarchy

  • parent, siblings, and children
Code 
fig = Figure()
ax = Axis(fig[1,1], title="Cell:Black, Parent:Blue, Siblings:Orange, Children:Green")
lines!(ax, GlobalGrids.parent(c), color = :blue)
poly!(ax, c, color = :black)
lines!(ax, GlobalGrids.siblings(c), color = :orange)
lines!(ax, GlobalGrids.children(c), color = :green)
fig

Rings and Disks

  • grid_ring(cell, i) returns the cells that are exactly grid distance i away from the cell.
  • grid_disk(cell, i) is similar to grid_ring, but returns all cells inside the ring as well.
Code 
fig = Figure()
ax = Axis(fig[1,1])
colors = [:blue, :green, :red, :magenta]

poly!(ax, c, color=:black)
for i in 1:4
    poly!(ax, GlobalGrids.grid_ring(c, i), color=colors[i])
end


fig

Quantizing Geometries to the H3 Grid

Point

Code 
ll = LonLat(-54.0, -12.0)

fig = Figure()
ax = Axis(fig[1,1], title="Cells Containing Point at Resolutions 0-3")
for res in 0:3
    lines!(ax, h3cells(ll, res))
end
scatter!(ax, ll)

fig

MultiPoint

Code 
coords = [LonLat(360rand() - 180, 180rand() - 90) for _ in 1:100]

x = H3Cell.(coords, 0)

f, a, p = lines(x, axis=(;type=GeoAxis, dest = "+proj=ortho +lon_0=-50 +lat_0=30"))
scatter!(a, coords, color=:black)
f

Line

  • By default, the cells will be an arbitrary shortest path between the cells that contain the two points (containment = :shortest_path).
  • The other containment options for line(string)s are :overlap and :overlap_bbox.
Code 
line = GI.Line([LonLat(0, 0), LonLat(1,1)])

fig = Figure(size=(1000, 300))

for (i, containment) in enumerate((:shortest_path, :overlap, :overlap_bbox))
    ax = Axis(fig[1, i], aspect=DataAspect(), title="containment = :$containment")
    x = h3cells(line, 6; containment)
    lines!(ax, x)
    lines!(ax, GI.coordinates(line), linewidth=3, color = :black)
end

fig

(Multi)Polygon

  • Containment modes for polygons:
    • :center: Include cell if its center is inside the polygon.
    • :full: Include cell if all vertices are inside the polygon.
    • :overlap: Include cell if any part intersects the polygon.
    • :overlap_bbox: Include cell if any part of the bounding box intersects the polygon.
Code 
using OSMGeocoder

obj = geocode(state="North Carolina")

# Create figure
function make_axis(i, j, containment)
    ax = Axis(fig[i, j], title="containment = :$containment")
    x = h3cells(obj.geometry[1], 4; containment)
    lines!(ax, x)
    lines!(ax, obj.geometry)
end
fig = Figure()
make_axis(1, 1, :center)
make_axis(1, 2, :full)
make_axis(2, 1, :overlap)
make_axis(2, 2, :overlap_bbox)
fig

Extents

Note

The same containment options for (multi)polygons are available for Extents. test

Code 
ex = GI.extent(obj)

# for plotting
x1, x2 = ex.X
y1, y2 = ex.Y
linestring = GI.LineString([(x1,y1), (x1,y2), (x2,y2), (x2, y1), (x1,y1)])

# Create figure
function make_axis(i, j, containment)
    ax = Axis(fig[i, j], title="containment = :$containment")
    x = h3cells(ex, 4; containment)
    lines!(ax, linestring)
    lines!(ax, x)
end
fig = Figure()
make_axis(1, 1, :center)
make_axis(1, 2, :full)
make_axis(2, 1, :overlap)
make_axis(2, 2, :overlap_bbox)
fig

Quantizing Rasters

Code 
import Scratch
ENV["RASTERDATASOURCES_PATH"] = Scratch.get_scratch!("RASTERDATASOURCES_PATH")
nothing
Code 
using GlobalGrids, Rasters, RasterDataSources, ArchGDAL

r = Raster(getraster(WorldClim{Elevation}).elev)

dict = h3cells(r, 2)

cells = collect(keys(dict))
vals = maximum.(values(dict))

# Do some work to avoid polygon issues on boundaries:
idx = findall(x -> !GlobalGrids.crosses_meridian(x, 180), cells)

cells2 = cells[idx]
vals2 = vals[idx]

fig = Figure()
ax1 = Axis(fig[1, 1], title="Without Polygon Filtering")
ax2 = Axis(fig[2, 1], title="With Polygon Filtering")

poly!(ax1, cells, color=vals)
poly!(ax2, cells2, color=vals2)
fig
Note

When plotting large numbers of cells, its more efficient to just plot the centroid of the cell:

Code 
dict = h3cells(r, 4)

cells = collect(keys(dict))

vals = maximum.(values(dict))

scatter(GI.centroid.(cells), color=vals, markersize=2)