Subsetting at the Protocol Level

A DAP2 client implemented over HTTPS, and the access-cost model it makes visible
Kahli Burke
Abstract

A DAP2 client — .dds, .das, and a reader for the XDR payload of .dods — run against NOAA's OISST v2 high-resolution analysis on a live Hyrax endpoint, and used to look at bytes per request, the latency floor and what concurrency does to it, and the difference between a subset touching one chunk and one touching five hundred. Later sections cover a DAP4 client returning the same values over the chunked .dap response, a DMR++ chunk map with a byte-range read against object storage, and a reduction in which two cloud workers read 133 MB of the record and return 506 KB.

A DAP2 client written against the spec, pointed at NOAA's OISST v2 quarter-degree daily analysis [1] on a public Hyrax endpoint. Later sections move to DAP4, to a DMR++ chunk map with byte-range reads against S3, and to a reduction split across two cloud workers.

The client

.dds for structure, .das for attributes, .dods?CE for the binary payload, .ascii?CE for the same values as text, which is convenient for checking a parser against.

OISST serves sst as Float32 in °C with no packing, so there is no scale_factor/add_offset handling below. Land and ice come back as the fill value, a large negative number.

import Downloads

# A DAP2 client over plain HTTPS. Every dataset exposes sibling responses on one base URL:
# .dds (structure), .das (attributes), .dods?CE (binary payload for a constraint expression).
const HYRAX = "https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.oisst.v2.highres/sst.day.mean.2023.nc"

# GET a URL as raw bytes, retrying transient faults (flaky HTTP/2 streams, timeouts) with a short
# backoff, so a single server fault does not surface as a failed cell.
function dget(url; tries = 4)
    local err
    for k in 1:tries
        try
            io = IOBuffer()
            Downloads.request(url; output = io, timeout = 120)
            return take!(io)
        catch e
            err = e
            k < tries && sleep(0.5 * k)
        end
    end
    throw(err)
end

"DAP client ready — base = OISST v2 high-res daily SST (NOAA PSL Hyrax)"
"DAP client ready — base = OISST v2 high-res daily SST (NOAA PSL Hyrax)"
# Parse a DAP2 .dods response for a Grid subset. After the "Data:\n" marker the payload is
# XDR big-endian: the data array (element count written twice, then the values), followed by
# each map vector the same way.
function parse_dods_grid(bytes, nt, ny, nx)
    marker = codeunits("\nData:\n")
    mi = findfirst(i -> @view(bytes[i:i+length(marker)-1]) == marker, 1:length(bytes)-length(marker))
    mi === nothing && error("no DAP Data: marker — malformed .dods response")
    io = IOBuffer(@view bytes[mi+length(marker):end])
    rdvec(T, n) = (read(io, UInt32); read(io, UInt32); [ntoh(read(io, T)) for _ in 1:n])
    flat = rdvec(Float32, nt * ny * nx)          # sst[time][lat][lon], C-order
    time = rdvec(Float64, nt)
    lat  = rdvec(Float32, ny)
    lon  = rdvec(Float32, nx)
    sst  = reshape(flat, nx, ny, nt)             # → [lon, lat, time] column-major
    (sst = sst, time = time, lat = lat, lon = lon)
end

struct Slab; sst; lat; lon; day; bytes::Int; ce::String; end

# A span of days in one request. Dimension-index ranges are 0-based and inclusive, per DAP.
# Retries GET+parse: a truncated body shows up as an EOF during parse.
function fetch_span(t0, t1, y0, y1, x0, x1)
    ce = "sst[$t0:$t1][$y0:$y1][$x0:$x1]"
    nt = t1 - t0 + 1
    for k in 1:3
        raw = dget(HYRAX * ".dods?" * ce)
        try
            g = parse_dods_grid(raw, nt, y1 - y0 + 1, x1 - x0 + 1)
            per = length(raw) ÷ nt
            slabs = [Slab(g.sst[:, :, n], g.lat, g.lon, Int(g.time[n]), per, ce) for n in 1:nt]
            return (slabs = slabs, bytes = length(raw), ce = ce)
        catch e
            k == 3 && rethrow()
            sleep(0.4k)
        end
    end
end

fetch_slab(t, y0, y1, x0, x1) = first(fetch_span(t, t, y0, y1, x0, x1).slabs)

# Check against the ASCII response for the same indices.
s = fetch_slab(0, 360, 362, 720, 723)
(ce = s.ce, bytes_on_wire = s.bytes, lat = s.lat, lon = s.lon,
 row0 = s.sst[:, 1], matches_ascii = isapprox(s.sst[1, 1], 27.51; atol = 0.01))
