eduaguilera / eduaguilera/whep

Land-use spatializer writes denormal float32 values (smallest is 1.4e-45)

Open
#985 1 comment 0 reactions 1 assignee Claimed by @lbm364dl View on GitHub
area:spatialize bug needs-expert priority:medium
Dominant language
R
Stars
1
Forks
5
Avg merge
1d 11h
Merged PRs (30d)
183

Description

The gridded land-use file written by `prepare_spatialize_all.R` contains crop
fractions far below any physically meaningful area, down to and including
denormal float32 values.

Found while profiling LPJmL runtime — a third of the crop stands LPJmL simulates
come from these entries.

## What is in the file

`LPJmL_inputs/whep/lpjml_inputs/landuse/cft_default_cft_aggregation_30min_1750-2023.nc`,
year 2023, 46,851 land cells, 426,600 nonzero crop entries:

| | |
|---|---|
| smallest nonzero value | **1.401e-45** (= `FLT_TRUE_MIN`, smallest positive denormal float32) |
| entries below 1e-15 (≈ µm² of the cell) | 1.11% |
| entries below 1e-9 (≈ 2.5 m²) | 3.16% |
| entries below 4e-6 (**less than one hectare**) | **21.9%**, holding 0.0004% of cropped area |
| median nonzero entry | 2.9e-04 |

Cells carry 9.1 crops each on average, out of 32.

## Where they come from

In `R/spatialize.R`, `.spatialize_year()` joins country crop areas to grid cells
cartesian — the comment says so directly:

```r
# Join country_areas (cartesian: each country-crop gets its cells)
dat <- grid_cp[ca, on = .(area_code, item_prod_code), nomatch = NA]
```

then allocates proportionally:

```r
allocated_rf = rf_potential / rf_pot_sum * rainfed_target
```

so every crop a country reports gets a share in every cropland cell of that
country. The only filter is

```r
result <- dat[allocated_rf > 0 | allocated_ir > 0, ...]
```

which drops exact zeros only. Anything that underflows toward zero without
reaching it survives into the file, and float32 storage then carries it down
into the denormal range.

## Why it matters

1. **It is not data.** A value of 1e-45 of a 2,500 km² cell is not a land area;
it is arithmetic residue. Reporting a crop as "present" in a cell on that
basis is misleading.
2. **It is expensive downstream.** LPJmL gives every crop in a cell its own soil
column and full daily water, nitrogen and phenology cycle. A third of all
crop stands in a WHEP run hold 0.005% of the land between them, and dropping
those below 1e-4 makes the model **18.7% faster** while moving global NPP,
soil carbon and vegetation carbon not at all and transpiration by 0.1%.
3. **Denormals are slow** on some hardware — they can trap to microcode — so
they are worth removing on their own account.

## Suggested fix

Threshold the allocation before writing, rather than filtering on `> 0`. Two
options:

- an absolute floor in hectares, e.g. drop allocations below one hectare
(`4e-6` of a 30-arcmin cell), which removes 21.9% of entries and 0.0004% of
area; or
- a relative floor, e.g. drop anything below 1e-6 of the cell.

Either way the dropped area should be folded back into the same country-crop's
remaining cells so national totals are preserved, rather than silently lost.

Happy to prepare a PR if the approach looks right — the choice of threshold is
yours since it is a data question, but the denormals specifically look like a
straightforward bug.

Contributor guide

Open the contributing guide

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.