Learning Geospatial Region Embedding with Heterogeneous Graph

To effectively encode spatial, environmental, and societal signals into intra-region representations, we extract multi-modal geo-contextual data including satellite imagery, POI information, and the location of the region of interest.

Spatial Position Embedding: We adopt a simple yet effective method to extract the spatial information of the region of interest, drawing inspiration from industrial conventions [42, 9]. Specifically, we divide the overall geospace into multiple $1\textrm{km}\times 1\textrm{km}$ spatial grids with coordinate system origin ($lon^{0},lat^{0}$) and use the abscissa $x_{R}$ and ordinate $y_{R}$ of the region of interest $R$ as spatial information $E_{pos}(x,y)$, illustrated in Figure 3. Notably, $E_{pos}$ is directly related to the geo-coordinate ($lon,lat$) at the centroid of the targeted region, for any region $E^{i}$ in the target space, we have:

$$E^{i}_{pos}=(lon^{i}-lon^{0},lat^{i}-lat^{0})\Delta*D^{-1}$$

(1)

where $\Delta$ is the transform vector decided by the coordinate system. $D$ denotes the scale of the grids, and for $1\textrm{km}\times 1\textrm{km}$ grids, $D=1$.

Environmental Feature Embedding: We utilize satellite imagery to mine the environmental features and expect the satellite imagery encoding to be efficient and concise. To achieve this, we draw inspiration from the European Space Agency (ESA) with ESA WorldCover Dataset [50], elaborate in Appendix A, and design a segmentation-based process to encode the satellite images instead of the most commonly used visual encoders such as CNNs. Given a satellite image $S_{i}$ of a region of interest $R_{i}$, we conduct semantic segmentation to obtain a series of geo-entities $\{entity_{1},entity_{2},...,entity_{j}\}_{j=1}^{J}$ within that region, such as developed areas, grass, and trees based on ESA Worldcover, leveraging their area proportion as environmental feature embedding as below:

$$E_{env}=\{p_{ent_{1}},p_{ent_{2}},...,p_{ent_{j}}\}_{j=1}^{J}$$

(2)

where $p_{ent_{j}}=\frac{A_{j}}{A_{R}}$, and $A_{j}$ is the area occupied by $entity_{j}$, $A_{R}$ is the area of region $R_{i}$, $J$ is the number of environmental entity types. We discuss the motivation for using the segmentation-based approach instead of visual encoders and additional efficiency experiments in Appendix B.

Societal Feature Embedding: To effectively represent the societal feature of a region, given POI dataset $D_{POI}$ and the region of interest $R_{i}$, we firstly find all points $D^{R_{i}}_{POI}$($poi_{1},poi_{2},...,poi_{j})$ which are located at $R_{i}$. Then, we count the proportion of different $K$ POI categories in the region $R_{i}$ to get $E_{soc}^{\prime}$ as below:

$$E_{soc}^{\prime}=\{p_{poi_{1}},p_{poi_{2}},...,p_{poi_{k}}\}_{k=1}^{K},\quad p_{\text{poi}_{k}}=|C_{k}\cap D_{POI}^{R_{i}}|/|D_{POI}^{R_{i}}|$$

(3)

where $K$ is the number of POI categories, $C_{k}=\{poi_{j}\in D_{POI}\mid\text{category}(poi_{j})=k\}$. However, directly employing the proportion of different POI categories $E_{soc}^{\prime}$ as societal features can still be problematic, since it cannot distinguish between areas with high POI density (usually with a higher degree of socialization) and those with relatively low POI density. Therefore, we transform $E_{soc}^{\prime}$ into final societal feature embedding $E_{soc}$ by multiply a social impact factor $f$ as below:

$$E_{soc}=f\cdot E_{soc}^{\prime},\quad f=log(D^{R_{i}}_{POI}+1)$$

(4)

Finally, the intra-region feature embedding is generated by concatenating the spatial position embedding $E_{pos}$, environmental feature embedding $E_{env}$ and societal feature embedding $E_{soc}$:

$$E_{intra}=concat[E_{pos},E_{env},E_{soc}]$$

(5)

