Add ASTER GED v3 surface-emissivity dataset

[xkykai]

Adds ASTERGEDv3 — the ASTER Global Emissivity Dataset v3 (NASA LP DAAC): a static 2000–2008 clear-sky-mean land-surface emissivity climatology on a WGS84 lat/lon grid at 1 km (:low_1km, default) or 100 m (:high_100m), as a DataWrangling submodule.

Usage

using NumericalEarth, ArchGDAL   # ArchGDAL + Earthdata credentials activate the HDF5 read/download

region = BoundingBox(longitude = (-112.8, -111.2), latitude = (35.2, 36.8))
metadatum = Metadatum(:emissivity; dataset = ASTERGEDv3(), region)
ε = Field(metadatum, grid)
land_surface = SurfaceRadiationProperties(albedo = α, emissivity = ε)

What it does

• Resolves the intersecting 1°×1° HDF5 tiles via NASA CMR, downloads them from Earthdata Cloud (EARTHDATA_USERNAME/EARTHDATA_PASSWORD; tiles cached for reuse), reads them through GDAL's HDF5 driver, and mosaics + clips to the BoundingBox.
• Collapses the five ASTER TIR bands to one broadband ε with Ogawa & Schmugge (2004) coefficients (they sum to 1, so it is a convex combination), decoding in Float32 (mean ×0.001, SDev ×0.0001). The decoded broadband field is cached as a regional NetCDF.
• Returns a reduced (Center, Center, Nothing) Field — safely stateindex-able at any k, so it plugs directly into SurfaceRadiationProperties. :emissivity_uncertainty propagates the per-band SDev to a broadband σ(ε).
• ASTERGEDv3 is exported from NumericalEarth; :high_100m/:low_1km map internally to NASA's AG100/AG1KM short names.

Gap filling — surface-aware inpainting

ASTER GED retrieves an emissivity almost everywhere, including over water (open-ocean broadband ε ≈ 0.98 — a real value, not a fill). Genuine gaps appear only where a clear-sky retrieval was never obtained (persistent cloud in the humid tropics, snow screened as cloud) and decode to NaN. Building a Field fills them with NearestNeighborInpainting(Inf), so the field is finite everywhere.

The fill is surface-aware: the tiles' land/water map is read and written alongside the emissivity, so a land gap is filled only from land neighbors and a water gap only from water — a gap never borrows across the coastline (an isolated gap with no reachable same-class donor falls back to nearest-valid so the field stays finite).

This is generic, opt-in DataWrangling machinery: inpaint_mask! / Field take an optional per-cell regions label field, supplied per dataset via the inpainting_regions hook (ASTER GED returns its land/water partition). With no regions (the default) the inpainting is byte-for-byte unchanged, so ocean datasets (ECCO/EN4/WOA) that deliberately extend values across the coastline are unaffected.

Example — 100 m broadband emissivity

Clear-sky desert (Grand Canyon, gap-free) vs. humid tropics (Congo, cloud gaps inpainted). Plots the Fields directly — no NetCDF round-trip:

using NumericalEarth
using NumericalEarth.DataWrangling: BoundingBox
using Oceananigans
using ArchGDAL      # activates the ASTER GED HDF5 read/download extension
using CairoMakie

dataset = ASTERGEDv3(resolution = :high_100m)

function emissivity_field(dataset, name, longitude, latitude; size = (512, 512))
    grid = LatitudeLongitudeGrid(CPU(); size, longitude, latitude,
                                 topology = (Bounded, Bounded, Flat))
    region = BoundingBox(; longitude, latitude)
    return Field(Metadatum(name; dataset, region), grid)
end

gc_emissivity  = emissivity_field(dataset, :emissivity,             (-112.8, -111.2), (35.2, 36.8))
gc_uncertainty = emissivity_field(dataset, :emissivity_uncertainty, (-112.8, -111.2), (35.2, 36.8))
congo_emissivity = emissivity_field(dataset, :emissivity, (18, 20), (-1, 1))

fig = Figure(size = (1650, 520), backgroundcolor = :white)
Label(fig[0, 1:6], "ASTER GED v3 broadband emissivity (100 m)"; fontsize = 19, font = :bold)

ax1 = Axis(fig[1, 1]; title = "Emissivity — Grand Canyon (4-tile mosaic)",
           xlabel = "longitude (°)", ylabel = "latitude (°)", aspect = 1)
hm1 = heatmap!(ax1, gc_emissivity; colormap = :viridis, colorrange = (0.90, 0.98))
Colorbar(fig[1, 2], hm1; label = "broadband ε")

ax2 = Axis(fig[1, 3]; title = "Uncertainty σ(ε) (SDev path)",
           xlabel = "longitude (°)", ylabel = "latitude (°)", aspect = 1)
hm2 = heatmap!(ax2, gc_uncertainty; colormap = :viridis, colorrange = (0, 0.02))
Colorbar(fig[1, 4], hm2; label = "broadband σ(ε)")

ax3 = Axis(fig[1, 5]; title = "Emissivity — Congo basin (cloud gaps inpainted)",
           xlabel = "longitude (°)", ylabel = "latitude (°)", aspect = 1)