(ce = "sst[0:0][360:362][720:723]", bytes_on_wire = 374, lat = Float32[0.125, 0.375, 0.625], lon = Float32[180.125, 180.375, 180.625, 180.875], row0 = Float32[27.51, 27.48, 27.46, 27.439999], matches_ascii = true)

A single-day regional subset

A box over the western North Atlantic on 1 January 2023, given in degrees and converted to index ranges before the request goes out. The Gulf Stream front separates shelf water from the warm core; land drops out at the fill value.

# Address the grid in degrees instead of indices. OISST: lon 0.125→359.875, lat −89.875→89.875,
# both on a 0.25° step. A geographic box becomes an index range → a constraint expression.
loni(lon) = clamp(round(Int, (mod(lon, 360) - 0.125) / 0.25), 0, 1439)
lati(lat) = clamp(round(Int, (lat + 89.875) / 0.25), 0, 719)
fetch_region(lonW, lonE, latS, latN; day = 0) =
    fetch_slab(day, lati(latS), lati(latN), loni(lonW), loni(lonE))

# The Gulf Stream on New Year's Day: warm water peeling off Cape Hatteras into the cold N. Atlantic.
gulf = fetch_region(280, 315, 25, 48; day = 0)
@info "pulled $(gulf.ce) → $(gulf.bytes) bytes on the wire ($(size(gulf.sst,1))×$(size(gulf.sst,2)) grid), $(slab_date(gulf))"
staticmap(gulf; title = "Gulf Stream SST — $(slab_date(gulf))")
[ Info: pulled sst[0:0][460:552][1120:1260] → 53692 bytes on the wire (141×93 grid), 2023-01-01
output

Constraint expressions

The full variable is sst[time = 365][lat = 720][lon = 1440]. sst[0:0][360:362][720:723] selects one day, three latitudes and four longitudes, and the server returns those twelve values with the coordinate variables that label them.

Requested as ASCII below so the Grid structure is legible. The values match what the binary parser produced above.

# What a constraint expression actually is: name the array and the index ranges you want.
# Here — one day, a 3×4 patch of ocean near the dateline — asked for as human-readable ASCII.
ce = "sst[0:0][360:362][720:723]"
print(String(dget(HYRAX * ".ascii?" * ce)))
Dataset {
    Grid {
     ARRAY:
        Float32 sst[time = 1][lat = 3][lon = 4];
     MAPS:
        Float64 time[time = 1];
        Float32 lat[lat = 3];
        Float32 lon[lon = 4];
    } sst;
} Datasets/noaa.oisst.v2.highres/sst.day.mean.2023.nc;
---------------------------------------------
sst.sst[1][3][4]
[0][0], 27.51, 27.48, 27.46, 27.439999
[0][1], 27.49, 27.48, 27.47, 27.46
[0][2], 27.49, 27.48, 27.48, 27.47

sst.time[1]
81449.0

sst.lat[3]
0.125, 0.375, 0.625

sst.lon[4]
180.125, 180.375, 180.625, 180.875

Transfer cost

Sea-surface temperature around the Gulf Stream on one day, answered three ways against the same endpoint: retrieve the file, retrieve a global field for that day and crop locally, or constrain to the region server-side.

# Three ways to answer one question, measured in bytes off the wire.
full_year  = 365 * 720 * 1440 * 4         # the whole sst variable, as you'd get by downloading the file
one_day    = 720 * 1440 * 4               # a global map for the day, then cropped locally
box_only   = gulf.bytes                   # the constraint expression: just the Gulf Stream box, one day
vals  = [full_year, one_day, box_only]
ratio = round(Int, full_year / box_only)

human(b) = b ≥ 2^30 ? string(round(b / 2^30, digits = 2), " GB") :
           b ≥ 2^20 ? string(round(b / 2^20, digits = 1), " MB") :
                      string(round(b / 2^10, digits = 1), " KB")
cats = ["Download the file", "Fetch one global day", "Constraint expression"]

bar(i) = (value = vals[i],
          label = (show = true, position = "top", formatter = human(vals[i]),
                   fontSize = 16, fontWeight = "bold", color = "#f2b880"))

ch = echart(;
    # A valid font stack (not the CSS keyword 'inherit') so the canvas honours every fontSize.
    textStyle = (fontFamily = "'Helvetica Neue', Helvetica, Arial, sans-serif",),
    title = (text = "Bytes to answer one question",
             subtext = "$(ratio)× less data with a constraint expression than downloading the file",
             left = "center", textStyle = (color = "#e5e7eb", fontSize = 17, fontWeight = "bold"),
             subtextStyle = (color = "#93a1b0", fontSize = 12.5)),
    grid  = (top = 82, left = 66, right = 30, bottom = 46, containLabel = true),
    xAxis = (type = :category, data = cats, axisLabel = (color = "#c2ccd8", fontSize = 13)),
    yAxis = (type = :log, min = 1000, name = "bytes off the wire (log scale)",
             nameTextStyle = (color = "#93a1b0", fontSize = 12), axisLabel = (color = "#93a1b0", fontSize = 12),
             splitLine = (lineStyle = (color = "#2a2f3a",),)),
    series = [(type = :bar, barWidth = "46%",
               itemStyle = (color = "#f2a65a", borderRadius = [4, 4, 0, 0]),
               data = [bar(1), bar(2), bar(3)])],
    tooltip = (show = false,))
ch

A month of daily subsets

Thirty consecutive days of the same box, one constraint expression each. The front meanders, pinches off warm-core rings, and advects them downstream.

# The time axis is contiguous, so the month is one hyperslab rather than thirty daily requests.
jan = fetch_span(0, 29, lati(25), lati(48), loni(280), loni(315))
month_slabs = jan.slabs
@info "$(jan.ce) → $(round(jan.bytes / 1024, digits = 1)) KB, $(length(month_slabs)) daily fields"
mapview(month_slabs; fps = 6, title = "Gulf Stream SST — January 2023 (daily)")
[ Info: sst[0:29][460:552][1120:1260] → 1538.1 KB, 30 daily fields

Arbitrary regions

Selecting a region and a day rewrites the constraint expression and redraws from the response.

The boxes all lie within 0–360°E so each maps to a contiguous longitude index range. One spanning the prime meridian would need two requests or a wrapped read.

[ Info: Gulf Stream (N. Atlantic) · day 0 → sst[0:0][460:552][1120:1260] → 53692 bytes off the wire
output

Other collections

The client knows nothing about the dataset beyond the variable name, so pointing it at other collections needs only the URL. Four PSL collections queried for their structure below, including a 4-D pressure-level field.

datasetvariableshapefull_GB
Sea-surface temperature — OISST daily 0.25°sst(365, 720, 1440)1.41 GB
Air temperature on 17 pressure levels — NCEPair(365, 17, 73, 144)0.24 GB
Outgoing longwave radiation — daily, 1974→olr(17746, 73, 144)0.69 GB
Land surface max temperature — CPC 0.5°tmax(365, 360, 720)0.35 GB

A reduction over the full record

The Oceanic Niño Index is the area-averaged SST of the Niño 3.4 box (5°S–5°N, 170°W–120°W) as an anomaly against the seasonal cycle. The constraint expression follows from the definition: that box, every month of the record, about seventeen megabytes of a two-gigabyte variable. Area weighting, climatology and the running mean are computed locally.

The request goes out as concurrent time bands rather than one response. A multi-decade series through a single box cuts across the file's chunking — see below — so it gains from overlapping the round trips, and the individual responses stay small enough to stream cleanly.

Record length comes from the .dds rather than a literal, so this stays correct as months are appended.

# Niño 3.4 (5°S–5°N, 170°W–120°W) is the ocean's ENSO thermostat. Pull ONLY that box across the
# entire monthly record — a sliver of a multi-gigabyte variable — and reduce it to one per month.
NINO_BASE = "https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.oisst.v2.highres/sst.mon.mean.nc"

# Read the record length from the .dds instead of hard-coding it, so this stays correct as NOAA
# keeps appending months to the file.
function dds_len(base, dim)
    dds = String(dget(base * ".dds"))
    parse(Int, match(Regex("\\b$dim\\[$dim = (\\d+)\\]"), dds).captures[1])
end
NMONTH = dds_len(NINO_BASE, "time")

# Ask for the record in concurrent TIME BANDS rather than as one long response. A multi-decade series
# through a single box is the access pattern that cuts most sharply across the file's chunking, so it
# is also the one that gains most from overlapping the round-trips — and it keeps each response small
# enough to stream cleanly. Same total bytes, same answer, a fraction of the wall clock.
function fetch_band(base, t0, t1, y0, y1, x0, x1)
    ce = "sst[$t0:$t1][$y0:$y1][$x0:$x1]"
    for k in 1:5
        raw = dget(base * ".dods?" * ce)
        try
            return (grid = parse_dods_grid(raw, t1 - t0 + 1, y1 - y0 + 1, x1 - x0 + 1), bytes = length(raw))
        catch e
            k == 5 && rethrow()                          # short body ⇒ EOF mid-parse ⇒ re-ask this band
            sleep(0.5k)
        end
    end
end

function fetch_box(base, nt, y0, y1, x0, x1; band = 60)
    spans = [(t, min(t + band - 1, nt - 1)) for t in 0:band:nt-1]
    parts = asyncmap(s -> fetch_band(base, s[1], s[2], y0, y1, x0, x1), spans; ntasks = length(spans))
    (sst   = cat((p.grid.sst for p in parts)...; dims = 3),
     time  = vcat((p.grid.time for p in parts)...),
     lat   = parts[1].grid.lat, lon = parts[1].grid.lon,
     bytes = sum(p.bytes for p in parts), requests = length(parts))
end

x0 = loni(190); x1 = loni(240); y0 = lati(-5); y1 = lati(5)
box = fetch_box(NINO_BASE, NMONTH, y0, y1, x0, x1)

# area-weighted (cos φ) spatial mean over the box → one SST per month
wts = cosd.(box.lat)
function monthly_mean(sst, w)
    nt = size(sst, 3); out = zeros(Float64, nt)
    for t in 1:nt
        s = 0.0; ws = 0.0
        for j in axes(sst, 2), i in axes(sst, 1)
            v = sst[i, j, t]
            v > -100f0 && (s += w[j] * v; ws += w[j])
        end
        out[t] = s / ws
    end
    out
end
sst_series = monthly_mean(box.sst, wts)

# calendar-month climatology → anomaly → 3-month running mean (a simplified ONI)
enso_dates = [Date(1800, 1, 1) + Day(round(Int, t)) for t in box.time]
clim = [mean(sst_series[month.(enso_dates) .== m]) for m in 1:12]
anom = [sst_series[k] - clim[month(enso_dates[k])] for k in eachindex(sst_series)]
oni  = [mean(anom[max(1, k-1):min(length(anom), k+1)]) for k in eachindex(anom)]

(months = length(sst_series), requests = box.requests,
 MB_on_wire = round(box.bytes / 2^20, digits = 1),
 full_variable_GB = round(NMONTH * 720 * 1440 * 4 / 2^30, digits = 1),
 span = (first(enso_dates), last(enso_dates)), latest = round(oni[end], digits = 2),
 biggest_el_nino = enso_dates[argmax(oni)] => round(maximum(oni), digits = 2))
(months = 538, requests = 9, MB_on_wire = 16.9, full_variable_GB = 2.1, span = (Dates.Date("1981-09-01"), Dates.Date("2026-06-01")), latest = 1.32, biggest_el_nino = Dates.Date("2015-12-01") => 2.66)
xs = [string(year(d), "-", lpad(month(d), 2, '0')) for d in enso_dates]   # "1997-12" — unique per month
ys = round.(oni, digits = 3)

# Everything the chart says about itself is derived from what actually came back, so the caption
# cannot drift out of step with the record as NOAA appends months.
span_txt = "$(year(first(enso_dates)))–$(year(last(enso_dates)))"
cost_txt = "reduced from $(round(box.bytes / 2^20, digits = 1)) MB of a " *
           "$(round(NMONTH * 720 * 1440 * 4 / 2^30, digits = 1)) GB record"

echart(:line, xs, ys;
    textStyle = (fontFamily = "'Helvetica Neue', Helvetica, Arial, sans-serif", color = "#c2ccd8"),
    title = (text = "Niño 3.4 SST anomaly — the ENSO signal, $(span_txt)",
             subtext = "area-weighted °C, 3-month running mean · $(cost_txt)",
             left = "center", textStyle = (color = "#e5e7eb", fontSize = 16),
             subtextStyle = (color = "#93a1b0", fontSize = 12)),
    grid = (top = 74, left = 36, right = 20, bottom = 30, containLabel = true),
    xAxis = (type = :category, data = xs, boundaryGap = false,
             axisLabel = (color = "#93a1b0", interval = 59, formatter = "{value}"),
             axisTick = (show = false,), axisLine = (lineStyle = (color = "#3a4150",),)),
    yAxis = (name = "°C", nameTextStyle = (color = "#93a1b0",), axisLabel = (color = "#93a1b0",),
             splitLine = (lineStyle = (color = "#232831",),)),
    showSymbol = false, lineStyle = (width = 1.6, color = "#f2a65a"),
    areaStyle = (color = "#f2a65a", opacity = 0.12),
    markLine = (silent = true, symbol = "none", label = (show = false,),
                lineStyle = (color = "#5a6472", type = "dashed"),
                data = [(yAxis = 0.5,), (yAxis = -0.5,)]),
    tooltip = (trigger = "axis",))

Latency and concurrency

A single modest subset is dominated by latency — round trip plus the server's time to evaluate the constraint and read the underlying bytes — and that component does not shrink by asking for less.

Sixteen 10°×10° tiles, one constraint expression each, fetched sequentially and then concurrently. Sequential pays the latency once per tile; concurrent overlaps them and approaches the slowest single request. The gap widens with tile count.

This is how a Dask or Zarr client reads a chunked array from object storage.

# A basin covered by a grid of 10°×10° tiles — each tile is one constraint expression / one request.
tiles = [(lonW, lonW + 10, latS, latS + 10) for lonW in 300:10:330 for latS in 20:10:50]

# Sequential: fetch the tiles one after another.
t_seq = @elapsed seqres = [fetch_region(t...; day = 0) for t in tiles]
# Concurrent: the same requests in flight together (async I/O — the Dask/Zarr way to read chunks).
t_par = @elapsed parres = asyncmap(t -> fetch_region(t...; day = 0), tiles; ntasks = length(tiles))

(tiles = length(tiles), total_MB = round(sum(s.bytes for s in seqres) / 2^20, digits = 2),
 latency_ms_per_request = round(Int, 1000 * t_seq / length(tiles)),
 sequential_s = round(t_seq, digits = 1), concurrent_s = round(t_par, digits = 2),
 speedup = round(t_seq / t_par, digits = 1))
(tiles = 16, total_MB = 0.11, latency_ms_per_request = 369, sequential_s = 5.9, concurrent_s = 0.72, speedup = 8.2)

Chunking and access patterns

OISST stores sst with a chunk shape of (1, 720, 1440) — one global field per timestep — which the .das reports and the cell below reads.

The two access patterns used above land very differently under that layout. A spatial box on one day touches a single chunk: 4 MB read to return a few kilobytes. The Niño 3.4 series touches one chunk per month, over five hundred of them, each read whole — on the order of two gigabytes read to return seventeen megabytes.

Cloud-optimized access attacks this from both ends: shape requests to touch fewer chunks, or take the chunk map and read the required byte ranges from object storage directly. The DMR++ section does the second.

# The .das carries the on-disk HDF5 chunk shape.
function chunk_shape(base; var = "sst")
    das = String(dget(base * ".das"))
    blk = match(Regex(var * raw"\s*\{.*?\}", "s"), das)          # the sst attribute block
    blk === nothing && return "?"
    m = match(r"_ChunkSizes\s+([\d,\s]+?);", blk.match)
    m === nothing ? "contiguous (unchunked)" : replace(strip(m.captures[1]), r"\s+" => "")
end
daily_chunk   = chunk_shape(HYRAX)        # sst.day.mean.2023.nc  → time,lat,lon
monthly_chunk = chunk_shape(NINO_BASE)    # sst.mon.mean.nc

# Chunks touched by each access pattern, and the bytes the server must read as a result.
gulf_chunks_touched = 1                   # one day, any spatial box → a single (1,720,1440) chunk
nino_chunks_touched = NMONTH              # one small box, every month → one chunk per month

chunk_read_GB(n) = round(n * 720 * 1440 * 4 / 2^30, digits = 2)
(daily_chunk = daily_chunk, monthly_chunk = monthly_chunk,
 gulf_map_touches_chunks = gulf_chunks_touched,
 gulf_server_reads_MB = round(720 * 1440 * 4 / 2^20, digits = 1),
 nino_series_touches_chunks = nino_chunks_touched,
 nino_server_reads_GB = chunk_read_GB(nino_chunks_touched),
 nino_sends_MB = round(box.bytes / 2^20, digits = 1))
(daily_chunk = "1,720,1440", monthly_chunk = "1,720,1440", gulf_map_touches_chunks = 1, gulf_server_reads_MB = 4.0, nino_series_touches_chunks = 538, nino_server_reads_GB = 2.08, nino_sends_MB = 16.9)

DAP4

DAP4 replaces the .dds/.das split with the DMR — one typed XML document with real types, groups and shared dimensions, describing the constrained response rather than the whole dataset. Data comes back over .dap as a chunked binary stream, each chunk headed by a flag byte and a 24-bit length, the first carrying the DMR and the rest the values, with a CRC32 per variable at the end. Constraint expressions are fully qualified and carried in dap4.ce.

The endpoint serves both protocols, so the reader below asks for the same twelve values as the DAP2 example and checks them against it.

Two details in the reader. Byte order is a property of the response, not of each chunk, and servers vary in which chunk carries the flag, so it is honoured if any chunk sets it. And the trailing checksum has to be kept out of the value count.

(constraint = "/sst[0][360:362][720:723]", chunks_in_response = 2, dmr_bytes = 4281, byte_order = "little-endian", checksum_crc32 = "5c6f5a29", values_degC = Float32[27.51, 27.48, 27.46, 27.44, 27.49, 27.48, 27.47, 27.46, 27.49, 27.48, 27.48, 27.47], same_numbers_as_dap2 = true)

DMR++ and byte-range access

A DMR++ records, per variable, the byte ranges holding each chunk and the filters needed to decode them, so a client can work out which chunks its subset needs and pull them from object storage directly.

Below, a real DMR++ for a MERRA-2 granule from OPeNDAP's test corpus, its chunk map extracted — 97 chunks as (offset, nBytes) — followed by the range read: one request, 4608 bytes out of a 472 MB object.

The DMR++ is public; the granule it points at is behind Earthdata Login, so the range read runs against a public object instead. The difference is access control, which in the deployed system is most of the architecture.

Moving the computation

Earthdata is not one bucket in one region, and once the client is small the analysis can sit next to whichever archive it reads.

The two cells below run on remote workers on different cloud providers, rega and regc. Each takes one hemisphere and reads the whole 45-year monthly record for its half at 1° — eighteen constraint expressions, tens of megabytes, none of it leaving that machine. Per ocean cell it removes the seasonal cycle, fits a linear trend, and correlates the anomaly against the Niño 3.4 index computed earlier.

The index goes up, about four kilobytes; the record it is regressed against stays put. Two small arrays per hemisphere come back and are joined below.

Top panel: the ENSO teleconnection, with the equatorial Pacific tongue, out-of-phase midlatitude lobes, and the Indian Ocean basin-wide response. Bottom: the linear trend over the same record, positive nearly everywhere, strongest in the Arctic and the western boundary currents, with the subpolar North Atlantic minimum and the Southern Ocean as exceptions.

import Downloads
using Statistics

# The shared kernel. Slate re-establishes a cell's imports and DEFINITIONS on every worker, so the
# region cells below call these functions locally and only their results cross the boundary — the
# functions themselves never travel as data (they live in the notebook's own module, which the far
# side can't decode).
#
# For a longitude band: read the whole 45-year monthly record at 1° and, per ocean cell, remove the
# seasonal cycle, fit a linear trend, and correlate the anomaly with an index supplied by the caller.
const MON = "https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.oisst.v2.highres/sst.mon.mean.nc"

function mesh_get(url; tries = 5)
    local err
    for k in 1:tries
        try
            io = IOBuffer(); Downloads.request(url; output = io, timeout = 180)
            return take!(io)
        catch e
            err = e; k < tries && sleep(0.6k)
        end
    end
    throw(err)
end

function mesh_parse(bytes, nt, ny, nx)
    mk = codeunits("\nData:\n")
    mi = findfirst(i -> @view(bytes[i:i+length(mk)-1]) == mk, 1:length(bytes)-length(mk))
    mi === nothing && error("no DAP Data: marker")
    io = IOBuffer(@view bytes[mi+length(mk):end])
    rd(T, n) = (read(io, UInt32); read(io, UInt32); [ntoh(read(io, T)) for _ in 1:n])
    flat = rd(Float32, nt * ny * nx); rd(Float64, nt)
    lat = rd(Float32, ny); lon = rd(Float32, nx)
    (sst = reshape(flat, nx, ny, nt), lat = lat, lon = lon)
end

# Banded and concurrent, as in the local client. The strided CE thins 0.25° to 1° server-side, so
# the decimation costs nothing on the wire.
function mesh_fetch(x0, x1, nt; s = 4, band = 30, ntasks = 8)
    spans = [(t, min(t + band - 1, nt - 1)) for t in 0:band:nt-1]
    ny = length(0:s:719); nx = length(x0:s:x1)
    parts = asyncmap(spans; ntasks = ntasks) do sp
        a, b = sp
        ce = "sst[$a:$b][0:$s:719][$x0:$s:$x1]"
        for k in 1:5
            raw = mesh_get(MON * ".dods?" * ce)
            try
                return (grid = mesh_parse(raw, b - a + 1, ny, nx), bytes = length(raw))
            catch e
                k == 5 && rethrow(); sleep(0.5k)
            end
        end
    end
    (sst = cat((p.grid.sst for p in parts)...; dims = 3),
     lat = parts[1].grid.lat, lon = parts[1].grid.lon,
     bytes = sum(p.bytes for p in parts), requests = length(parts))
end

# `index` is computed by the caller and sent up; only the two maps that get drawn come back.
function ocean_response(x0, x1, index; s = 4, band = 30)
    nt = length(index)
    t_fetch = @elapsed g = mesh_fetch(x0, x1, nt; s = s, band = band)
    nx, ny = size(g.sst, 1), size(g.sst, 2)

    ix = Float64.(index) .- mean(index)            # centred index
    den_i = sum(abs2, ix)
    tt = collect(0.0:nt-1); tt .-= mean(tt)        # centred time axis, in months
    den_t = sum(abs2, tt)

    trend = fill(NaN32, nx, ny)                    # °C per decade
    corr  = fill(NaN32, nx, ny)
    a = zeros(Float64, nt)

    t_comp = @elapsed for j in 1:ny, i in 1:nx
        ok = true
        @inbounds for m in 1:nt
            g.sst[i, j, m] < -100f0 && (ok = false; break)      # land / ice fill value
        end
        ok || continue
        @inbounds for m in 1:nt; a[m] = g.sst[i, j, m]; end
        # Deseasonalise by subtracting the mean of every 12th sample. Any consistent 12-phase
        # grouping works; it need not start in January.
        for p in 1:12
            idx = p:12:nt
            μ = mean(@view a[idx])
            @inbounds for m in idx; a[m] -= μ; end
        end
        sa = sum(abs2, a)
        sa == 0 && continue
        trend[i, j] = Float32(sum(tt[m] * a[m] for m in 1:nt) / den_t * 120)   # per month → per decade
        corr[i, j]  = Float32(sum(ix[m] * a[m] for m in 1:nt) / sqrt(sa * den_i))
    end

    (trend = trend, corr = corr, lat = g.lat, lon = g.lon,
     bytes_in = g.bytes, requests = g.requests, months = nt,
     secs_fetch = round(t_fetch, digits = 1), secs_compute = round(t_comp, digits = 1))
end
"mesh kernel ready — fetch + deseasonalize + trend + correlation, primed on every worker"
"mesh kernel ready — fetch + deseasonalize + trend + correlation, primed on every worker"
# rega's half. `oni` is computed on the laptop and is the only thing sent up; the 45-year record for
# this hemisphere is read and reduced on this machine. Names are region-private so nothing but
# `west` has to cross the boundary.
tw = @elapsed rw = ocean_response(0, 719, oni)               # Western hemisphere: 0–180°E
west = (cloud = "rega", host = gethostname(), os = string(Sys.KERNEL, " ", Sys.MACHINE),
        hemisphere = "Western (0–180°E)",
        MB_in = round(rw.bytes_in / 2^20, digits = 1),
        KB_out = round(2 * length(rw.trend) * 4 / 1024, digits = 1),
        index_KB_up = round(length(oni) * 8 / 1024, digits = 1),
        requests = rw.requests, secs_fetch = rw.secs_fetch, secs_compute = rw.secs_compute,
        secs = round(tw, digits = 1),
        trend = rw.trend, corr = rw.corr, lat = rw.lat, lon = rw.lon)
"rega · $(west.os) · $(west.hemisphere) · $(west.index_KB_up) KB index up → $(west.MB_in) MB read here → " *
"$(west.KB_out) KB back · $(west.requests) requests · fetch $(west.secs_fetch)s + compute $(west.secs_compute)s"
"rega · Linux x86_64-linux-gnu · Western (0–180°E) · 4.2 KB index up → 66.5 MB read here → 253.1 KB back · 18 requests · fetch 4.5s + compute 3.4s"
# regc's half, on a different cloud provider, running the same kernel against the same index.
te = @elapsed re = ocean_response(720, 1439, oni)            # Eastern hemisphere: 180–360°E
east = (cloud = "regc", host = gethostname(), os = string(Sys.KERNEL, " ", Sys.MACHINE),
        hemisphere = "Eastern (180–360°E)",
        MB_in = round(re.bytes_in / 2^20, digits = 1),
        KB_out = round(2 * length(re.trend) * 4 / 1024, digits = 1),
        index_KB_up = round(length(oni) * 8 / 1024, digits = 1),
        requests = re.requests, secs_fetch = re.secs_fetch, secs_compute = re.secs_compute,
        secs = round(te, digits = 1),
        trend = re.trend, corr = re.corr, lat = re.lat, lon = re.lon)
"regc · $(east.os) · $(east.hemisphere) · $(east.index_KB_up) KB index up → $(east.MB_in) MB read here → " *
"$(east.KB_out) KB back · $(east.requests) requests · fetch $(east.secs_fetch)s + compute $(east.secs_compute)s"
"regc · Linux x86_64-linux-gnu · Eastern (180–360°E) · 4.2 KB index up → 66.5 MB read here → 253.1 KB back · 18 requests · fetch 4.8s + compute 4.4s"
# Assembled locally. `west` and `east` are read at top level so the DAG sees the dependency on both
# region cells. Each half is [lon, lat] and they differ in longitude, so join along lon.
gcorr  = vcat(west.corr,  east.corr)         # [lon = 360, lat = 180]
gtrend = vcat(west.trend, east.trend)
glon   = vcat(west.lon,   east.lon)          # 0.5 … 359.5°E
glat   = west.lat                            # −89.5 … 89.5°N, shared

MB_in  = west.MB_in + east.MB_in
KB_out = west.KB_out + east.KB_out
years  = round(Int, NMONTH / 12)

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

ax1 = Axis(fig[1, 1]; ylabel = "latitude (°N)", limits = (0.0, 360.0, -90.0, 90.0),
           title = "Correlation of SST anomaly with the Niño 3.4 index")
hm1 = heatmap!(ax1, glon, glat, gcorr; colormap = :balance, colorrange = (-0.8, 0.8))
coast!(ax1; color = (:white, 0.35), linewidth = 0.5)
Colorbar(fig[1, 2], hm1; label = "correlation")
ax1.aspect = DataAspect(); hidexdecorations!(ax1; grid = false)

ax2 = Axis(fig[2, 1]; xlabel = "longitude (°E)", ylabel = "latitude (°N)", limits = (0.0, 360.0, -90.0, 90.0),
           title = "Linear SST trend over the same $(years)-year record")
hm2 = heatmap!(ax2, glon, glat, gtrend; colormap = :vik, colorrange = (-0.4, 0.4))
coast!(ax2; color = (:white, 0.35), linewidth = 0.5)
Colorbar(fig[2, 2], hm2; label = "°C / decade")
ax2.aspect = DataAspect()

Label(fig[0, :], "Computed on rega + regc from $(MB_in) MB read on those machines\n" *
      "$(west.index_KB_up) KB of index up, $(KB_out) KB of results back",
      fontsize = 13, color = (:white, 0.72), padding = (0, 0, 4, 0))

@info "assembled on $(gethostname()) ($(Sys.KERNEL)) · rega [$(west.host), $(west.MB_in) MB in, $(west.secs)s] " *
      "+ regc [$(east.host), $(east.MB_in) MB in, $(east.secs)s] · " *
      "$(round(Int, MB_in * 1024 / KB_out))× reduction"
figpng(fig)
[ Info: assembled on heinlein.local (Darwin) · rega [localhost, 66.5 MB in, 8.7s] + regc [main, 66.5 MB in, 10.1s] · 269× reduction
output

Summary

A page of code covers enough of DAP2 to do real work against a production endpoint, and having the requests in hand made the cost of each step measurable while the analysis was assembled.

A 45-year climate index came out of seventeen megabytes of a two-gigabyte variable. Latency, not bandwidth, set the floor for a single subset, and concurrency collapsed it across sixteen tiles. Chunk layout accounted for the difference between a subset touching one chunk and one touching over five hundred. DAP4 returned the same values over a more clearly specified wire format. A DMR++ chunk map and one range request retrieved 4608 bytes from a 472 MB object with no server evaluating a constraint. Two cloud workers read 133 MB and returned 506 KB, having been sent four kilobytes of index.

Dataset: NOAA OISST [1]; format: netCDF [2]; protocol and server: DAP4 and Hyrax [3][4].

A companion notebook, Reading DAP Through the netCDF Library, works the same ground from the other side: a four-dimensional reanalysis through NCDatasets.jl with no protocol-level code, and a look at the requests the library issues.

References

  1. Huang, Boyin and Liu, Chunying and Banzon, Viva and Freeman, Eric and Graham, Garrett and Hankins, Bill and Smith, Tom and Zhang, Huai-Min. Improvements of the Daily Optimum Interpolation Sea Surface Temperature (DOISST) Version 2.1. 2021.
  2. Rew, Russ and Davis, Glenn. NetCDF: An Interface for Scientific Data Access. 1990.
  3. OPeNDAP, Inc.. DAP4 Data Access Protocol, Specification Volume 1. 2018.
  4. Gallagher, James and Potter, Nathan and Sgouros, Tom. Hyrax: An OPeNDAP Data Server. 2007.