Our objective is to model pairwise second-order relationships between regions and their adjacent neighbors and higher-order dependencies with distant region groups exhibiting similar environmental or societal characteristics. Inspired by hypergraph theories [21, 4], elaborated in Appendix C.1, we employ a heterogeneous graph formulation GeoHG to represent these relationships explicitly. As shown in Figure 4, we take heterogeneous nodes in the graph structure as transfer nodes to achieve high-order message passing in the hypergraph. Given a pair of $1\textrm{km}\times 1\textrm{km}$ satellite image $S_{i}$ and POI data $P_{i}$ in $D_{satellite\_{POI}}$, we construct an undirected weighted heterogeneous graph $G_{i}$, containing three types of nodes: regional nodes, environmental entity nodes, and societal entity nodes. Specifically, regional nodes correspond to $1\textrm{km}\times 1\textrm{km}$ grid cells from the geospace; environmental entity nodes represent the $J$ distinct geo-entity classes detected from satellite image $S_{i}$, which are consistent with intra-region environmental features in Section 3.1; societal entity nodes constitute the $K$ different POI categories. The overall structure of GeoHG is shown in Figure 2 and the details are discussed in Appendix 4. We then elaborate on how second-order and high-order information is obtained through GeoHG $\mathcal{G}_{i}=(\mathcal{V},\mathcal{E})$, respectively.

Second-Order Relation Representation: Second-order relations capture the pairwise spatial adjacency between regional nodes (grid cells) in $\mathcal{G}_{i}$. We construct undirected edges $\mathcal{E}_{RNR}$, Region Nearby Region (RNR) between regional nodes whose corresponding grid cells are spatially adjacent in a $3\times 3$ grid, thereby encoding these local second-order dependencies.

$$\mathcal{E}_{RNR,mn}=1,\text{if }\mathcal{V}_{m}\text{ and }\mathcal{V}_{n}\text{ are {adjacent} in geospace}$$

(6)

High-Order Relation Representation: Differ from traditional multi-view graphs [7, 59] represent complex correlations through additional second-order region connections, we explicitly model higher-order relations by leveraging the environmental entity nodes $\mathcal{V}_{env}$ and societal entity nodes $\mathcal{V}_{soc}$ in $\mathcal{G}_{i}$. Specifically, we construct heterogeneous weighted edges $\mathcal{E}_{ELR}$ and $\mathcal{E}_{SLR}$ through $\mathcal{V}_{env}$ and $\mathcal{V}_{soc}$ to associate each regional node with its constituting environmental entities (e.g., water, vegetation) and societal entities (e.g., educational, commercial POIs). For $p_{ent_{it}}\geq\theta_{Env}$ or $f\cdot p_{poi_{ik}}\geq\theta_{Soc}$:

$$\mathcal{E}_{ELR,ij}=1\cdot p_{ent_{ij}},\quad\mathcal{E}_{SLR,ik}=1\cdot f\cdot p_{poi_{ik}}$$

(7)

where $p_{ent_{ij}}$ represents the area proportion of geo-entity $j$ at region $i$, $p_{poi_{ij}}$ denotes the quantity proportion of POI category $k$ at region $i$ and $f$ is the social impact factor elaborated in Equation 4. $\theta_{Env}$ and $\theta_{Soc}$ serve as hyperparameters designed to optimize the graph structure for conciseness and efficiency. This edge formulation allows encoding higher-order relationships between regions exhibiting similar environmental/societal attributes, even if they lack spatial proximity. By combining the second-order spatial adjacency edges with these higher-order hyperedge associations, $\mathcal{G}_{i}$ holistically represents the mixed-order relational patterns within the region.

Having derived intra-region features capturing regions’ intrinsic attributes and inter-regional dependencies encoded in the heterogeneous graph $\mathcal{G}_{i}$, we employ a model-agnostic graph neural network framework to jointly reason over the intra-regional and inter-regional information. In the graph, each regional node $v$ in $\mathcal{G}_{i}$ is represented with an initial node feature $\mathbf{x}_{v}$ corresponding to its intra-regional feature $E_{intra}$ elaborated in Sec 3.1. Note that for simplicity, instead of utilizing hypergraph neural networks [4, 44, 39], we model hyperedges with the form of complete subgraphs in $\mathcal{G}_{i}$. Therefore, we apply Heterogeneous Graph Neural Networks (HGNNs) over $\mathcal{G}_{i}$ to update the node representations by iteratively aggregating and transforming multi-hop neighborhood information:

$\displaystyle\mathbf{x}_{v}^{(l+1)}$