hm3 = heatmap!(ax3, congo_emissivity; colormap = :viridis, colorrange = (0.90, 0.98))
Colorbar(fig[1, 6], hm3; label = "broadband ε")

save("asterged_emissivity_map_100m.png", fig)

Notes / follow-ups

• The broadband coefficients are a Sahara-desert regression applied globally; the static clear-sky mean also misses snow (ε 0.03–0.08 higher) and soil-moisture effects.
• Docs @example page: follow-up.

🤖 Generated with Claude Code

[codecov]

Codecov Report

❌ Patch coverage is 48.17814% with 128 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
ext/NumericalEarthArchGDALExt.jl	0.00%	101 Missing ⚠️
src/DataWrangling/ASTERGED/ASTERGED.jl	76.53%	23 Missing ⚠️
src/DataWrangling/inpainting.jl	89.18%	4 Missing ⚠️

📢 Thoughts on this report? Let us know!

[xkykai]

Running the script above gives

[glwagner]

Why did you "save for plotting" to NetCDF? is there some issue plotting the Field?

[glwagner]

should this be done on GPU? Or at least permitted?

[glwagner]

It might make sense to add a section in the docs tutorial for land data + computations that are done here. It can be a future PR though.

[glwagner]

Another note to self is that I am wondering how this will all work in a distributed context, eg if we run 100m resolution globally, we want to do the calculations in here independently on every node. It's ok to put optimizations in a future PR but we should think about it a bit.

[xkykai]

Why did you "save for plotting" to NetCDF? is there some issue plotting the Field?

no there's nothing wrong with it. I was just stitching together two vibe-coded scripts because it'll take a while to rerun. but a script like

using NumericalEarth
using NumericalEarth.DataWrangling: BoundingBox, Metadatum
using Oceananigans
using ArchGDAL  # activates the HDF5 tile read/download extension
using CairoMakie

dataset = ASTERGEDv3()

#####
##### 1. Multi-tile static emissivity map — US Southwest (Grand Canyon / Colorado Plateau)
#####

southwest = BoundingBox(longitude = (-112.8, -111.2), latitude = (35.2, 36.8))
southwest_grid = LatitudeLongitudeGrid(CPU(), Float32;
                                       size = (320, 320),
                                       longitude = (-112.8, -111.2),
                                       latitude = (35.2, 36.8),
                                       topology = (Bounded, Bounded, Flat))

emissivity  = Field(Metadatum(:emissivity; dataset, region = southwest), southwest_grid)
uncertainty = Field(Metadatum(:emissivity_uncertainty; dataset, region = southwest), southwest_grid)

#####
##### 2. Coverage-hole probe — Congo basin (persistent cloud)
#####

congo = BoundingBox(longitude = (18, 20), latitude = (-1, 1))
congo_grid = LatitudeLongitudeGrid(CPU(), Float32;
                                   size = (200, 200),
                                   longitude = (18, 20),
                                   latitude = (-1, 1),
                                   topology = (Bounded, Bounded, Flat))

congo_emissivity = Field(Metadatum(:emissivity; dataset, region = congo), congo_grid)

#####
##### Plot the Fields directly — Oceananigans' Makie extension supplies the coordinates
#####

fig = Figure(size = (1650, 520))
Label(fig[0, 1:6], "ASTER GED v3 AG100 static broadband emissivity"; fontsize = 19, font = :bold)

ax1 = Axis(fig[1, 1]; title = "Emissivity — US Southwest (4-tile mosaic)",
           xlabel = "longitude (°)", ylabel = "latitude (°)", aspect = 1)
hm1 = heatmap!(ax1, emissivity; colormap = :viridis, colorrange = (0.90, 0.98))
Colorbar(fig[1, 2], hm1; label = "broadband ε")

ax2 = Axis(fig[1, 3]; title = "Uncertainty (SDev path)",
           xlabel = "longitude (°)", ylabel = "latitude (°)", aspect = 1)
hm2 = heatmap!(ax2, uncertainty; colormap = :viridis, colorrange = (0, 0.02))
Colorbar(fig[1, 4], hm2; label = "broadband σ(ε)")

ax3 = Axis(fig[1, 5]; title = "Congo basin — coverage holes inpainted, water filled",
           xlabel = "longitude (°)", ylabel = "latitude (°)", aspect = 1)
hm3 = heatmap!(ax3, congo_emissivity; colormap = :viridis, colorrange = (0.90, 0.98))
Colorbar(fig[1, 6], hm3; label = "broadband ε")

save("asterged_emissivity_static_map.png", fig)

would work fine.

[xkykai]

Another note to self is that I am wondering how this will all work in a distributed context, eg if we run 100m resolution globally, we want to do the calculations in here independently on every node. It's ok to put optimizations in a future PR but we should think about it a bit.

Suggestions by Fable 5:

1. Each rank should build Metadatum with its own BoundingBox
2. Derive each rank's BoundingBox automatically from its local grid (plus an interpolation margin).
3. Each rank must fetch its own dataset.
4. A persistent granule-keyed tile cache instead so neighboring ranks sharing a boundary tile don't re-download it.
5. Gap inpainting is window-local: a rank whose window lies entirely inside a retrieval gap (humid-tropics cloud holes) has no valid source — either exchange halos between inpainting iterations or pre-fill gaps from the "1km" product.