Areal Embedding: Methods & Applications

• Areal embedding is a technique that maps geographic regions, such as polygons and tiles, into a metric space to capture spatial dependencies and latent structures.
• It employs methods like set-indexed Gaussian processes, region-based integration, and graph-based vector representations to replace traditional adjacency matrices.
• Applications span environmental monitoring, real estate appraisal, and transportation demand prediction, with empirical results showing improved predictive performance and model coherence.

Searching arXiv for recent and directly relevant papers on areal embedding and adjacent formulations. Areal embedding denotes a representation of a spatial support—such as a polygon, grid cell, tile, patch, or region—in a space where dependence, similarity, or latent structure can be modeled directly. In the literature, the term appears in several technically distinct forms: set-indexing of Gaussian processes for areal and point data, region-wise integrals of latent Gaussian processes, pooled geospatial foundation-model features for polygons and tiles, graph- or coordinate-based vectors for downstream prediction, and scalar-field constructions over 2D embedded canvases (Godoy et al., 2022, Tanaka et al., 2019, Fang et al., 19 Jan 2026, Cheng et al., 2016). This suggests that areal embedding is not a single formalism but a family of constructions whose common aim is to make areas first-class representational objects.

1. Conceptual scope and terminology

Representative usages of the term differ by research tradition.

Literature	Embedded object	Representative formulation
Spatial statistics	Polygons and points as supports	Set-indexed GP via Hausdorff distance; region-integrated latent GP
Geospatial foundation models	Tiles, patches, pixels, polygons	Pooling or storing fixed-length vectors for areas
Predictive modeling	Coordinates, grid cells, postal areas	Learned vectors from contrastive, graph, or road-network models
Visualization	2D embedded point clouds	Continuous scalar field over an embedded canvas

In spatial statistics, areal embedding is explicitly defined as representing spatial supports such as polygons directly in a metric space so that a single spatial process can model both areal and point data. The Hausdorff–Gaussian Process embeds supports in the space of non-empty, bounded, compact subsets of $D \subset \mathbb{R}^d$, denoted $B(D)$, equipped with the Hausdorff metric (Godoy et al., 2022). In Spatially Aggregated Gaussian Processes, an areal embedding is a vector representation of a region constructed by integrating latent Gaussian processes or task-specific mixed outputs over the region (Tanaka et al., 2019).

In geospatial foundation-model work, areal embeddings are aggregated representations of geographic areas—polygons, grid cells, or tiles—constructed from underlying embeddings, or they are directly released as embedding products aligned to chip or pixel footprints (Fang et al., 19 Jan 2026, Czerkawski et al., 2024). In application-specific predictive models, the term covers coordinate-to-vector maps that can be averaged over polygons, node states in multi-feature spatio-temporal graph models, and road-network-informed vectors attached to grid cells for house price interpolation (Holvoet et al., 22 Nov 2025, Han et al., 2024, Han et al., 2023). In visualization, the term is used differently: a 2D embedding of high-dimensional data is treated as a spatial domain over which a continuous scalar field is estimated (Cheng et al., 2016).

2. Set-indexed areal embedding in spatial statistics

The Hausdorff–Gaussian Process (HGP) operationalizes areal embedding by treating each support $S$ as an element of $B(D)$ and using the Hausdorff distance as the metric input to a covariance kernel. For sets $A,B \subset \mathbb{R}^d$, the Hausdorff distance is

$$d_H(A,B) = \max\left\{ \sup_{a\in A} \inf_{b\in B} ||a-b||,\ \sup_{b\in B} \inf_{a\in A} ||a-b|| \right\}.$$

If $A$ and $B$ are singletons, then $d_H(A,B)=||a-b||$, so the construction reduces immediately to ordinary point-referenced Gaussian process modeling. This yields a unified formulation for point and areal supports (Godoy et al., 2022).

For supports $S_i \in B(D)$, HGP defines a Gaussian process $B(D)$0 with covariance

$B(D)$1

where $B(D)$2 is a marginal standard deviation function and $B(D)$3 is an isotropic correlation function on Hausdorff distance. The recommended default is the powered exponential correlation

