2  H3 Grid System

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...
   2076.2 ms  ✓ QuartoNotebookWorkerJSON3Ext (serial)
  1 dependency successfully precompiled in 2 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.

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

2.1.1 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

2.1.2 From Base Cell and Digits

H3Cell(121, [1, 5, 0])

  ⬡ H3Cell 3 121-150

2.1.3 From Coordinates and Resolution:

coords = LonLat(-75.0, 54.0)

H3Cell(coords, 10)

  ⬡ H3Cell 10 7-1203346514

2.1.4 From Hexadecimal String

c = H3Cell(coords)

str = GlobalGrids.h3_string(c)
@info str

H3Cell(str)

[ Info: 8a0e506e6a67fff

  ⬡ H3Cell 10 7-1203346514

2.2 Inspecting H3Cells

@show GlobalGrids.resolution(c)
@show GlobalGrids.h3_base_cell(c)
@show GlobalGrids.is_pentagon(c)
@show GlobalGrids.h3_digits(c)
@show GlobalGrids.h3_face_numbers(c)
nothing

GlobalGrids.resolution(c) = 10
GlobalGrids.h3_base_cell(c) = 7
GlobalGrids.is_pentagon(c) = false
GlobalGrids.h3_digits(c) = [1, 2, 0, 3, 3, 4, 6, 5, 1, 4]
GlobalGrids.h3_face_numbers(c) = Int32[2]

2.3 GeoInterface Integration:

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

2.4 Grid Hierarchy

• parent, siblings, and children

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

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

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

2.6 “Rasterizing” Geometries to the H3 Grid

2.6.1 Point

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

2.6.2 MultiPoint

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

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

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

2.6.4 (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.

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

2.6.5 Extents

Note

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

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

2.7 “Rasterizing” Rasters

• For data that is associated with a location, i.e. rasters, cells creates a Dict{H3Cell, Vector{T}} where T is the eltype of the data.
• Note that plotting polygons across a left/right boundary on a given projection can have strange behavior. See https://en.wikipedia.org/wiki/180th_meridian#Software_representation_problems.

import Scratch
ENV["RASTERDATASOURCES_PATH"] = Scratch.get_scratch!("RASTERDATASOURCES_PATH")
nothing

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

[ Info: Starting download for https://geodata.ucdavis.edu/climate/worldclim/2_1/base/wc2.1_10m_elev.zip
┌ Info: Downloading
│   source = "https://geodata.ucdavis.edu/climate/worldclim/2_1/base/wc2.1_10m_elev.zip"
│   dest = "/home/runner/.julia/scratchspaces/00000000-0000-0000-0000-000000000000/RASTERDATASOURCES_PATH/WorldClim/Elevation/zips/wc2.1_10m_elev.zip"
│   progress = 1.0
│   time_taken = "0.13 s"
│   time_remaining = "0.0 s"
│   average_speed = "9.775 MiB/s"
│   downloaded = "1.271 MiB"
│   remaining = "0 bytes"
└   total = "1.271 MiB"

Note

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

dict = h3cells(r, 4)

cells = collect(keys(dict))

vals = maximum.(values(dict))

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