diff --git a/landingpage/shared-sources.toml b/landingpage/shared-sources.toml index 945b000..f0f9738 100644 --- a/landingpage/shared-sources.toml +++ b/landingpage/shared-sources.toml @@ -11,7 +11,7 @@ [source.shared] repo = "freva-org/grid-doctor" ref = "main" -commit = "6d891de1176a075c426f41a20776ec5d94ffd3eb" +commit = "a96e532b799df46e86542c1491a96f607cddc1e7" source_dir = "docs/shared" target_dir = "shared" diff --git a/landingpage/shared.lock b/landingpage/shared.lock index f9cae8d..0b0aed1 100644 --- a/landingpage/shared.lock +++ b/landingpage/shared.lock @@ -4,8 +4,8 @@ # tell a hand edit from a sync without going to the network. [shared] -commit = "6d891de1176a075c426f41a20776ec5d94ffd3eb" +commit = "a96e532b799df46e86542c1491a96f607cddc1e7" [shared.files] -"high-level-access.md" = "0ab36941e00f9fca3f36c2e0ca286cc760bff9910d8fa55c21e1e2b5cce6bc1d" -"technical-decisions.md" = "e2ef22035790186b2bb22177275f7f8d329e73fe8e65231b2e1ec33010411757" +"high-level-access.md" = "9dda1b30221272bd71247b61f0a420763d81502176feec4c2028e9c008d02469" +"technical-decisions.md" = "415aab02d8dda8ba1636f48c4ad06f916c5cc1db91c956157eeb495572e5931c" diff --git a/landingpage/shared/high-level-access.md b/landingpage/shared/high-level-access.md index 5396a59..60a961d 100644 --- a/landingpage/shared/high-level-access.md +++ b/landingpage/shared/high-level-access.md @@ -1,10 +1,10 @@ - # Accessing High-Resolution Regional Data -Some datasets in the hub cover a small area (a city, an island, a +Some datasets cover a small area (a city, an island, a catchment) at very high HEALPix levels (16 is roughly 100 m, 20 roughly 6 m cells). These stores look different from the global pyramids in two ways, and both change how you load them: @@ -143,6 +143,17 @@ coordinates and the global HEALPix indices in `berlin["cell"]`: berlin["t2m"].mean() ``` +!!! note "Fields with missing values" + For variables with NaNs (ocean-only fields, observation gaps), cells + on coarser levels can be only partly valid, and a plain mean over + cells is biased. If the store has a `_valid_fraction` + variable, use it as the weight: + + ```python + ds.sst.weighted(ds.sst_valid_fraction.fillna(0)).mean("cell") + ``` + + ### Circular regions Select a bounding box around the centre as above, then mask by @@ -206,7 +217,7 @@ era5_region = era5_level9.isel( ``` No interpolation, no index join. The hierarchies coincide by -construction, which is one of the reasons the hub uses nested ordering +construction, which is one of the reasons to use nested ordering everywhere. ## Coarser overview levels @@ -214,7 +225,7 @@ everywhere. High-level stores are still pyramids. For a quick overview, open a coarse level of the *same* store. Levels at or below the coordinate threshold are written with materialised coordinates and behave exactly -like every other dataset in the hub, including in the browser viewer: +like every other global pyramid, including in the browser viewer: ```python overview = xr.open_zarr(".../city-example.zarr/level_8.zarr", chunks=None) diff --git a/landingpage/shared/technical-decisions.md b/landingpage/shared/technical-decisions.md index 1e1007e..d79a23e 100644 --- a/landingpage/shared/technical-decisions.md +++ b/landingpage/shared/technical-decisions.md @@ -1,4 +1,4 @@ - @@ -23,6 +23,13 @@ arithmetic mean of all pixels. No latitude-dependent cosine weighting is needed, and an entire class of subtle statistics bugs, forgotten or misapplied area weights, simply cannot occur. +!!! note "Fields with missing values" + This holds for complete fields. For fields with NaNs (ocean-only + variables, sea ice, observation gaps), a coarse cell can be partly + valid, and its value is the mean over that valid part only. Averages + over coarse levels then need the valid fraction of each cell as a + weight; + ![Equal area comparison](assets/healpix-equal-area.png#only-dark) ![Equal area comparison](assets/healpix-equal-area-light.png#only-light) /// caption @@ -75,179 +82,71 @@ that mode every pyramid level is independently remapped from the source grid, which is substantially more expensive. -## Spherical geometry - -All HEALPix cell boundaries and centre coordinates are computed on a -**perfect sphere**, not on the WGS84 ellipsoid[^wgs84]. The remapping -software we use, -[ESMF](https://earthsystemmodeling.org/regrid/) (Earth System Modeling -Framework, a widely used library for regridding Earth-system model -output), computes conservative overlap areas exclusively on the sphere, -with great-circle cell edges; it has no ellipsoidal geometry -mode[^esmf]. The HEALPix target mesh, the source mesh, and the overlap -calculation therefore all use one and the same geometry. - -**Plain-language summary.** It is tempting to think that a spherical -pipeline is "off by up to 21 km" compared to the real, slightly -flattened Earth, and that switching the grid definition to the WGS84 -ellipsoid would fix this. Neither is true. As long as *every* step of -the pipeline uses the same convention, the sphere approximation cancels -out of the remapping almost entirely: the residual error is a few parts -per million. The 14 to 21 km error only becomes real when conventions -are *mixed*, that is, when an ellipsoid-defined grid is fed into a -sphere-based remapping tool (or vice versa). Adopting ellipsoidal -HEALPix today, while ESMF and the rest of the regridding ecosystem -remain spherical, would create exactly that mixed state. It would -*introduce* a systematic geolocation error of up to about 14 km, many -pixel widths at kilometre-scale resolution, into a change whose stated -motivation is geolocation accuracy. Consistency beats fidelity. - -### Why the sphere approximation cancels - -Everyone's data carries geodetic (WGS84) latitude $\varphi$. Our -pipeline interprets $\varphi$ directly as spherical latitude, a smooth -and invertible convention applied identically to the source grid and -the HEALPix target. Because both grids are transformed the same way, -their overlap topology in $(\lambda, \varphi)$ space is preserved -exactly; a geodetic position round-trips to the same HEALPix cell -without loss. - -The only distortion is the *area measure*. On the sphere the area -element is $dA_s = R^2 \cos\varphi \, d\varphi \, d\lambda$; on the -ellipsoid it is $dA_e = M(\varphi) N(\varphi) \cos\varphi \, d\varphi -\, d\lambda$, with the meridional and prime-vertical radii of -curvature[^snyder] - -$$ -M(\varphi) = \frac{a (1 - e^2)}{(1 - e^2 \sin^2\varphi)^{3/2}}, -\qquad -N(\varphi) = \frac{a}{(1 - e^2 \sin^2\varphi)^{1/2}}, -$$ - -and WGS84 eccentricity $e^2 = f(2-f) \approx 0.006694$[^wgs84]. Write -the distortion ratio $\rho(\varphi) = dA_e / dA_s$. A first-order -conservative remapping weight is a **ratio of areas**[^jones], - - -$$ -w_{ij} \;=\; \frac{A(S_i \cap T_j)}{A(T_j)}, -$$ - - -where $S_i$ is cell $i$ of the source grid (for example an ICON -triangle or a satellite pixel footprint), $T_j$ is target HEALPix cell -$j$, and $A(\cdot)$ denotes area. The weight $w_{ij}$ is the fraction -of the target cell covered by that source cell, and the remapped field -is the weighted average $\tilde f_j = \sum_i w_{ij} f_i$. Numerator -and denominator sample $\rho$ at the *same latitude*, so the distortion -cancels to first order: - +## Sphere or ellipsoid + +HEALPix can be defined on a perfect sphere, the original +construction[^gorski], or on the WGS84 ellipsoid[^wgs84], where the cell +boundaries are laid out in authalic latitude so that cells stay +equal-area on the ellipsoid[^gibb][^hpgeo]. grid-doctor uses the sphere +by default. + +Both definitions give the same results with +[ESMF](https://earthsystemmodeling.org/regrid/), even though ESMF +computes all overlap areas on the sphere[^esmf]. healpix-geo returns the +cell vertices of either definition in geodetic coordinates, the same +coordinates the source data uses, so source and target mesh always +share one coordinate space; the ellipsoidal definition only moves where +the cell boundaries lie. The remaining distortion of the sphere +approximation affects the numerator and denominator of every +conservative weight[^jones] almost equally and cancels to a relative +error of about $10^{-5}$ at level 10, smaller at finer levels. + +What does not work is converting a remapped product from one definition +to the other. Re-binning already averaged cells onto displaced cells +adds a second interpolation step and roughly doubles the error. +Products for either definition should always be remapped from the +native data. + +| Route | Product | RMS error | × local contrast | +|:--|:--|--:|--:| +| A | native → spherical HEALPix | 0.0840 K | 0.94 | +| B | native → ellipsoidal (WGS84) HEALPix | 0.0841 K | 0.94 | +| C | A converted to ellipsoidal | 0.1639 K | 1.84 | -$$ -w_{ij}^{\mathrm{sphere}} -= \frac{\displaystyle\int_{S_i \cap T_j} \rho^{-1}\, dA_e} - {\displaystyle\int_{T_j} \rho^{-1}\, dA_e} -= w_{ij}^{\mathrm{ellipsoid}} - \left( 1 + \mathcal{O}\!\left(\tfrac{\rho'}{\rho}\, - \Delta\varphi_{\mathrm{cell}} \right) \right). -$$ - -The residual scales with the *variation* of $\rho$ across a single -cell, $\rho'/\rho \sim 2 e^2 \sin 2\varphi \approx 0.013$ per radian. -For a level-10 HEALPix cell ($\Delta\varphi_{\mathrm{cell}} \approx 5.4 -\times 10^{-4}$ rad) the relative weight error is - -$$ -\varepsilon_{\mathrm{consistent}} -\;\sim\; 2 e^2 \, \Delta\varphi_{\mathrm{cell}} -\;\approx\; 7 \times 10^{-6}, -$$ - -and it *shrinks* with increasing resolution. What does not cancel is -the absolute area of a cell: integrating with spherical cell areas over -ellipsoidal data carries a smooth latitude-dependent bias -$\rho(\varphi) - 1$ of roughly $-0.4\,\%$ at the equator to $+0.9\,\%$ -at the poles. This affects diagnostics, not remapping, and is fully -correctable in post-processing by integrating with true ellipsoidal -cell areas. +/// caption +One day of MODIS-Aqua night-time SST (4 km) remapped to HEALPix +level 10 (about 6.4 km) with ESMF conservative weights, about 2.66 +million ocean cells. Each route is scored against a reference that bins +the native pixels directly into cells of the same definition. The local +contrast (median spread within groups of four neighbouring cells, +0.089 K) is the relevant scale because remapping errors live at the cell +scale. A and B agree within 1 % between 60° S and 45° N and within 5 % +poleward of that. +/// -### Why mixing conventions does not cancel +### Different definitions are different grids -An ellipsoidal HEALPix grid keeps its cells equal-area on the ellipsoid -by defining cell boundaries through an auxiliary latitude[^snyder]. -Suitable choices are the authalic latitude $\xi$ (the construction used -by the rHEALPix DGGS on the ellipsoid[^gibb] and by healpix-geo's WGS84 -support[^hpgeo]) or, for pure geolocation, the geocentric latitude -$\psi$: +Two datasets on HEALPix at the same level are not on the same grid +unless they use the same definition. The same cell index denotes a +different patch of the Earth under the two definitions, displaced by +the difference between geodetic and authalic latitude[^snyder]: $$ \varphi - \xi \approx \tfrac{e^2}{3} \sin 2\varphi -\quad (\text{max } 7.7' \approx 14\ \mathrm{km}), -\qquad -\varphi - \psi \approx \tfrac{e^2}{2} \sin 2\varphi -\quad (\text{max } 11.5' \approx 21\ \mathrm{km}), -$$ - -both maximal at $\varphi = 45°$. Feed such a grid into a spherical -remapper together with a source grid whose coordinates are plain -geodetic latitudes, and the two grids are now expressed in *different* -latitude conventions: the target is displaced relative to the source by -$\delta\varphi(\varphi)$. Conservative overlaps are then computed -between systematically misregistered polygons, and the remapped field -is, to first order, the true field evaluated at a shifted position: - -$$ -\tilde f(\mathbf{x}) \;\approx\; f(\mathbf{x} + \boldsymbol{\delta}), -\qquad |\boldsymbol{\delta}| \lesssim 14\text{ to }21\ \mathrm{km}. -$$ - -Unlike the consistent case, this error is **first order** in $e^2$ and -is a *displacement*, the worst kind of error for high-resolution data: -HEALPix cells measure roughly 6.4 km at level 10, 1.6 km at level 12, -and 0.8 km at level 13, so the misregistration amounts to **10 to 20 -cell widths** at kilometre scale. Every front, coastline, and swath -edge lands in the wrong cells. Comparing the two regimes at level 12: - -$$ -\frac{\varepsilon_{\mathrm{mixed}}}{\varepsilon_{\mathrm{consistent}}} -\;\sim\; -\frac{e^2/3}{2 e^2 \, \Delta\varphi_{\mathrm{cell}}} -\;=\; -\frac{1}{6\,\Delta\varphi_{\mathrm{cell}}} -\;\approx\; 10^3 . +\quad (\text{max } 7.7' \approx 14\ \mathrm{km\ at}\ 45°), $$ -Switching the grid definition without switching the remapping stack -makes the error roughly **three orders of magnitude larger** than -staying consistently spherical. - -### Migration path - -Ellipsoidal HEALPix is not wrong, it is premature. The correct -construction defines the grid through the authalic latitude, following -the rHEALPix ellipsoidal DGGS[^gibb] and healpix-geo[^hpgeo]. Because -the authalic mapping is equal-area by construction[^snyder], a -spherical conservative remapping performed in authalic coordinates is -*exactly* area-preserving on the ellipsoid; the only residual is the -difference between great-circle edges and images of ellipsoidal -geodesics, which is $\mathcal{O}(e^4)$. A correct ellipsoidal pipeline -therefore requires the transform $\varphi \to \xi(\varphi)$ to be -applied to the source mesh before weight generation, either inside the -remapping tool (an ellipsoid-aware ESMF or equivalent) or consistently -in every tool that touches the data, including downstream cell lookups. - -A half-migrated ecosystem is strictly worse than either pure -convention. Our commitment is therefore: - -1. **Today**: all geometry (HEALPix boundaries[^gorski], source - meshes, ESMF overlap calculation) stays on the perfect sphere, - with geodetic latitude interpreted as spherical latitude end to end. -2. **When** an ellipsoid-aware conservative remapping tool is available - and validated (a new ESMF version or an equivalent that applies the - authalic transform internally), we switch the pipeline **and remap - all published datasets** in one coordinated step, so that no dataset - ever mixes conventions with another. +with WGS84 eccentricity $e^2 \approx 0.006694$. That is about 2 cell +widths at level 10, 9 at level 12 and 18 at level 13; at level 10, +99.8 % of all cells have a different index under the two definitions. +Data on different definitions must therefore not be compared, overlaid +or differenced cell by cell, and cells must not be located with the +other definition (for example `lonlat_to_healpix` with a different +ellipsoid than the data was produced with). Otherwise the field appears +shifted by up to 14 km: fronts, coastlines and swath edges land in the +wrong cells. To compare such datasets, remap both from their native +data onto the same definition. ## Target level selection @@ -427,16 +326,15 @@ excluded because they do not commute with mean coarsening in the pyramid; if an RMS quantity is required, bin the *squared* field with the mean and apply the square root at read time. -### Sphere, not ellipsoid — also for satellite data +### Sphere, also for satellite data -Satellite geolocation is geodetic (WGS84). Binning nevertheless indexes -cells on the **perfect sphere**, exactly like all remapped datasets. The -geodetic-vs-spherical latitude discrepancy (up to ~0.19° at 45° latitude, -~21 km) is accepted so that every dataset in the hub shares a single -indexing geometry. Indexing one dataset on the ellipsoid while all others -use the sphere would shift it by dozens of pixels at typical swath target -levels, producing visible misregistration (offset coastlines) when -datasets are overlaid. +Satellite geolocation is geodetic (WGS84). Binning indexes cells on the +**perfect sphere**, like remapping, by assigning each geodetic position +to the cell that contains it; this is as accurate as binning on the +ellipsoid (see [Sphere or ellipsoid](#sphere-or-ellipsoid)). Data +binned on one definition is not comparable cell by cell with data on the +other: the cells are displaced by up to ~0.13° (~14 km) at 45° latitude, +several to dozens of pixels at typical swath target levels. ### Minimum sample count and coverage tracking @@ -463,6 +361,9 @@ and ``"nearest"`` do for remapped data. The multi-resolution pyramid is built by first remapping the source dataset to the finest HEALPix level, then deriving all coarser levels by iterated coarsening, one level at a time, always a factor of 4. +Every level builds on the previous one, so the whole pyramid is a single +lazy computation (see +[Chunked application and single-pass writes](#chunked-application-and-single-pass-writes)). ### Why coarsen rather than remap at each level? @@ -477,7 +378,26 @@ coarsening hierarchy. ### Mean coarsening (continuous fields) For continuous fields (those remapped with conservative weights), each -parent cell's value is the NaN-aware mean of its 4 children. +coarse cell's value is the mean over **all valid finest-level cells** +beneath it. This is not the same as the mean of its 4 children once +children are only partly valid: a child built from a single valid +finest-level cell would otherwise count as much as a child built from +four. Repeated over several levels, that mean of means drifts away from +the finest level. + +The pyramid therefore carries two running totals per cell instead of a +mean: the sum of the valid finest-level values and their count. Both +coarsen exactly by summing groups of 4, and the mean is only formed when +a level is written: + +$$ +\bar{x}_\text{parent} = \frac{\sum_{c} s_c}{\sum_{c} n_c}, +\qquad s_c = \sum_{\text{valid } i \in c} x_i, \quad n_c = \#\{\text{valid } i \in c\} +$$ + +Every level is therefore identical to coarsening directly from the +finest level, while each step stays a local reduction of 4 neighbouring +cells. ### Mode coarsening (categorical fields) @@ -491,17 +411,28 @@ stored in the dataset attributes: `grid_doctor_method = "nearest"` triggers mode coarsening, `"conservative"` triggers mean coarsening. An explicit `coarsen_mode` parameter is available to override this. +Unlike the mean, a mode of modes is not the mode of all finest-level +cells, and mode coarsening is applied level by level, with the minimum +valid fraction checked against the 4 children at each step. + ### Minimum valid fraction (default: 50%) -A parent cell is set to NaN when fewer than half of its children are -valid (at least 2 of 4). +For mean coarsening, a cell is set to NaN when fewer than half of the +**finest-level** cells beneath it are valid. The threshold is applied to +this cumulative fraction, not to the 4 children of each step. + +**Representativeness.** A cell's value should represent the majority +of its area. Because the threshold uses the finest-level count, the +value is guaranteed to cover at least half of the cell's area at every +level. A per-step rule could not guarantee this: 2 valid children out +of 4, each built from 2 valid children out of 4, cover only a quarter of +the area two levels up. Below the threshold, a few pixels' values would +"speak for" an area that has no data, which is indistinguishable from +interpolation into unknown territory. -**Representativeness.** A parent cell's value should represent the -majority of its area. With at least 2 of 4 valid children, the value -is guaranteed to cover at least half the parent cell's area. Below -that, a single child pixel's value would "speak for" 3 other pixels -that have no data, which is indistinguishable from interpolation into -unknown territory. +Cells masked by the threshold keep their running sums and counts, so +their valid data still contributes to coarser levels where it is part of +a majority. **Cascade prevention.** Each coarsening step is a factor-of-4 reduction. If only 1 of 4 valid children were sufficient, a single @@ -516,6 +447,55 @@ the very first step because 1/4 < 1/2. A feature can only survive coarsening if it covers at least half the area at every scale, which is exactly when it is a real, resolvable feature at that resolution. +The threshold is configurable with `min_valid_fraction`. Set it to `0` +to keep every cell that contains any valid data, for instance when the +coarse levels are only used for weighted aggregates. + +### Averaging over coarse levels: valid fractions + +Because every coarse value is a mean over its valid area, a plain mean +over coarse cells weights a coastal cell that is 10% ocean as much as an +open-ocean cell. The level-to-level consistency of the pyramid only +holds for means weighted by the number of valid finest-level cells: +with two coarse cells, one fully valid at 10 °C (4 valid cells) and one +with a single valid cell at 20 °C, the finest-level mean is +(4 × 10 + 20) / 5 = 12 °C, while the plain mean of the two coarse values +is 15 °C. Both coarse values are correct; the plain average is not. + +The weights cannot be folded into the stored values without losing their +meaning (a coastal cell would no longer hold the SST of its ocean +part). Instead, the pyramid can store them alongside: + +```python +pyramid = gd.create_healpix_pyramid(ds, valid_fraction=True) +``` + +adds `_valid_fraction` (float32, fraction of valid finest-level +cells, linked through the CF `ancillary_variables` attribute) to every +level, including the finest one. Global or regional means then agree +across all levels: + +```python +ds = pyramid[3] +ds.sst.weighted(ds.sst_valid_fraction.fillna(0)).mean("cell") +``` + +The option is opt-in because the fraction has the shape of its variable: + +- `True` stores a fraction for every cell variable, with its full shape. + This is correct for masks that change over time or height (sea ice, + clouds, orography on pressure levels). +- `"static"` stores a single `cell`-only fraction taken from the first + time step / level. It is tiny, but only correct when the mask never + changes, as for a land-sea mask. +- A list of names, or a mapping such as `{"sst": "static", "ice": True}`, + restricts the fractions to selected variables. + +Fractions are mostly exactly 0 or 1 and compress well; in a test with 8 +time steps the full-shape fractions added about 7% to the store, the +static ones under 1%. This mirrors the `ocean_fraction` variables of +the nextGEMS HEALPix output. + ## Output metadata and CRS convention @@ -530,7 +510,9 @@ Every output dataset carries a standardised set of metadata. | `healpix_order` | `nested` or `ring` | | `grid_doctor_version` | Package version that produced the data | | `grid_doctor_method` | `conservative` or `nearest` | -| `grid_doctor_coarsened_from_level` | Immediate parent level (coarsened levels only) | +| `grid_doctor_coarsened_from_level` | Level the values were coarsened from: the finest level for pyramids, the input level for `coarsen_healpix` (coarsened levels only) | +| `grid_doctor_coarsen_mode` | `mean` or `mode` (coarsened levels only) | +| `grid_doctor_min_valid_fraction` | Threshold below which cells were set to NaN: cumulative over finest-level cells for `mean`, per coarsening step for `mode` (coarsened levels only) | ### CRS variable @@ -593,7 +575,7 @@ returned. Combined with `fill_value = NaN` and empty-chunk elision on the data variables, storage and access cost are proportional to the domain, not the globe, while cross-dataset and cross-level alignment stay a bit-shift (`parent = id >> 2k`) exactly as for every other -dataset in the hub. +nested HEALPix dataset. Levels at or below the threshold within the same pyramid are written with materialised coordinates as usual, so coarse overview levels stay @@ -625,6 +607,28 @@ amortised over the batch. The backend can be overridden explicitly via `backend="scipy"`, `"numba"`, or `"cupy"` in any remapping call. +### Chunked application and single-pass writes + +For dask-backed input, the regridded field is chunked along `cell` +(`cell_chunks`, by default 4^10 cells). Rows of the weight matrix are +independent, so each output chunk is computed from its own block of +matrix rows; the blocks share the matrix memory instead of copying it. +No task holds the whole HEALPix field, and the coarser levels keep the +same chunking because every coarsening step only combines groups of 4 +neighbouring cells. + +Inside dask tasks the Numba kernels run single-threaded: dask already +runs many tasks in parallel, and a multi-threaded kernel per task would +start one thread team per task and exhaust the process thread limit on +many-core nodes. Without dask, the multi-threaded kernels are used. + +`save_pyramid` writes all levels in a single computation, so the +regridding of each finest-level chunk runs exactly once and is shared by +all coarser levels. Computing levels separately (for example +`pyramid[5].compute()` followed by `pyramid[4].compute()`) regrids the +source again each time; `persist()` the finest level when exploring +interactively. + ---