$B(D)$4

with $B(D)$5 and $B(D)$6 reparameterized as the distance at which correlation drops to $B(D)$7. A useful heteroscedastic construction is

$B(D)$8

with $B(D)$9 allowing different marginal variances for point and polygon supports. Embedded into a GLMM, the latent effect enters as

$S$0

This replaces adjacency matrices entirely; spatial dependence is determined by Hausdorff distances rather than by a binary neighborhood graph (Godoy et al., 2022).

The construction is motivated by two limitations of conventional areal models. First, CAR, ICAR, MCAR, and DAGAR encode dependence through adjacency and therefore do not distinguish polygons by size, shape, or orientation. Second, aggregated-GP data fusion models define polygon-level effects by integrals such as

$S$1

which require numerical integration over polygons. HGP avoids such integrals by defining the latent variables directly on supports. Inference in the paper uses Bayesian posterior sampling with NUTS in Stan, with convergence assessed by split-$S$2, and exact GP inference retains the usual $S$3 scaling of dense covariance models (Godoy et al., 2022).

Empirically, HGP was competitive with specialized areal and fusion models. In areal simulations under aggregated-GP truth, HGP outperformed DAGAR strongly: $S$4 in $S$5 of replicates, while DAGAR was favored with $S$6 in only $S$7. Under BYM truth, HGP remained competitive, with $S$8 within $S$9 for $B(D)$0 of replicates. In the respiratory disease hospitalization application over 134 intermediate zones, HGP had the best LOOIC, $B(D)$1, versus $B(D)$2 for DAGAR and $B(D)$3 for BYM; its posterior for $B(D)$4 implied that correlation decays to $B(D)$5 by $B(D)$6 km. In the PM2.5 fusion application with points plus 184 tiles, HGP achieved the lowest RMSP, $B(D)$7, nominal CPP of $B(D)$8, and the best interval score, while also estimating larger point-level SD than areal SD ($B(D)$9, $A,B \subset \mathbb{R}^d$0) (Godoy et al., 2022).

3. Aggregation-based latent embeddings of regions

A second major formulation treats an areal unit as the integral of a continuous latent field over a region. In Spatially Aggregated Gaussian Processes (SAGP), independent latent GPs $A,B \subset \mathbb{R}^d$1 are linearly mixed into output functions

$A,B \subset \mathbb{R}^d$2

or, in vector form, $A,B \subset \mathbb{R}^d$3. An areal observation for output $A,B \subset \mathbb{R}^d$4 and region $A,B \subset \mathbb{R}^d$5 is

$A,B \subset \mathbb{R}^d$6

with the uniform-weight special case

$A,B \subset \mathbb{R}^d$7

The covariance between two aggregated observations is the double integral of the mixed kernel over both regions (Tanaka et al., 2019).

Within this framework, two areal embeddings are explicit. The latent-space areal embedding is

$A,B \subset \mathbb{R}^d$8

and the task-specific embedding for output $A,B \subset \mathbb{R}^d$9 is

$$d_H(A,B) = \max\left\{ \sup_{a\in A} \inf_{b\in B} ||a-b||,\ \sup_{b\in B} \inf_{a\in A} ||a-b|| \right\}.$$0

These constructions place regions from different outputs or domains in a common latent space and allow transfer learning because latent GPs are shared across datasets. Posterior means and variances of both point-level functions and region-integrated embeddings follow from standard GP conditioning, with cross-covariance terms expressed as integrals of kernels against region weights (Tanaka et al., 2019).

The SAGP viewpoint differs from the HGP viewpoint in where the embedding occurs. HGP embeds the support itself in a metric space of sets. SAGP keeps a continuous spatial field on the underlying domain and obtains areal representations by integration. This leads to different computational burdens. For arbitrary polygons, exact closed forms for

$$d_H(A,B) = \max\left\{ \sup_{a\in A} \inf_{b\in B} ||a-b||,\ \sup_{b\in B} \inf_{a\in A} ||a-b|| \right\}.$$1