$\displaystyle=\text{HGNN}^{(l)}\left(\mathbf{x}_{v}^{(l)},\mathcal{N}(v),\mathcal{G}_{i}\right),\mathbf{x}_{v}^{(l)}\in\mathbb{R}^{d}$

(8)

where $\mathcal{N}(v)$ denotes the neighbors (spatial and higher-order) of node $v$ in $\mathcal{G}_{i}$. After L layers of updates, the final node embedding $\mathbf{x}_{v}^{(L)}$ integrates cues from the intra-regional features as well as mixed-order inter-region dependencies. This final embedding $\mathbf{x}_{v}^{(L)}$ serves as the comprehensive representation integrating intra-regional and inter-regional information for the given region node. Our framework is model-agnostic, enabling flexible integration of any HGNN variant for neighborhood aggregation, and allowing easy extension to more expressive HGNN models.

To effectively empower various geospatial tasks with the devised representations, we present two training strategies: self-supervised pre-training that learns generalizable region representations without task-specific labels, enabling efficient transfer to diverse downstream tasks; and end-to-end training that directly optimizes task-specific objectives for peak performance.

Self-Supervised Pre-training: Inspired by CLIP [29], we pre-train our model using a contrastive learning paradigm where we maximize the similarity between a region $\mathbf{R}_{i}$ and its corresponding positive region groups $\mathbf{C}_{i}$, while minimizing its similarity with other region groups in the batch. Specifically, we sample regional nodes adjacent to $\mathbf{R}_{i}$ as well as nodes exhibiting similar intra-regional features to construct $\mathbf{C}_{i}$. During training, we leverage a GNN model $\mathcal{M}_{pretrain}$ to obtain embeddings for $\mathbf{R}_{i}$ and all nodes in $\mathbf{C}_{i}$. The embeddings of the sampled positive regional nodes are further pooled to obtain the unified positive region representation $\mathbf{e}_{j}$. The process of generating embeddings for $\mathbf{R}_{i}$ and $\mathbf{C}_{i}$ is conducted as follows:

$$\mathbf{e}_{i}=\mathcal{M}_{pretrain}(\mathbf{X}_{i}),\quad\mathbf{e}_{j}=\operatorname{Pool}\left(\mathcal{M}_{pretrain}(\mathbf{X}_{j}),{j\in\mathcal{C}_{i}}\right)$$

(9)

We further define a similarity function $f_{\mathbf{score}}$ to measure the similarity between the representations $\mathbf{e}_{i},\mathbf{e}_{j}$. $f_{\mathbf{score}}$ can be a simple dot product or a more complex metric and then employ an InfoNCE-based loss [28] to conduct contrastive learning:

$$\mathcal{L}_{pretrain}=\mathbb{E}[-\log\frac{\exp\left(f_{\mathbf{score}}\left(\mathbf{e}_{i},\mathbf{e}_{j}\right)\right)}{\sum_{\forall i,n}\exp\left(f_{\mathbf{score}}\left(\mathbf{e}_{i},\mathbf{e}_{n}\right)\right)}]$$

(10)

where $f_{\mathbf{score}}(\mathbf{e}_{i},\mathbf{e}_{j})$ represents the score of positive pairs while $f_{\mathbf{score}}\left(\mathbf{e}_{i},\mathbf{e}_{n}\right)$ refers to the scores of negative pairs. After obtaining the pre-trained model $\mathcal{M}_{pretrain}$ through self-supervised learning, we fine-tune it for different downstream tasks by adding three trainable linear layers. The weights of these additional linear layers are updated during training for the specific downstream task while the weights of $\mathcal{M}_{pretrain}$ itself remain fixed.

End-to-End Training: For a specific downstream task, we directly optimize the whole GeoHG’s parameters from scratch using the task’s supervised signal. Taking the regional regression task as an illustrative example, we adopt the HGNN model $\mathcal{M}$ to obtain embeddings for given region $i$, then feed it into a three-linear layer regression head. $\mathcal{M}$’s parameters are updated by minimizing the Mean Square Error (MSE) loss on the training data, jointly learning task-specific region representations.

