The alarming rate of climate change and of nature loss is increasing the pressure on countries to rethink how their land is managed. Land use already accounts for about a quarter of greenhouse-gas emissions, is the largest driver of biodiversity loss, and decides the security of food and water [1, 2]. Deciding what to farm, protect or restore needs models of the landscape that are detailed, cover whole countries, and are honest about their uncertainty, and every such model starts from how the land is represented.

There are two fundamental ways to represent land data in geography. The first is the vector form, where points, lines and polygons describe what is on the ground. It is common for registries, such as parcels, buildings and roads, and for interactive maps. The second is the raster form, where the land is cut into a regular grid of pixels. Rasters offer a universal way to combine data from different sources, since every layer aligns on the same grid, and they are the natural format of satellite imagery. As a result, most modeling of land use and land cover (LULC) is done on rasters [3].

But raster pixels are not independent observations. As Tobler stated in his first law of geography [4], near things are more related than distant things. At the scale of pixels, this means that neighboring pixels often belong to the same field, forest or building, and therefore carry highly redundant information [5]. This motivated a third representation, object-based image analysis (OBIA), where land is described by meaningful spatial units rather than individual pixels [6]. Objects are a natural representation, since they are closer to how land is used and managed. Despite the idea being proposed a couple of decades ago [7], objects remain far less common than rasters because they are especially tricky to obtain from raw imagery. Classical OBIA constructs them by segmenting images, a process that can depend strongly on the sensor, season and hand-tuned parameters, and on the type of land [5, 8].

At the same time, the amount and diversity of data available about the landscape has grown drastically in recent years. Satellites now provide repeated observations of the Earth’s surface at global scale and across multiple spatial, spectral and temporal resolutions [9], which creates a continuous stream of information about the land. This, in turn, has enabled foundation models for Earth observation, which learn general-purpose representations from large remote-sensing datasets [10]. AlphaEarth Foundations [11], for example, summarizes each 10 m pixel of the land every year as a 64-dimensional embedding, a compact vector that integrates optical, radar, elevation and climate data.

Pan/Zoom

These learned representations raise a question that is largely hidden by the raster representation. What should the spatial unit of land be? With a fixed pixel size, the same spatial scale is imposed everywhere, despite the strong heterogeneity of the landscape. We know that some features require fine-grained representations, while others form coherent patterns over much larger areas [13]. A single resolution, however, ignores this variation, and while image segmentation can let object size vary, its settings have to be tuned by hand, by trial and error, for each sensor and landscape [14, 5]. We therefore ask whether meaningful spatial units can instead emerge directly from the learned embeddings, without having to define them manually.

Here we use Earth embeddings to cluster the land into adaptive land units, producing a compact graph of the landscape in which pixels that belong to the same object are grouped together. We use Bayesian hierarchical clustering [15], constrained so that only neighboring clusters can merge, and parallelize it to run on a GPU. By processing overlapping tiles with a ghost area around each tile [16], the clustering scales to an entire country with few seam artifacts.

The resulting representation adapts its spatial resolution to the landscape. Units become large where the land is uniform and small where it changes, rather than having a fixed imposed size (Fig. 1), which makes the representation compact enough for a whole country to fit on a single consumer-grade GPU. We show that, at the same number of cells and without using any labels, these units describe the land better than a regular grid or than units cut from the satellite image itself, and that they work as modeling units, both for mapping land use and for spatial models of a whole country that do not fit in memory on pixels.

2.1 Data and study area

2.2 Parallel Bayesian hierarchical clustering

Each land pixel is a node of the planar 4-neighbor pixel graph and carries its 64-dimensional embedding. We merge them into adaptive land units using Bayesian hierarchical clustering (BHC) [15], restricted so that only clusters sharing an edge of the graph can merge, that is, every merged unit is a connected region of the original pixel map. For each candidate merge, BHC compares the probability that the pixels of both clusters come from a single distribution with that of the separate clusters they already form, and merges when the posterior probability of the merge exceeds one half (\( r > 0.5\) ). More precisely, each unit is modeled as a Gaussian around its mean embedding, and a single parameter \( s\) , the within-unit variance relative to the national variance, sets how homogeneous a unit ends up being (Supplementary Section S1).