are generally unavailable for common kernels such as squared exponential and Matérn, so the paper uses rasterization or quadrature, with complexity reduced by pooling kernel values by distinct inter-grid distances under stationarity (Tanaka et al., 2019).

The reported results show that these integrated embeddings can refine coarse-grained areal data effectively. On NYC and Chicago datasets spanning partitions of varying granularities, SAGP outperformed GPR, a two-stage GP, and SLFM. For NYC poverty, PM2.5, and crime, and Chicago poverty, SAGP achieved MAPE $$d_H(A,B) = \max\left\{ \sup_{a\in A} \inf_{b\in B} ||a-b||,\ \sup_{b\in B} \inf_{a\in A} ||a-b|| \right\}.$$2, $$d_H(A,B) = \max\left\{ \sup_{a\in A} \inf_{b\in B} ||a-b||,\ \sup_{b\in B} \inf_{a\in A} ||a-b|| \right\}.$$3, $$d_H(A,B) = \max\left\{ \sup_{a\in A} \inf_{b\in B} ||a-b||,\ \sup_{b\in B} \inf_{a\in A} ||a-b|| \right\}.$$4, and $$d_H(A,B) = \max\left\{ \sup_{a\in A} \inf_{b\in B} ||a-b||,\ \sup_{b\in B} \inf_{a\in A} ||a-b|| \right\}.$$5, respectively. Transfer learning across NYC and Chicago further improved Chicago poverty refinement, to approximately $$d_H(A,B) = \max\left\{ \sup_{a\in A} \inf_{b\in B} ||a-b||,\ \sup_{b\in B} \inf_{a\in A} ||a-b|| \right\}.$$6 versus approximately $$d_H(A,B) = \max\left\{ \sup_{a\in A} \inf_{b\in B} ||a-b||,\ \sup_{b\in B} \inf_{a\in A} ||a-b|| \right\}.$$7 for SLFM transfer (Tanaka et al., 2019).

4. Geospatial foundation models and Earth as an embedding product

In geospatial foundation-model work, areal embeddings are usually fixed-length vectors attached to spatial footprints and aggregated over larger areal units as needed. A formal definition is given by the embedding function

$$d_H(A,B) = \max\left\{ \sup_{a\in A} \inf_{b\in B} ||a-b||,\ \sup_{b\in B} \inf_{a\in A} ||a-b|| \right\}.$$8

or, with time,

$$d_H(A,B) = \max\left\{ \sup_{a\in A} \inf_{b\in B} ||a-b||,\ \sup_{b\in B} \inf_{a\in A} ||a-b|| \right\}.$$9

For an areal unit $A$0, the aggregated embedding can be written as continuous mean pooling,

$A$1

discrete averaging,

$A$2

or weighted pooling,

$A$3

Patch-level and pixel-level products are distinguished, and the ecosystem is organized into a three-layer taxonomy of Data, Tools, and Value (Fang et al., 19 Jan 2026).

This product-oriented perspective is developed further in the AlphaEarth and Major TOM lines of work. AlphaEarth produces a 64-dimensional vector per location and year that summarizes annual multi-modal Earth observation signals into a physically structured feature space. Using approximately 12.1 million location-year samples across the Continental United States from 2017 to 2023, the study reports that 12 of 26 environmental variables exceed $A$4, while temperature and elevation approach $A$5. The strongest dimension-variable relationships include A57 with annual precipitation ($A$6), A40 with daytime LST ($A$7), A48 with EVI ($A$8), A26 with tree cover ($A$9), A00 with evapotranspiration ($B$0), and A50 with mean air temperature ($B$1). The relationships remain robust under spatial block cross-validation with mean $B$2 and temporally stable with mean inter-year correlation $B$3 (Rahman, 10 Feb 2026).

Major TOM, by contrast, emphasizes global and dense release of embedding products. It extends Major TOM Core datasets by standardizing fragmentation, preprocessing, and packaging, and releases four dense embedding datasets—SSL4EO-S2, SSL4EO-S1, SigLIP on Sentinel-2 RGB, and DINOv2 on Sentinel-2 RGB—aligned to the Major TOM grid. These datasets were computed from more than 62 TB of raw imagery, distilling roughly 9.368 trillion pixels into 169+ million embedding vectors. Embeddings and metadata are stored in GeoParquet with columns including unique_id, embedding, grid_cell, timestamp, geometry, utm_footprint, utm_crs, pixel_bbox, centre_lat, and centre_lon, making them directly searchable and spatially indexable (Czerkawski et al., 2024).