We adopt Mean Absolute Error (MAE), rooted mean squared error (RMSE), and coefficient of determination ($R^{2}$) as evaluation metrics. An increase in $R^{2}$, along with a reduction in MAE and RMSE values is indicative of enhanced model accuracy. We report the performances of each model in Table 1 and the table with RMSE metrics is provided in Appendix F. From these tables, we have three key findings: 1) Both GeoHG and GeoHG-SSL outperform all competing baselines over the 5 datasets for 4 cities. For instance, GeoHG surpasses the previous SOTA performance in Beijing, achieving $R_{2}$ improvements of +16.7%, +14.9%, +6.1%, +37%, +33% in Carbon, Population, GDP, Night Light and $PM_{2.5}$ tasks respectively. Consistent trends are observed across other cities, underscoring GeoHG’s stable accuracy, versatility, and strong generalization capabilities for geospatial region embedding. 2) The end-to-end trained GeoHG outperforms its pre-trained version GeoHG-SSL in most tasks while GeoHG-SSL shows a good performance in multi-task transferring. 3) Multisource and multimodal approaches, i.e., UrbanVLP [9] and GeoStructual [20], largely surpass traditional satellite imagery-based models by their more comprehensive embedding views.

Effects of Core Components. To examine the effectiveness of each core component in our proposed framework, we conducted an ablation study based on the following variants for comparison: a) w/o Env, which excludes environment information from the satellite imagery for embedding. b) w/o Soc, which omits POI entities for embedding. c) w/o Pos, which does not utilize the corresponding location information of the regions. The $R^{2}$ results for five tasks in two cities, Guangzhou and Shanghai, are displayed in Figure 5. We can observe that removing location information markedly degrades performance across all tasks. Meanwhile, environmental features severely impact model accuracy on the PM2.5, Night Light, and GDP tasks. In contrast, excluding POI data slightly reduces performance on the Carbon task. Appendix F.3 shows ablation results for the other two cities.

Effect of Graph Construction. Our framework utilizes heterogeneous graph and hyperedges between regions to reflect the geospatial high-order relations. To validate the effectiveness of high-order relation representation, we compare GeoHG against two variants: GeoHG-MLP which discards graph structure and relies only on intra-region feature representation, employing a 3-layer MLP for regression; GeoHG-Mono which keeps edges in adjacent regions while discarding hyperedges for high-order relations. The results presented in Figure 5 indicated that our GeoHG significantly benefits from efficient representation of complex mixed-order relations in geospace and integration of intra-region and inter-region representations, thereby resulting in enhanced performance.

Interpretation of Mixed-order Geospatial Correlations. To investigate the power of heterogeneous graph structure in capturing mixed-order geospatial relations, inspired by GNNExplainer [46], we depict the trend of learned geospatial dependency among Carbon Emission datasets in Guangzhou. Notably, we do not incorporate any external features (e.g., road network, human mobility data) to construct aided edges in the graph. By selecting the top $N$ important nodes by their weight in the GNN prediction process, as illustrated in Figure 6, we observe that GeoHG effectively captures the high-order correlations between the target region and remote region groups through message passing within the heterogeneous graph structure. The details about GNNExplainer and experiment results are introduced in Appendix F.5.

Applying geospatial models globally across tasks often requires large labeled datasets for supervised training - a process that is time-consuming, computationally expensive, and hindered by data scarcity. We therefore evaluate GeoHG’s performance under limited data regimes of 5%, 10%, and 20% available training samples, with results shown in Figure 7. Results demonstrate GeoHG’s strong data efficiency - with only 5% data, GeoHG outperforms previous SOTA methods like UrbanVLP [9] and StructuralGeo [20] trained on the entire training dataset across all tasks. Utilizing only 20% of the data, GeoHG suffers minor performance degradation of 4%, 0.5%, 2.01%, 1.7%, and 3.2% on the Carbon, GDP, Light, PM2.5, and Population tasks, respectively. These findings illustrate GeoHG’s promising potential for data-efficient global deployment across diverse geospatial prediction tasks.

Geospatial data interpolation methods like Inverse Distance Weighting (IDW) [58, 24] and Universal Kriging (UK) are employed to generate predictions for data-scarce regions. We evaluate GeoHG against these methods on population prediction in an enlarged area of Shenzhen, with only 5% visible data (863 points), tasked with inferring the remaining 16,401 regions (16,401 $km^{2}$), as illustrated in Figure 7. GeoHG accurately captures the true distribution, substantially outperforming IDW and UK which deviate significantly. This highlights that traditional methods, solely relying on spatial relationships for interpolation, fail to model the intricate socio-environmental dependencies critical for characterizing population distributions. In contrast, GeoHG’s effective geospatial representations and modeling of high-order relationships enable accurate data pattern learning from limited samples.