Reading DAP Through the netCDF Library
A companion to Subsetting at the Protocol Level, which implemented a DAP2 client directly. Here the access goes through NCDatasets.jl and the netCDF-C library beneath it, delivered as a JLL. The dataset is the NCEP/NCAR Reanalysis on pressure levels — (lon, lat, level, time) — used to diagnose the January 2021 sudden stratospheric warming: the WMO wind-reversal criterion, 10 hPa temperature over the polar cap, and the downward propagation that follows. The diagnostic costs about 1.4 MB of a 249 MB variable, with no constraint expression written. The second half recovers the requests the library issues. Over DAP2 it translates array indexing into a hyperslab; over DAP4, from identical calling code, it asks for the variable without index ranges and subsets locally — 48 bytes against 14.6 MB for the same twelve values, on both a THREDDS and a Hyrax endpoint.
Most DAP traffic is not generated by DAP clients. It comes from netCDF, which has carried DAP support since the late 1990s and sits under nearly every analysis tool — xarray via netcdf4-python, ncdump, NCDatasets.jl here. Hand a URL to an API expecting a path and the dataset presents as an array.
The first half puts that to work: a four-dimensional field, an aggregation across files, CF conventions applied on the way through, and a stratospheric diagnostic computed without protocol-level code. The second half turns on libcurl logging and compares the request streams for DAP2 and DAP4 against the same data.
import NetworkOptions
# netCDF reads its runtime config from an rc file. Two settings here: Julia's CA bundle for TLS, and
# libcurl verbose logging, which a later cell reads back. Must be set before netCDF initialises.
const NCRC = joinpath(tempdir(), "slate_ncrc")
write(NCRC, """
HTTP.VERBOSE=1
HTTP.SSL.CAINFO=$(something(NetworkOptions.ca_roots(), NetworkOptions.ca_roots_path()))
""")
ENV["NCRCENV_RC"] = NCRC
using NCDatasets
# NCDatasets is a typed Julia layer over netCDF-C, which arrives as a JLL — the same C library
# `ncdump` and xarray use, and where DAP lives. The zeroed build date is a reproducible build.
const NETCDF_VERSION = first(split(unsafe_string(
ccall((:nc_inq_libvers, NCDatasets.NetCDF_jll.libnetcdf), Cstring, ())), " of "))
(NCDatasets = string(pkgversion(NCDatasets)),
netcdf_c = NETCDF_VERSION,
artifact = basename(NCDatasets.NetCDF_jll.libnetcdf_path))(NCDatasets = "0.14.15", netcdf_c = "4.10.1", artifact = "libnetcdf.22.dylib")
Opening a remote dataset
NCDataset takes the OPeNDAP URL where a filename would go, and the result behaves as usual: named dimensions, shapes, element types, attributes.
The dataset is the NCEP/NCAR Reanalysis on 17 standard pressure levels, daily for 2021 [1]. The vertical coordinate is what the diagnostics below need.
_FillValue comes back as Union{Missing,Float32} and time as DateTime, both applied by the library.
const PSL = "https://psl.noaa.gov/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/"
# A URL, where the API expects a filename. Nothing else about the code changes.
uwnd_ds = NCDataset(PSL * "pressure/uwnd.2021.nc")
air_ds = NCDataset(PSL * "pressure/air.2021.nc")
u = uwnd_ds["uwnd"] # zonal wind [lon, lat, level, time]
T = air_ds["air"] # temperature [lon, lat, level, time]
# The library hands back real Julia types. `Union{Missing,Float32}` is CF's _FillValue already
# decoded; `time` is already a DateTime with the calendar applied. None of that was our work.
levels = Array(uwnd_ds["level"])
lats = Array(uwnd_ds["lat"])
lons = Array(uwnd_ds["lon"])
times = Array(uwnd_ds["time"])
(variable = u.attrib["long_name"], units = u.attrib["units"],
dims = NCDatasets.dimnames(u), size = size(u),
eltype = string(eltype(u)),
full_variable_MB = round(prod(size(u)) * 4 / 2^20, digits = 1),
levels_hPa = levels, lat_range = (first(lats), last(lats)),
time_type = string(eltype(times)), first_day = first(times), last_day = last(times))(variable = "mean Daily U-wind", units = "m/s", dims = ("lon", "lat", "level", "time"), size = (144, 73, 17, 365), eltype = "Union{Missing, Float32}", full_variable_MB = 248.8, levels_hPa = Float32[1000.0, 925.0, 850.0, 700.0, 600.0, 500.0, 400.0, 300.0, 250.0, 200.0, 150.0, 100.0, 70.0, 50.0, 30.0, 20.0, 10.0], lat_range = (90.0f0, -90.0f0), time_type = "Dates.DateTime", first_day = Dates.DateTime("2021-01-01T00:00:00"), last_day = Dates.DateTime("2021-12-31T00:00:00"))The January 2021 warming
A major sudden stratospheric warming is defined by the zonal-mean zonal wind at 10 hPa, 60°N reversing from westerly to easterly [2]. The 2020–21 winter has one.
The read it requires is a single latitude and every longitude, across a winter spanning two yearly files. NCDatasets takes a vector of URLs with aggdim, giving one dataset with a 731-day time coordinate; the zonal mean is taken locally.
Keeping all 17 levels rather than just 10 hPa costs little and supplies the time–height section below from the same request.
using Statistics, Dates
# Two remote files presented as one dataset, aggregated along time by the library.
winter_ds = NCDataset([PSL * "pressure/uwnd.2020.nc", PSL * "pressure/uwnd.2021.nc"];
aggdim = "time", deferopen = false)
wt = Array(winter_ds["time"])
LAT60 = findmin(abs.(lats .- 60f0))[2] # 60°N
HPA10 = findmin(abs.(levels .- 10f0))[2] # 10 hPa
winter = findall(d -> Date(2020, 11, 1) ≤ Date(d) ≤ Date(2021, 3, 31), wt)
# A single indexing expression: every longitude, the 60°N row, all 17 levels, the winter's days.
# Retaining all levels also supplies the time–height section below from this one request.
block = winter_ds["uwnd"][:, LAT60, :, winter] # [lon, level, time]
# Zonal mean around the latitude circle → [level, time]
ubar = [mean(skipmissing(view(block, :, k, n))) for k in axes(block, 2), n in axes(block, 3)]
wdates = Date.(wt[winter])
u10 = ubar[HPA10, :] # the WMO diagnostic
# Major warming: the first day the 10 hPa / 60°N zonal-mean wind turns easterly.
reversed = findall(<(0), u10)
central = isempty(reversed) ? nothing : wdates[first(reversed)]
(read_values = length(block), read_KB = round(length(block) * 4 / 1024, digits = 1),
lat_used = lats[LAT60], level_used = levels[HPA10],
winter_days = length(winter),
min_wind = round(minimum(u10), digits = 1),
easterly_days = length(reversed),
ssw_central_date = central)(read_values = 369648, read_KB = 1443.9, lat_used = 60.0f0, level_used = 10.0f0, winter_days = 151, min_wind = -7.0f0, easterly_days = 12, ssw_central_date = Dates.Date("2021-01-05"))10 hPa temperature over the polar cap
10 hPa temperature north of 20°N on two dates: 22 December 2020 with the vortex intact, and 5 January 2021, by which the cold air has separated into two lobes with a warm anomaly over the pole — a split rather than a displacement.
LATMIN = 20
G = polar_grid(lons, lats, LATMIN)
CX, CY = polar_coast(LATMIN)
# Two dates: a fortnight before the reversal, and the reversal itself.
dates_shown = [Date(2020, 12, 22), Date(2021, 1, 5)]
# One field per date: [lon, lat] at 10 hPa. `air_ds` covers 2021 only, so use an aggregated
# temperature dataset spanning both years — the same trick as the wind.
temp_ds = NCDataset([PSL * "pressure/air.2020.nc", PSL * "pressure/air.2021.nc"];
aggdim = "time", deferopen = false)
tt = Date.(Array(temp_ds["time"]))
fields = map(dates_shown) do d
n = findfirst(==(d), tt)
f = temp_ds["air"][:, G.js, HPA10, n] # [lon, lat-cap] — a single day at 10 hPa
vcat(f, f[1:1, :]) # close the seam
end
bytes_per_field = length(fields[1]) * 4
lo, hi = extrema(vcat((collect(skipmissing(f)) for f in fields)...))
CMAP = Reverse(:RdBu) # warm = red, cold = blue
fig = Figure(size = (980, 520))
for (i, (d, f)) in enumerate(zip(dates_shown, fields))
ax = polar_axis(fig, (1, i), LATMIN; title = string(d))
surface!(ax, G.X, G.Y, zeros(size(G.X)); color = Float32.(coalesce.(f, NaN32)),
colormap = CMAP, colorrange = (lo, hi), shading = NoShading)
lines!(ax, CX, CY; color = (:white, 0.5), linewidth = 0.6)
end
Colorbar(fig[1, 3]; colormap = CMAP, limits = (lo, hi), label = "10 hPa temperature (K)")
Label(fig[0, :], "The polar vortex before and during the January 2021 sudden stratospheric warming",
fontsize = 15, color = (:white, 0.85))
@info "2 maps × $(round(bytes_per_field/1024, digits=1)) KB — out of a $(round(prod(size(T))*4/2^20, digits=0)) MB variable"
figpng(fig)[ Info: 2 maps × 16.4 KB — out of a 249.0 MB variable
Downward propagation
What makes these events useful at the two-to-six week range is that the circulation anomaly descends, arriving in the troposphere as a weakened jet, a negative Arctic Oscillation, and equatorward cold excursions [3].
The section below plots zonal-mean zonal wind at 60°N against pressure and time, from the same array read as the diagnostic. Each level is standardised against its own winter mean and spread, since the climatological westerlies are an order of magnitude stronger at 10 hPa than near the surface and the tropospheric signal does not survive a shared scale. The black contour is the zero line of the unstandardised field, so the reversal itself stays locatable.
xs = 1:length(wdates)
tickat = findall(d -> day(d) == 1, wdates) # first of each month
xtk = (tickat, string.(monthabbr.(month.(wdates[tickat]))))
cday = findfirst(==(central), wdates)
# Absolute wind is dominated by the climatological westerlies, which are an order of magnitude
# stronger aloft than near the surface — so the descent is invisible in raw m/s. Standardise each
# LEVEL against its own winter mean and spread; the anomaly is then comparable top to bottom.
uanom = similar(ubar)
for k in axes(ubar, 1)
row = @view ubar[k, :]
uanom[k, :] = (row .- mean(row)) ./ std(row)
end
fig = Figure(size = (1000, 660))
# TOP — time–height section. Pressure decreases upward, so the axis is log-scaled and reversed:
# 10 hPa (stratosphere) at the top, 1000 hPa (surface) at the bottom.
ax1 = Axis(fig[1, 1]; ylabel = "pressure (hPa)", yscale = log10, yreversed = true,
yticks = ([10, 30, 100, 300, 1000], ["10", "30", "100", "300", "1000"]),
title = "Zonal-mean zonal wind anomaly at 60°N — the easterly signal descends",
xticks = xtk)
hm = heatmap!(ax1, xs, Float64.(levels), permutedims(uanom);
colormap = Reverse(:RdBu), colorrange = (-2.5, 2.5))
# The black line is the ACTUAL zero-wind contour — the physical boundary, not the normalisation.
contour!(ax1, xs, Float64.(levels), permutedims(ubar); levels = [0.0],
color = :black, linewidth = 1.8)
vlines!(ax1, [cday]; color = (:white, 0.55), linestyle = :dash)
Colorbar(fig[1, 2], hm; label = "σ from winter mean")
hidexdecorations!(ax1; ticks = false, grid = false)
# BOTTOM — the WMO diagnostic itself: 10 hPa, 60°N. Below zero is a major warming.
ax2 = Axis(fig[2, 1]; xlabel = "winter 2020–21", ylabel = "m/s",
title = "The WMO criterion: zonal-mean zonal wind at 10 hPa, 60°N", xticks = xtk)
band!(ax2, xs, min.(u10, 0), zeros(length(xs)); color = (:steelblue, 0.45)) # only where easterly
lines!(ax2, xs, u10; color = :orange, linewidth = 2)
hlines!(ax2, [0]; color = (:white, 0.6), linewidth = 1)
vlines!(ax2, [cday]; color = (:white, 0.55), linestyle = :dash)
text!(ax2, cday + 3, maximum(u10) * 0.72;
text = "reversal · $(central)\n$(length(reversed)) easterly days",
color = (:white, 0.85), fontsize = 12, align = (:left, :top))
ylims!(ax2, minimum(u10) - 6, maximum(u10) + 6)
linkxaxes!(ax1, ax2)
rowsize!(fig.layout, 2, Relative(0.3))
figpng(fig)
The requests the client emits
No constraint expression appears above; netCDF translated the array indexing into requests.
HTTP.VERBOSE in netCDF's rc file routes libcurl's log to standard error. The C library writes to file descriptor 2, below where Julia's redirect_stderr operates, so the descriptor is exchanged for the duration of the call and read back after.
Neither server sends content-length — these responses are chunked — so transfer sizes are counted from the constraint expression instead, which states how many elements it names. A CE with no index ranges names the whole variable.
For DAP2 the sequence is .dds, .das, the coordinate variables, then .dods carrying the hyperslab for the indices used.
# netCDF's libcurl log goes to C stderr (fd 2), below the level Julia's redirect_stderr operates on.
# Exchange the descriptor for the duration of the call, then restore it.
unpct(s) = replace(s, r"%[0-9A-Fa-f]{2}" => m -> string(Char(parse(UInt8, m[2:3], base = 16))))
shorten(u) = (p = split(u, "?", limit = 2); length(p) == 1 ? basename(p[1]) : basename(p[1]) * "?" * p[2])
function trace_requests(f)
path, io = mktemp()
saved = ccall(:dup, Cint, (Cint,), 2)
ccall(:dup2, Cint, (Cint, Cint), Base.cconvert(Cint, Base.fd(io)), 2)
local value
try
value = f()
ccall(:fflush, Cint, (Ptr{Cvoid},), C_NULL)
finally
ccall(:dup2, Cint, (Cint, Cint), saved, 2)
ccall(:close, Cint, (Cint,), saved)
close(io)
end
gets = [unpct(strip(replace(l, r"^> GET " => "", r" HTTP/[\d.]+$" => "")))
for l in eachline(path) if startswith(l, "> GET")]
(value = value, requests = gets)
end
# These responses are chunked, so there is no content-length to read. Element counts come from the
# constraint expression itself, which is exact. A CE with no index ranges names the whole variable.
function ce_elements(ce, fullshape)
rngs = [m.match for m in eachmatch(r"\[[^\]]*\]", ce)]
isempty(rngs) && return prod(fullshape)
n = 1
for r in rngs
parts = split(strip(r, ['[', ']']), ':')
n *= length(parts) == 1 ? 1 :
length(parts) == 2 ? parse(Int, parts[2]) - parse(Int, parts[1]) + 1 :
(parse(Int, parts[3]) - parse(Int, parts[1])) ÷ parse(Int, parts[2]) + 1
end
n
end
# A small read — one day, one level, a 4×3 corner — to keep the logged CE short enough to read.
DAY = "https://psl.noaa.gov/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/pressure/air.2021.nc"
t2 = trace_requests() do
ds = NCDataset(DAY)
patch = ds["air"][1:4, 1:3, HPA10, 5]
close(ds)
size(patch)
end
data_req = last(t2.requests)
(read_shape = t2.value, n_requests = length(t2.requests),
elements_requested = ce_elements(data_req, size(T)),
urls = shorten.(t2.requests))(read_shape = (4, 3), n_requests = 5, elements_requested = 12, urls = ["air.2021.nc.dds", "air.2021.nc.das", "air.2021.nc.dds", "air.2021.nc.dods?level,lat,lon,time,time_bnds", "air.2021.nc.dods?air.air[4][16][0:2][0:3]"])
DAP2 and DAP4 compared
netCDF selects DAP4 from a #dap4 fragment or a dap4:// scheme, with nothing else in the calling code changing — so the two are directly comparable on the same library, dataset and indexing expression.
They translate differently. DAP2 sends the hyperslab. DAP4 sends dap4.ce=/air, the fully-qualified variable name with no index ranges, transfers the variable, and subsets locally. Both return the same twelve values.
The table runs it against a THREDDS endpoint and a Hyrax endpoint, on variables of 14.6 MB and 35.6 MB. Sizes come from the emitted constraint expression; elapsed time is shown alongside but says little, since a fast link moves tens of megabytes in about a second.
Requesting a 14.6 MB variable as one DAP4 response intermittently fails against THREDDS with an HTTP/2 framing error. The probe traces outside the try, so the request is recorded either way and any failure lands in the last column.
# The identical read, over both protocols, against two different servers — NOAA's THREDDS and
# NASA's Hyrax, OPeNDAP's own server. The ONLY thing that changes is the URL.
HYRAX = "https://oceandata.sci.gsfc.nasa.gov/opendap/MODISA/L3SMI/2023/0101/AQUA_MODIS.20230101.L3m.DAY.SST.sst.9km.nc"
SURF = "https://psl.noaa.gov/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/surface/air.sig995.2021.nc"
SERVERS = [
(name = "NOAA PSL · THREDDS", var = "air",
dap2 = SURF, dap4 = replace(SURF, "/dodsC/" => "/dap4/") * "#dap4"),
(name = "NASA OB.DAAC · Hyrax", var = "sst",
dap2 = HYRAX, dap4 = HYRAX * "#dap4"),
]
human(b) = b ≥ 2^20 ? string(round(b / 2^20, digits = 1), " MB") :
b ≥ 2^10 ? string(round(b / 2^10, digits = 1), " KB") : string(b, " B")
# Trace OUTSIDE the try, so we still see what was requested even when the transfer fails.
function protocol_probe(url, var)
tr = trace_requests() do
patch = nothing; err = nothing; shape = ()
t = @elapsed try
ds = NCDataset(url)
shape = size(ds[var])
patch = ndims(ds[var]) == 3 ? ds[var][1:4, 1:3, 1] : ds[var][1:4, 1:3]
close(ds)
catch e
err = first(sprint(showerror, e), 55)
end
(secs = round(t, digits = 1), got = patch === nothing ? 0 : length(patch),
shape = shape, err = err)
end
req = filter(r -> occursin("?", r) && occursin(var, r), tr.requests)
ce = isempty(req) ? "—" : shorten(last(req))
(tr.value..., request = ce,
asked = ce == "—" ? 0 : ce_elements(ce, tr.value.shape))
end
rows = NamedTuple[]
for s in SERVERS, (proto, url) in (("DAP2", s.dap2), ("DAP4", s.dap4))
p = protocol_probe(url, s.var)
push!(rows, (server = s.name, protocol = proto,
variable = human(prod(p.shape) * 4),
values_wanted = 12,
values_requested = p.asked,
on_the_wire = human(p.asked * 4),
secs = p.secs,
values_returned = p.got,
note = something(p.err, "")))
end
slate_table(rows; align = (server = :left, note = :left))| server | protocol | variable | values_wanted | values_requested | on_the_wire | secs | values_returned | note |
|---|---|---|---|---|---|---|---|---|
| NOAA PSL · THREDDS | DAP2 | 14.6 MB | 12 | 12 | 48 B | 0.9 | 12 | |
| NOAA PSL · THREDDS | DAP4 | 14.6 MB | 12 | 3836880 | 14.6 MB | 30.9 | 0 | NetCDF error: NetCDF: libcurl failure (NetCDF error cod |
| NASA OB.DAAC · Hyrax | DAP2 | 35.6 MB | 12 | 12 | 48 B | 1.4 | 12 | |
| NASA OB.DAAC · Hyrax | DAP4 | 35.6 MB | 12 | 9331200 | 35.6 MB | 2.2 | 12 |
Summary
The analysis in the first half is ordinary array code. The wind reversal of 5 January 2021 came out of about 1.4 MB of a 249 MB variable, with fill values and calendar handling supplied by the library and two yearly files presented as one time coordinate. netCDF's DAP2 client turns array indexing into constraint expressions faithfully, underneath tooling with no protocol awareness at all.
The DAP4 path translates the same calls differently: it requests the whole variable and subsets locally — 48 bytes against 14.6 MB on THREDDS and 35.6 MB on Hyrax, for the same twelve values. It reproduces on both servers, so it sits in netCDF's DAP4 handler rather than in either server, and since the values returned are correct it is only visible in a request log.
Worth knowing when planning access patterns through this stack, and the rc file makes it observable without leaving the library.
References
- Kalnay, E. and Kanamitsu, M. and Kistler, R. and Collins, W. and Deaven, D. and Gandin, L. and Iredell, M. and Saha, S. and White, G. and Woollen, J.. The NCEP/NCAR 40-Year Reanalysis Project. 1996.
- Butler, Amy H. and Seidel, Dian J. and Hardiman, Steven C. and Butchart, Neal and Birner, Thomas and Match, Aaron. Defining Sudden Stratospheric Warmings. 2015.
- Baldwin, Mark P. and Dunkerton, Timothy J.. Stratospheric Harbingers of Anomalous Weather Regimes. 2001.