An important limitation of BHC is that it merges one pair at a time, which takes \( n - 1\) sequential steps for \( n\) pixels, that is around a million for a single 10 km tile. But because merges are restricted to graph neighbors, merges in different parts of the map do not interact, and we propose a parallel alternative, parallel BHC (PBHC), in which all clusters can merge at the same time. In each round, every cluster points to the neighbor with the most probable merge, provided \( r > 0.5\) , and groups of clusters linked by these pointers merge in a single step (Borůvka-style, Fig. 2a), using the BHC recursion extended to several children as in Bayesian rose trees [18] (Supplementary Algorithm S1). The rounds stop when no such merge remains, and a 10 km tile takes 8–16 rounds and 1–3 s on one GPU instead of a million sequential merges.

Pan/Zoom

2.3 Tiling

In this section, we first evaluate how well the units describe the land. A good unit holds a single kind of land, such as one field, one forest or one lake, rather than a mix of them. We measure this by the purity of a partition against a reference map, here Basemap05 [12] or the declared crop of every field. Each unit is given the class that covers most of its pixels, and the purity is the share of labeled pixels that fall in a unit of their own class. It is therefore also the highest accuracy that any map built on these units could reach when evaluated against those labels. Impurity \( \varepsilon\) is one minus purity and is zero when each unit contains a single pixel. Splitting units cannot increase impurity, so we compare partitions at the same unit density \( \rho\) , the number of units per km\( ^2\) of land. At a given density, lower impurity indicates better agreement with the reference labels.

First, we compare the units with other partition algorithms, from regular grids to the superpixel segmenters of object-based analysis (Section 3.1). We also cluster the pixels of the satellite image directly instead of the embeddings, to measure which is the new information the embeddings add (Section 3.2). Then, we evaluate how they perform as modeling units, first for mapping land use over Denmark from a few labeled points (Section 3.3), and then for estimating the share of cropland in every 10 km cell with its uncertainty, from a model of the whole country that does not fit in memory on pixels (Section 3.4).

3.1 Units that follow the landscape

Pan/Zoom

Table 1. Comparison between partitions (10 tiles of 20 km; squares of 180 m). Variance is the share of the embedding variance lost, \( \varepsilon_1\) and \( \varepsilon_2\) the impurity at Basemap05 levels 1 and 2, and \( \varepsilon_\mathrm{c}\) the crop impurity against the declared crop of every field. The reported time is per tile

Clustering	\( \rho\) (km\( ^{-2}\) )	Variance \( \downarrow\)	\( \varepsilon_1\)  \( \downarrow\)	\( \varepsilon_2\)  \( \downarrow\)	\( \varepsilon_\mathrm{c}\)  \( \downarrow\)	Time (s) \( \downarrow\)
PBHC (GPU)	32.4	0.20	0.132	0.175	0.094	9.6
SLIC [19]	25.1	0.21	0.148	0.193	–	94.5
Felzenszwalb [20]	22.0	0.48	0.252	0.307	–	27.9
Quadtree [21]	22.0	0.46	0.219	0.261	–	6.2
Squares	31.6	0.32	0.196	0.240	0.178	0.1
PBHC, Sentinel-2	36.7	–	0.175	0.234	0.133	5.4
Ward (Reference) [22]	32.4	0.16	0.128	0.170	–	1275

On 10 random tiles of 20 km, we compare PBHC, at the same number of units, with the partitions that could be used instead, namely grids (regular squares and an adaptive quadtree), superpixel segmenters as used in object-based analysis (SLIC [19] and Felzenszwalb’s graph segmentation [20]), and Ward’s agglomerative clustering [22], the classical alternative to a Bayesian criterion. PBHC cuts the land into 32 units per km\( ^2\) with an impurity of 0.13 (Fig. 3a). To be as pure, regular squares need 4.8 times as many units (2.3–6.4 across tiles). Over the whole country, the 1.27M units reach a purity of 0.881 on held-out regions, against 0.825 for the same number of squares (cells of 180 m). A quadtree, which splits squares wherever the embeddings vary, needs 4.5 times as many units. Adapting the size of a cell is therefore not enough, and its boundaries have to follow the land. The advantage holds against the declared crop of every field. Squares need 3.7 times as many cells to match the units, and unit boundaries are almost three times as likely as square boundaries to lie on a field boundary.

