A round trip through a DGGS store

This tutorial shows how to persist a cell-indexed cube, reopen it lazily, and select only the stored cells needed for a region. It also compares the two cell-id encodings applicable to this regional axis. Complete levels can also use implicit encoding; see Store IO.

dggwrite and dggread are provided by the Zarr.jl extension, loaded by using Zarr.

julia

import DiscreteGlobalGrids as DGG
import DimensionalData as DD
using Zarr
using DiscreteGlobalGridsVisualization: dggpoly, dggpoly!
using GLMakie, GeoMakie
GLMakie.activate!(inline = true)

A cube to write

The example uses two level-1 IGEO7 cells and their level-4 descendants as a small regional store. CellVector names the cells, CellLookup turns them into a one-level axis, and Cells makes that axis a cube dimension.

julia

sys = DGG.IGeo7System()
roots = DGG.CellVector(DGG.levelgrid(sys, 1))[1:2]
cells = sort!(reduce(vcat, [collect(DGG.descendants(sys, c, 4)) for c in roots]))
lookup = DGG.CellLookup(DGG.CellVector(sys, 4, cells))

CellLookup(IGeo7System, level=4, ncells=629, 1 windows)

The two layers share the cell axis. Their values encode axis positions, which makes the later selection checks easy to read.

julia

n = length(cells)
elevation = Float32.(1:n)
slope = Float32.(0.5 .* (1:n))
cube = DD.DimStack((; elevation, slope), (DGG.Cells(lookup),))

┌ 629-element DimStack ┐
├──────────────────────┴────────────────────────────────────── dims ┐
  ↓ Cells CellLookup(IGeo7System, level=4, ncells=629, 1 windows)
├─────────────────────────────────────────────────────────── layers ┤
  :elevation eltype: Float32 dims: Cells size: 629
  :slope     eltype: Float32 dims: Cells size: 629
└───────────────────────────────────────────────────────────────────┘

Write the cube and read it back

dggwrite returns the path that dggread opens. chunks = 128 gives this small cube five chunks; :auto would place the whole example in one chunk.

julia

path = DGG.dggwrite(joinpath(mktempdir(), "demo.zarr"), cube; chunks = 128)
store = DGG.dggread(path)

┌ 629-element DimStack ┐
├──────────────────────┴────────────────────────────────────────── dims ┐
  ↓ Cells ChunkedCellLookup(IGeo7System, level=4, ncells=629, ranges)
├─────────────────────────────────────────────────────────────── layers ┤
  :elevation eltype: Float32 dims: Cells size: 629
  :slope     eltype: Float32 dims: Cells size: 629