Together, these works frame areal embeddings as first-class geospatial datasets rather than only as internal model features. The standardized access work extends TorchGeo with dataset classes, spatial and temporal slicing, dataset intersections, and samplers such as GridGeoSampler, so that areal embeddings can be loaded, aligned, and pooled over polygons in a common interface (Fang et al., 19 Jan 2026).

5. Learned areal embeddings for downstream prediction

Several recent models learn areal embeddings explicitly for prediction tasks. In multi-view contrastive risk modeling, the core object is a coordinate-to-embedding map

$B$4

trained by aligning a coordinate encoder with fused satellite-imagery and OSM views. The coordinate encoder uses spherical harmonics with $B$5 followed by a 2-layer SIREN MLP with 128 hidden units per layer. Polygon-level embeddings are obtained either by centroid evaluation, $B$6, or by spatial averaging,

$B$7

The model is trained with a symmetric InfoNCE loss using cosine similarity and $B$8. After training, inference requires coordinates only and is reported at $B$9 ms per coordinate. On 200,000 French real estate transactions, replacing raw latitude and longitude with these embeddings improved out-of-sample MSE across GLM, GAM, and GBM; the best reported variant, EU64_GS64, achieved test MSE $d_H(A,B)=||a-b||$0, $d_H(A,B)=||a-b||$1, and $d_H(A,B)=||a-b||$2, versus $d_H(A,B)=||a-b||$3, $d_H(A,B)=||a-b||$4, and $d_H(A,B)=||a-b||$5 for raw lat–lon (Holvoet et al., 22 Nov 2025).

In transportation demand prediction, areal embeddings arise as hidden node states in the Multi-Feature aware GCGRU. Each areal feature type is represented by a matrix $d_H(A,B)=||a-b||$6, projected by learned encoders to similarity scores, and normalized with a sentinel attention mechanism

$d_H(A,B)=||a-b||$7

The resulting graphs, together with the identity and proximity graphs, drive graph-convolutional GRU updates whose hidden states $d_H(A,B)=||a-b||$8 are the learned areal embeddings. On BusDJ and TaxiBJ, the ST-MFGCRN model outperformed the state-of-the-art baselines by up to $d_H(A,B)=||a-b||$9 on BusDJ and $S_i \in B(D)$0 on TaxiBJ; reported RMSE and MAE were $S_i \in B(D)$1 and $S_i \in B(D)$2 on BusDJ, and $S_i \in B(D)$3 and $S_i \in B(D)$4 on TaxiBJ (Han et al., 2024).

In real estate appraisal, areal embeddings are attached to grid cells built from a road-network graph. The study partitions the region into $S_i \in B(D)$5 cells and constructs a weighted adjacency

$S_i \in B(D)$6

where $S_i \in B(D)$7 is the set of grid cells traversed by road $S_i \in B(D)$8. Node2Vec then learns a $S_i \in B(D)$9 embedding $B(D)$00 for each areal unit by optimizing the standard skip-gram-with-negative-sampling objective over random-walk contexts. A simpler alternative is a fixed 2D sinusoidal positional encoding over grid indices. In the AMMASI model, the areal embedding of the house’s cell, $B(D)$01, is concatenated with house features and outputs of masked geographic and similar-feature attention:

$B(D)$02

Across the best HA versus HA+P choices per dataset for ASI, AMMASI achieved an average MAPE reduction of about $B(D)$03. The paper also reports that Node2Vec is not uniformly superior to sinusoidal 2D positional encoding, with region-specific variation in which representation performs best (Han et al., 2023).

These application-specific formulations share a practical pattern: areal embedding is used to compress spatial context into a vector that can be consumed by standard predictive architectures. The underlying constructions, however, differ substantially—contrastive coordinate encoders, learned graph states, and road