Ward’s agglomerative clustering on the same graph, the classical criterion for compact clusters, finds slightly purer partitions (Table 1). Nevertheless, it must be told how many units to make, has no parallel or tiled form, and takes 130 times longer per tile, while PBHC reaches nearly the same purity, decides its own number of units, and scales to the country. Note that Ward’s criterion is also, apart from a shape term, essentially the merge rule of the multiresolution segmentation in eCognition [14]. Of the superpixel segmenters, SLIC on the embeddings comes closest, but it takes ten times longer per tile and needs up to 100 GB of memory, which kept it from running above 65 units per km\( ^2\) . Felzenszwalb’s graph segmentation collapses at coarse scales.

3.2 The benefits of embeddings over satellite images

3.3 Mapping land use from sparse labels

3.4 Whole-country estimates with uncertainty

Pan/Zoom

Pan/Zoom

Other than for storage efficiency, a big incentive to create accurate yet compact representations that fits on one GPU is that they make joint models of the whole country possible. To illustrate that, we estimate the cropland share of every 10 km cell, alongside its uncertainty, from a survey of up to 10{,}000 points over the labeled land (Supplementary Section S5). Let each unit \( u\) have a probability \( p_u\) of being cropland, modeled as

\[ \begin{equation} \operatorname{logit} p_u = \mathbf{x}_u^\top \boldsymbol{\beta} + g_u , \end{equation} \]

(1)

where \( \boldsymbol{\beta}\) are the regression coefficients on the unit features \( \mathbf{x}_u\) , and the estimated share of a cell \( \hat \theta\) is the average of \( p_u\) over its pixels. While the regression alone gets most units roughly right, nearby units tend to be wrong in the same direction. The spatial field \( g_u\) corrects for this with a term that couples connected units in the graph. Because it links every unit to its neighbors, and through them to the rest of the country, it has to be fitted to all of Denmark at once, which means inverting a sparse matrix with one row per unit. On regular squares, even with the field, the estimates are less accurate and their intervals too narrow (Fig. 4a). On units without the field, the estimates are as accurate as with it but overconfident, and almost half of the 90 % intervals miss the true share (Fig. 4b). Finally, when accounting for the field in the units, the intervals are calibrated, covering the true share in 92 % of the cells (Fig. 4c). On pixels the model runs out of memory in most 60 km windows even on a 24 GB GPU. For window sizes that still fit in memory, it is up to 150 times slower, while on units, the entirety of Denmark takes 1.2 GB and 45 s to run (Fig. 5).

Adaptive land units compress the landscape while retaining its underlying structure. At the same number of cells they describe land use and crops better than a grid, lose only a few points of mapping accuracy against pixels, and bring a whole country onto one GPU. While we have shown this for two illustrative cases, its use is more general. Many land-use models are iterative, simulating change year after year, calibrating parameters or comparing scenarios, and each pass over a compact national graph takes seconds instead of a pass over hundreds of millions of pixels.

Much of this comes from clustering the embeddings rather than the image. Since the embeddings already describe what is on the ground, the clustering only has to find where it the changes are significant, and we see that compared to the raw satellite images, it needs twice less units to do so. The resulting graph neighbors also signal stronger information then when unclustered. On a grid, two cells are neighbors because of where the raster happens to fall, while two units are neighbors because two parts of the landscape meet there, which is why the same spatial field gives calibrated intervals on units but not on squares. Pixels remain the better choice when labels are plentiful. In most applications, however, much of the land is uniform over large areas, a field, a forest or a lake, and describing it pixel by pixel spends memory and computation where nothing changes. Units are therefore most useful when labels are sparse, models need to be fast, or the whole country has to be modeled at once.

Like any compression, the units trade some fine detail for compactness. They are also built one year at a time, and temporal applications may need a partition tuned for consistency over time (Supplementary Section S3). Finally, purity and accuracy measure agreement with Basemap05, which carries errors of its own.