├───────────────────────────────────────────────────────────────────────┴ metadata ┐
  Dict{String, Any} with 5 entries:
  "source"      => "/tmp/jl_2JJiYl/demo.zarr"
  "encoding"    => "ranges"
  "description" => StoreDescription(igeo7/z7int, level 4, RangesEncoding, coord…
  "conventions" => ["zarr-conventions/dggs", "xdggs"]
  "attrs"       => Dict{String, Any}("dggs"=>Dict{String, Any}("coordinate"=>"c…
└──────────────────────────────────────────────────────────────────────────────────┘

julia

DD.metadata(store)["description"]

StoreDescription(igeo7/z7int, level 4, RangesEncoding, coordinate "cell_id_ranges", 2 variables)

dggread reconstructs the grid description from store attributes. The returned arrays stay lazy until a value is requested; collect makes an explicit in-memory comparison with the original cube:

julia

axis = DD.lookup(store[:elevation], DGG.Cells)

ChunkedCellLookup(IGeo7System, level=4, ncells=629, ranges)

julia

collect(axis) == cells, collect(parent(store[:elevation])) == elevation

(true, true)

Selecting a region out of a store

The axis is a ChunkedCellLookup. It supports the same selectors as a CellLookup and uses the chunk manifest to locate the required id data. Three selectors name a single cell:

At(cell) — the cell itself;
Contains(cell) — the same cell;
Contains((lon, lat)) — the cell holding a point.

julia

c = cells[300]
at = store[:elevation][DGG.Cells(DD.At(c))]
contains_cell = store[:elevation][DGG.Cells(DD.Contains(c))]
contains_point = store[:elevation][DGG.Cells(DD.Contains((67.5, 66.7)))]
at, contains_cell, contains_point

(300.0f0, 300.0f0, 300.0f0)

Covering(target) selects every stored cell reached by the coverage of target. Here target is the extent of a level-2 ancestor, so the selection demonstrates the small spill beyond a region's exact boundary.

julia

target = DGG.cell_extent(DGG.levelgrid(sys, 2), DGG.ancestor(sys, cells[300], 2))
region = store[:elevation][DGG.Cells(DGG.Covering(target))]

┌ 89-element DimArray{Float32, 1} elevation ┐
├───────────────────────────────────────────┴───────────────── dims ┐
  ↓ Cells CellLookup(IGeo7System, level=4, ncells=89, 20 windows)
├───────────────────────────────────────────────────────── metadata ┤
  Dict{String, Any} with 1 entry:
  "_ARRAY_DIMENSIONS" => Any["cell_ids"]
└───────────────────────────────────────────────────────────────────┘
 Z7Cell("001000")  287.0
 Z7Cell("001001")  288.0
 Z7Cell("001002")  289.0
 ⋮                 
 Z7Cell("001614")  592.0
 Z7Cell("001615")  593.0
 Z7Cell("001631")  603.0

The selected result has a compressed in-memory CellLookup; the source remains a ChunkedCellLookup.

Which chunks a selection touches

A ChunkManifest describes the chunk grid in cell-axis positions. It records each chunk's bounds and maps an axis position to its chunk.

julia

manifest = DGG.chunkmanifest(axis, 128)

ChunkManifest(5 chunks of 128, 629 cells)

julia

DGG.nchunks(manifest), length(manifest)

(5, 629)

julia

DGG.chunkbounds(manifest, 5)

513:629

The elevation values equal their positions, so selected values reveal which positions were fetched. chunkof maps each position to its chunk:

julia

selected = Int.(collect(region))
sort(unique(DGG.chunkof.(Ref(manifest), selected)))

3-element Vector{Int64}:
 3
 4
 5

The figure shows stored cells coloured by chunk, the target as a dashed box, and the selected cells outlined. A compact spherical region can span several file chunks.

julia

chunk = DGG.chunkof.(Ref(manifest), 1:n)
corners = [(target.X[1], target.Y[1]), (target.X[2], target.Y[1]),
           (target.X[2], target.Y[2]), (target.X[1], target.Y[2])]
box = [c1 .+ t .* (c2 .- c1) for (c1, c2) in zip(corners, circshift(corners, -1))
       for t in range(0, 1; length = 30)]

fig = Figure(size = (780, 470))
ax = GeoAxis(fig[1, 1]; dest = "+proj=laea +lon_0=44 +lat_0=64",
    limits = ((-16.0, 104.0), (45.0, 83.0)),
    xticks = 0:20:100, yticks = 50:10:80,
    title = "Stored cells by chunk; the cells Covering selects, outlined")
plt = dggpoly!(ax, cube[:elevation]; color = chunk,
    colormap = cgrad(:Set2, 5; categorical = true), colorrange = (0.5, 5.5))
lines!(ax, GeoMakie.coastlines(); color = ("#212529", 0.55), linewidth = 0.6)
dggpoly!(ax, region; color = :transparent, strokecolor = :black, strokewidth = 0.8)
lines!(ax, box; color = :black, linewidth = 2, linestyle = :dash)
Colorbar(fig[1, 2], plt; label = "chunk", ticks = 1:5)
fig

How to make that outline land in fewer chunks is the subject of Subzone layout.

Choosing how the cell ids are stored

An encoding describes how the store lays out cell ids. encoding = :auto chooses from the axis shape:

`encoding`	Stores	`:auto` picks it when
`:ranges`	`(n, 2)` inclusive `[start, stop]` id intervals	the axis is sorted, unique and one level
`:dense`	one id per cell	otherwise; also the interop choice for readers without interval support

The Python package xdggs reads the dense layout only. target = :xdggs chooses it and checks the rest of what xdggs needs, so xr.open_dataset(path, engine="zarr").pipe(xdggs.decode) opens the result; see Writing a store for xdggs.

A ranges axis opens without reading coordinate data: length, chunk boundaries and selectors use rank/select arithmetic over its intervals. This store uses RangesEncoding, with one row per interval:

julia

size(Zarr.zopen(path)["cell_id_ranges"], 2)

merge chooses what one interval may span:

`merge`	A run is	Rows	Read back correctly by
`:step` (default)	ids adjacent as integers	more (91 here)	any reader that counts ids, grid-aware or not
`:rank`	consecutive cells	fewest (1 here)	a rank-aware reader such as this package

This axis is a single run of consecutive cells, so under :rank it is one row:

julia

ranks = DGG.dggwrite(joinpath(mktempdir(), "ranks.zarr"), cube; chunks = 128,
                     merge = :rank)
size(Zarr.zopen(ranks)["cell_id_ranges"], 2)

Reading a store by URL

dggread also opens a public gs://, s3:// or https:// store in place. An s3:// URL additionally requires using AWSS3 to activate Zarr's S3 support. A selection then fetches the chunks it needs. For example:

julia

pori = DGG.dggread("https://storage.googleapis.com/geo-assets/igeo7-zarr/pori_z7_r10.zarr")

dggwrite writes to a local path or an open Zarr.ZGroup; publishing a remote store means uploading the directory it produced.

Out of core sweeps a kernel over a store chunk by chunk, starting from a store like the one written here.

This page was generated using Literate.jl.

A round trip through a DGGS store ​

A cube to write ​

Write the cube and read it back ​

Selecting a region out of a store ​

Which chunks a selection touches ​

Choosing how the cell ids are stored ​

Reading a store by URL ​

A round trip through a DGGS store

A cube to write

Write the cube and read it back

Selecting a region out of a store

Which chunks a selection touches

Choosing how the cell ids are stored

Reading a store by URL