We asked whether good land units can emerge from the data instead of being fixed in advance. Clustering Earth embeddings with a parallel Bayesian hierarchical clustering, we turned the 402M land pixels of Denmark into 1.27M adaptive land units without any labels. They describe land use and crops better than a grid with the same number of cells, and they make whole-country models with reliable uncertainty possible on a single GPU. As Earth embeddings now cover the planet every year, adaptive land units open a new way to represent land, in which the spatial unit of every model is learned from the data rather than imposed by a grid. Code and data are openly available [24].

This work was supported by a research grant (VIL78584) from Villum Fonden (project DK-Future, Villum Synergy Programme).

Our parallel version of Bayesian hierarchical clustering (BHC) stays close to the original formulation of Heller and Ghahramani [15]. It differs in two ways, as only clusters that share an edge of a planar graph can merge, and many merges happen at once instead of one pair at a time.

A.1 Algorithm

A.2 Model

To cover the country, we cluster 10 km tiles, each with a 250 m ghost area on every side [16], and keep only the core of each tile. Clusters wholly inside the core are kept as they are, but a cluster cut by the core edge may have been shaped by the tile border, and we dissolve its pixels within 250 m of that edge back into single pixels while the rest stays one piece. A final round of the same algorithm over the stitched pieces clusters these seam zones with context from both sides. The national build repeats the scheme one level up, on 60 km blocks of 36 tiles, each block with its own 250 m ghost area.

To measure what is left of the seams, we compare how often a pixel edge on a seam separates two units with how often any other edge does. On one 60 km block, independent tiles give a ratio of 12.3. Stitching without a ghost area gives 6.3, and the full scheme 1.9. Over the whole country, seam edges remain 1.4 to 2.0 times as likely as other edges to separate units, but the extra boundaries on seams make up about 0.04 % of all unit boundaries.

All comparisons use the same 10 random 20 km tiles as the main text, with PBHC at \( s = 1\) as the reference. To be as pure as PBHC at levels 1 and 2, squares need 4.8 and 4.2 times as many units (2.3–6.4 over tiles at level 1), the quadtree 4.5 and 3.9 times (3.9–5.1), and PBHC on the Sentinel-2 image 2.1 and 2.3 times (1.7–2.4).

Ward’s agglomerative clustering [22] runs on the same embeddings and 4-neighbor graph, on each connected component of land, with the tree cut at PBHC’s number of units. It needs 0.86 times as many units as PBHC for the same purity (0.44–1.01 over tiles), but takes 11–28 min and 11–16 GB per tile on the CPU, against about 10 s on the GPU, needs the number of units as input, and has no stitching scheme for a national build.

Against the declared crops of 2024 (14 crop groups, field pixels only), PBHC is purer than both squares and the units cut from the Sentinel-2 image from about 7 to 100 units per km\( ^2\) . At coarser densities its large units merge neighboring fields and squares become purer, and at finer densities the units cut from the Sentinel-2 image are purer. At the matched density, 72 % of the field area lies in fields that are at least 80 % covered by a single unit, against 3 % for squares, and 66 % of the unit boundary edges lie within one pixel of a field boundary, against 24 % for squares, about what chance alone gives.

We also built and scored the units separately for every year from 2017 to 2025, and their crop purity stays between 0.90 and 0.92. Between consecutive years, however, about half of the land lies in units with no match in the other year (intersection over union below 0.5). Units regroup no more often on fields whose crop changed than on fields that kept it (43 % against 48 % of field pixels), and the changes therefore do not follow crop rotation.

We map the 11 Basemap05 level-1 land classes of 2024 over all of Denmark from \( K\) labeled points, with a model on units, on as many squares, or on pixels (Table 2).

We divide Denmark into 10 km blocks and group them into 40 km groups, which we assign to training, validation and test in proportions of 60, 20 and 20 % of the land. Training points within 5 km of a validation or test group are dropped, and a unit or square belongs to the group of its majority pixel.

For each of 3 seeds, a budget of \( K\) points takes the first \( K\) of \( 10^6\) points drawn uniformly over the training land, and all models train on the same points, a unit or square model on the units or squares that hold at least one point, labeled by their majority class. The column All of Table 2 is a single run with \( 10^7\) points, which labels every unit.

Units and the 1.26M squares of 180 m, as many as there are units, are classified by a multilayer perceptron (MLP) on their mean embedding and the log of their area. Pixels are classified by an MLP on their own embedding or by a convolutional network (CNN) of three \( 3 \times 3\) layers, which sees a \( 7 \times 7\) pixel window. All models use hidden size 128, dropout 0.1, unweighted cross-entropy, and Adam with learning rate \( 10^{-3}\) and weight decay \( 10^{-4}\) , with early stopping on the validation loss. A two-layer graph neural network on the land-unit graph is as accurate as the MLP (0.825 against 0.826 at \( K = 10^5\) ). Accuracy is the share of test pixels predicted correctly, where a unit or square passes its prediction to all of its pixels.

Inference is timed for all of Denmark on one RTX 4090 GPU (24 GB). Units and squares fit on the GPU as one graph (0.47 and 0.43 GB), while the pixel models read the raster from a local int8 copy in 580 tiles, 38 GB in total.

We estimate the cropland share of every 10 km cell of Denmark, with a 90 % interval, from a survey of \( K\) points (Figs. 4 and 5 ). The cells are only used for reporting, and the model is fitted on the 1.27M adaptive land units, with survey points drawn uniformly over the labeled land for \( K = 10^2\) , \( 10^3\) and \( 10^4\) and 3 seeds. Each point records whether its pixel is Cropland in Basemap05 level 1, and we count for each unit \( u\) the \( m_u\) points it holds and the \( k_u\) of them that are cropland.

The probability \( p_u\) that a pixel of unit \( u\) is cropland is modeled as

and the survey points of each unit give a binomial likelihood, \( k_u \sim \mathrm{Bin}(m_u, p_u)\) . Here \( \mathbf{x}_u\) holds an intercept and the standardized mean embedding and log area of the unit, as in Section D, and \( g_u\) is a spatial field shared by neighboring units. The coefficients have the prior \( \boldsymbol{\beta} \sim \mathcal{N}(0, I)\) , except for a standard deviation of 10 on the intercept. The field follows a proper conditional autoregressive (CAR) prior on the land-unit graph,

where \( W_{uv}\) is the length of the border between units \( u\) and \( v\) , divided by its mean over all bordering pairs, and \( D\) is the diagonal matrix of the row sums of \( W\) . \( Q\) is nonzero only for the 3.5M pairs of units that share a border, and the term \( \delta (D + I)\) , with \( \delta = 10^{-3}\) , makes it invertible, where a smaller \( \delta\) lets the field vary over longer distances.

We choose \( \sigma_g\) from \( \{0.5, 1, 2, 4\}\) by fitting on four fifths of the survey points and scoring the log-likelihood of the remaining fifth, and then refit on all points. The posterior of \( \boldsymbol{\beta}\) and \( g\) is approximated by a Gaussian at its mode (Laplace approximation), which we find by Newton’s method. The Hessian of the negative log posterior is as sparse as \( Q\) , and each Newton step is solved by conjugate gradients, which only needs sums over the edges of the graph, and memory grows only with the number of units and edges. We never form the posterior covariance and instead draw 100 joint samples by perturb-and-MAP [26].

Units and cells do not line up, and each unit counts toward cell \( c\) by its number of labeled pixels in that cell, \( n_{uc}\) . For each posterior sample, the cropland share of the cell is \( \hat\theta_c = \sum_u n_{uc} \, p_u / \sum_u n_{uc}\) . The estimate is the mean of the 100 samples, and the 90 % interval runs from their 5th to their 95th percentile. We score them against the true share \( \theta_c\) , the share of labeled pixels of the cell that are cropland in Basemap05, on the 406 cells with at least 50 km\( ^2\) of labeled land. Coverage is the share of cells whose interval contains \( \theta_c\) , with a target of 0.90, and the mean absolute error (MAE) compares the estimates with \( \theta_c\) .

Table 3 compares the model with three baselines, starting with the same model on units without the field (\( g = 0\) ). On squares, the field links 180 m squares that share a side, with an independent effect per square. The survey alone estimates each cell by the share of its points that are cropland, with a Wilson 90 % interval, on the cells that hold at least one point.

E.1 Memory and time