Positional Encoder Graph Neural Networks for Geographic Data

Graph neural networks (GNNs) provide a powerful and scalable solution for modeling continuous spatial data. However, in the absence of further context on the geometric structure of the data, they often rely on Euclidean distances to construct the input graphs. This assumption can be improbable in many real-world settings, where the spatial structure is more complex and explicitly non-Euclidean (e.g., road networks). In this paper, we propose PE-GNN, a new framework that incorporates spatial context and correlation explicitly into the models. Building on recent advances in geospatial auxiliary task learning and semantic spatial embeddings, our proposed method (1) learns a context-aware vector encoding of the geographic coordinates and (2) predicts spatial autocorrelation in the data in parallel with the main task. On spatial regression tasks, we show the effectiveness of our approach, improving performance over different state-of-the-art GNN approaches. We also test our approach for spatial interpolation, i.e., spatial regression without node features, a task that GNNs are currently not competitive at. We observe that our approach not only vastly improves over the GNN baselines, but can match Gaussian processes, the most commonly utilized method for spatial interpolation problems. The code for this study can be accessed via: https://github.com/konstantinklemmer/pe-gnn


INTRODUCTION
Geographic data describes all data that can be mapped onto the earth, for example using latitude and longitude coordinates. As such, geographic data is characterized by a natural geometric structure, which often strongly affects the observed spatial pattern. While traditional neural network approaches do not have an intuition to account for spatial dynamics, graph neural networks (GNNs) can represent spatial structures graphically. The recent years have seen many applications leveraging GNNs for modeling tasks in the geographic domain, such as inferring properties of a point-ofinterest [35] or predicting the speed of traffic at a certain location [5].
Nonetheless, as we show in this study, GNNs are not necessarily sufficient for modeling complex spatial effects: spatial context can be different at each location, which may be reflected in the relationship with its spatial neighborhood. The study of spatial context and dependencies has attracted increasing attention in the machine learning community, with studies on spatial context embeddings [24,34] and spatially explicit auxiliary task learning [22].
In this study, we seek to merge these streams of research. We propose the positional encoder graph neural network (PE-GNN), a flexible approach for better encoding spatial context into GNN-based predictive models. PE-GNN is highly modular and can work with any GNN backbone. PE-GNN contains a positional encoder (PE) [24,30], which learns a contextual embedding for point coordinates throughout training. The embedding returned by PE is then concatenated with other node features to provide the training data for the GNN operator. PE-GNN further predicts the Moran's I metric of spatial autocorrelation of the output as an auxiliary task in parallel to the main objective, expanding the approach proposed by [22] to continuous spatial coordinates. Lastly, we train PE-GNN by constructing a novel training graph, based on -nearest-neighborhood, from a randomly sampled batch of points at each training step. This forces PE to learn generalizable features, as the same point coordinate might have different spatial context (neighbors) at different training steps. Similarly, this training approach leads us to compute a "shuffled" Moran's I, implicitly nudging the model to learn a general representation of spatial autocorrelation which works across varying neighbor sets. Over a range of spatial regression tasks, we show that PE-GNN improves performance of different GNN backbones in most settings, except when very little training data is available.
The contributions of our study are as follows: • We propose PE-GNN, a novel GNN architecture including a positional encoder learning spatial context embeddings for each point coordinate to improve predictions. • We propose a novel way of training the positional encoder, distinct from [24]. We use the output of PE concatenated with other node features to predict an outcome variable. PE learns through backpropagation on the main regression loss in an end-to-end fashion. Training PE thus takes into account not only the eventual variable of interest, but also further contextual information at the current location-and its relation to other points. • We expand the Moran's I auxiliary task learning framework proposed by [22] for continuous spatial data. • Our training strategy involves the creation of a new training graph at each training step from the current, random point batch. This enables learning of a more generalizable PE embedding and allows computation of a "shuffled" Moran's I, which accounts for different neighbors at different training steps, thus overcoming the well-known scale sensitivity of Moran's I.
• To the best of our knowledge, PE-GNN is the first GNN based approach that is competitive with Gaussian Processes on pure spatial interpolation tasks, i.e., predicting a (continuous) output based solely on spatial coordinates, as well as substantially improving GNN performance on both spatial regression and spatial interpolation tasks.
2 RELATED WORK 2.1 Traditional and neural-network-based spatial regression modeling Our work considers the problem of modeling geospatial data. This poses a distinct challenge, as standard regression models (such as OLS) fail to address the spatial nature of the data, which can result in spatially correlated residuals. To address this, spatial lag models [3] add a spatial lag term to the regression equation that is proportional to the dependent variable values of nearby observations, assigned by a weight matrix. Likewise, kernel regression takes a weighted average of nearby points when predicting the dependent variable. Nevertheless, the most popular off-the-shelf methods for modeling continuous spatial data are based on Gaussian processes [6].
With the emergence of machine learning, there has been a rise of research on applications of neural network models for spatial modeling tasks. More specifically, graph neural networks (GNNs) are often used for these tasks with the spatial data represented graphically. Particularly, they offer flexibility and scalability advantages over traditional spatial modeling approaches. Specific GNN architectures including Graph Convolutional Networks [20], Graph Attention Networks [31] and GraphSAGE [12] are powerful methods for inference and representation learning with spatial data.
Recently, GNN approaches tailored to the specific complexities of geospatial data have been developed. The authors of Kriging Convolutional Networks [4] propose using GNNs to perform a modified kriging task. [12] apply GNNs for a spatio-temporal Kriging task, recovering data from unsampled nodes on an input graph.
We look to extend this line of research by providing stronger, explicit capacities for GNNs to learn spatial structures. Additionally, our proposed method is highly modular and can be combined with any GNN backbone.

Spatial context embeddings for geographic data
Through many decades of research on spatial patterns, a myriad of measures, metrics, and statistics have been developed to cover a broad range of spatial interactions. All of these measures seek to transform spatial locations, with optional associated features, into some meaningful embedding, for example, a theoretical distribution of the locations or a measure of spatial association. A set of point locations in 2 -space may be distributed randomly, or may follow some underlying spatial pattern. Spatial point processes are some of the oldest methods of analyzing spatial data and are a useful tool when point locations depend on one another, e.g., in the case of epidemiological or ecological data. Point processes are particularly popular for modeling geographic data, for example Cox processes [1] and Hawkes processes [29]. Empirical observations can be embedded using these point processes, for example, to make predictions on future spatial spread. If we now assume that point locations also come with a continuous feature value, the range of potential spatial effects and associated metrics is vastly expanded.
The most common metric for continuous geographic data is the Moran's I statistic, developed by [2]. Moran's I measures local and global spatial autocorrelation and acts as a detector of spatial clusters and outliers. The metric has also motivated several methodological expansions, like local spatial heteroskedasticity [26] and local spatial dispersion [33]. Measures of spatial autocorrelation have already been shown to be useful for improving neural network models through auxiliary task learning [22], model selection [21], embedding losses [23] and localized representation learning [10].
Beyond these traditional metrics, recent years have seen the emergence of neural network based embeddings for geographic information. [32] use kernel embeddings to learn social media user locations. [10] devise an approach using local point-of-interest (POI) information to learn region embeddings and integrate similarities between neighboring regions to learn mobile check-ins. [34] develop GPS2Vec, an embedding approach for latitude-longitude coordinates, based on a grid cell encoding and spatial context (e.g., tweets and images). [24] developed Space2Vec, another latitudelongitude embedding without requiring further context like tweets or POIs. Space2Vec transforms the input coordinates using sinusoidal functions and then reprojects them into a desired output space using linear layers. All of these studies highlight the potential of spatial coordinate embeddings (with or without spatial context enrichment) for improving predictive models for geographic data.

METHOD 3.1 Graph Neural Networks with Geographic Data
We elaborate on our method using Graph Convolutional Networks (GCNs) as an example backbone for our novel PE-GNN approach. Let us first define a datapoint = { , x i , c i }, where is a continuous target variable (scalar), x i is a vector of predictive features and c i is a vector of point coordinates, mapping the datapoint into geographic space (e.g., latitude and longitude values). We use the great-circle distance = ℎ (c i , c j ) between point coordinates to create a graph of all points = { 1 , . . . , } in the set. This can be done either by using a -nearest-neighbor approach or by drawing a boundary radius around each coordinate to define each point's neighborhood and create an adjacency matrix A. The graph = ( , ) consists of a set of vertices (or nodes) = { 1 , . . . , } and a set of edges = { 1 , . . . , } as assigned by A. Each vertex ∈ has respective node features x i and target variable . While the adjacency matrix A usually comes as a binary matrix (with values of 1 indicating adjacency and values of 0 otherwise), one can account for different distances between nodes and use point distances or kernel transformations thereof [4] to weight A. Now, let D define the degree matrix and I the identity matrix of our graph . We can then define the normalized adjacency matrixĀ as: As proposed by [20], a GCN layer can now be defined as where describes an activation function (e.g., ReLU) and W ( ) is a weight matrix parametrizing GCN layer . The input for the first GCN layer H (0) is given by the feature matrix X containing all node feature vectors x 1 , . . . , x n . The assembled GCN predicts the output Y = (X, Θ ) parametrized by Θ .

Context-aware spatial coordinate embeddings
In the traditional GCN setting, the only intuition for spatial context stems from connections between nodes which allow for graph convolutions, akin to pixel convolutions with image data. This restricts the capacity of the GCN to capture spatial patterns in several ways: • While defining good neighborhood structures can be crucial for GCN performance, this often comes down to somewhat arbitrary choices like selecting the nearest neighbors of each node. Without prior knowledge on the underlying data, the process of setting the right neighborhood parameters can be difficult and require extensive testing. Furthermore, a single value of might not be best for all nodes: different locations might be more or less dependent on their neighbors. • Assuming that no underlying graph connecting point locations is known, one would typically construct a graph using the distance (Euclidean or other) between pairs of points. In many real world settings (e.g., points-of-interest along a road network) this assumption is unrealistic and may lead to poorly defined neighborhoods. • GCNs contain no intrinsic tool to transform point coordinates into a different (latent) space that might be more informative for representing the spatial structure, with respect to the particular problem the GCN is trying to solve.
As such, GCNs struggle with tasks that explicitly require learning of complex spatial dependencies. For example, we show in our experiments that simple GCNs are not able to solve simple spatial interpolation tasks, i.e., predicting a continuous outcome variable from the point coordinates only. We propose a novel approach to overcome these difficulties, by devising a new positional encoder module, and learning a flexible spatial context encoding for each geographic coordinate, motivated by recent advances in transformers and spatial representation learning.
Given a batch of datapoints, we create the spatial coordinate matrix C from individual point coordinates c 1 , . . . , c n and define a positional encoder and a fully-connected neural network (Θ ), parametrized by Θ . Following the intuition of transformers [30] for geographic coordinates [24], the sinusoidal transform is a concatenation of scale-sensitive sinusoidal functions at different frequencies, so that with being the total number of grid scales and and setting the minimum and maximum grid scale (comparable to the lengthscale parameter of a kernel). The scale-specific encoder )] processes the spatial dimensions (e.g., latitude and longitude) of C separately, so that ; sin where = . The output from is then fed through the fully connected neural network (Θ ) to transform it into the desired vector space shape, creating the coordinate embedding matrix C = (C, , , Θ ).

Auxiliary learning of spatial autocorrelation
A further particularity of geographic data is that it often exhibits spatial autocorrelation: observations are related, in some shape or form, to their geographic neighbors. Measures of spatial autocorrelation have recently been successfully integrated into neural networks for discrete spatial data (images) [22,23], however no such approach exists for continuous spatial coordinates.
Spatial autocorrelation can be measured using the Moran's I metric of local spatial autocorrelation [2]. Moran's I captures localized homogeneity and outliers and thus functions as a detector of both spatial clustering and spatial change patterns. In the context of our problem, the Moran's I measure of spatial autocorrelation for our outcome variable is defined as: where , ∈ A denotes adjacency of observations and . As proposed by [22], predicting the Moran's I metric of the output can be used as auxiliary task during training. Auxiliary task learning [28] is a special case of multi-task learning, where one learning algorithm tackles two or more tasks at once. In auxiliary task learning, we are only interested in the predictions of one task; however, adding additional, auxiliary tasks to the learner might improve performance on the primary problem: the auxiliary task can add context to the learning problem that can help solve the main problem. This approach is commonly used, for example in reinforcement learning [9] or computer vision [14,15].
Translated to our GCN setting, we seek to predict the outcome Y and its local Moran's I metric (Y) using the same network, so that [Ŷ,( Y)] = (X). As [22] note, the local Moran's I metric is scale-sensitive and, due to its restriction to local neighborhoods, can miss out on longer-distance spatial effects [7,25]. But while [22] propose to compute the Moran's I at different resolutions, the GCN setting allows for a different, novel approach to overcome this issue.
Training of GCNs is often conducted using minibatches of points at each training step. Rather than constructing the graph of training points a priori, we opt for a procedure where in each training step, ℎ points are sampled from the training data as batch . A graph with corresponding adjacency matrix A is constructed for the batch, along with the Moran's I metric of the outcome variable (Y ). This approach brings a unique advantage: When training with (randomly shuffled) batches, points may have different neighbors in different training iterations. The Moran's I for point can thus change throughout iterations, reflecting a differing set of more distant or closer neighbors. This also naturally helps to overcome Moran's I scale sensitivity without the need to compute the metric at different resolutions. Altogether, we refer to this altered Moran's I as "shuffled Moran's I".

Positional Encoder Graph Neural Network (PE-GNN)
We now assemble the different modules of our method and introduce the Positional Encoder Graph Neural Network (PE-GNN). The modular approach allows for a seamless integration of the positional encoder and the spatial autocorrelation auxiliary task with any graph neural network architecture. The whole modeling pipeline of PE-GNN compared to a naive GNN approach is pictured in Figure 1. Taking GCN as the backbone architecture, PE-GCN is constructed as follows: Assuming a batch of randomly sampled points 1 , . . . , ℎ ∈ , a spatial graph is constructed from point coordinates c 1 , . . . , c ℎ using -nearest-neighborhood, resulting in adjacency matrix A . The point coordinates are then subsequently fed through the positional encoder (Θ ), consisting of the sinusoidal transform and a single fully-connected layer with sigmoid activation, embedding the 2 coordinates in a customizable latent space, returning vector embeddings c 1 , . . . , c ℎ = C . The neural network allows for explicit learning of spatial context, reflected in the vector embedding. We then concatenate the positional encoder output with the node features, to create the input for the first layer of our GCN as: The subsequent layers follow according to Equation 2.
To integrate the Moran's I auxiliary task, we compute the metric (Y ) for our outcome variable Y at the beginning of each training step according to Equation 5, using spatial weights from A . Prediction is then facilitated by creating two prediction heads, here linear layers, while the graph operation layers (e.g. GCN layers) are shared between tasks. Finally, we obtain predicted valuesŶ and(Y ). This auxiliary task learning process is also illustrated in Figure 1.
The loss of PE-GCN can be computed with any regression criterion, for example mean squared error (MSE): where denotes the auxiliary task weight. The final model is denoted as Θ ,Θ . Algorithm 1 describes a training cycle. Sample minibatch of ℎ points with features X , coordinates C and outcome Y .

5:
Construct a spatial graph with adjacency matrix A from coordinates C using -nearest neighbors 6: Using spatial adjacency A , compute Moran's I of output as (Y ) We begin training by initializing our model , for example a PE-GCN, with random weights on all layers and potential hyperparameters such as the PE embedding dimension. We also define our optimizer, setting hyper-parameters such as the learning rate. We now start the training cycle: At each training step, we first sample a minibatch of points from our training data. These points come as features X , point coordinates C and outcome variables Y . We then construct a graph from spatial coordinates C usingnearest-neighborhood, obtaining an adjacency matrix A . Next we use A as spatial weight matrix to compute local Moran's I values (Y ) from Y . As minibatches are randomly sampled, this creates a "shuffled" version of the metric. We then run inputs X , C , A through the two-headed model Θ ,Θ obtaining predictionŝ Y , (Ŷ ). We then compute the loss L (Y , (Y ),Ŷ , (Ŷ ), ), weighing the Moran's I auxiliary task according to weight parameter . Lastly, we use the loss L to update our model parameters Θ , Θ according to stochastic gradient descent. Training is conducted for (or until convergence) and after it is finished, the final model is returned.
PE-GNN, with any GNN backbone, helps to tackle the particular challenges of geographic patterns, outlined above: • While our approach still includes the somewhat arbitrary choice of -nearest neighbors to define the spatial graph, the proposed positional encoder network is not bound by this restriction, as it does not operate on the graph. This enables a separate learning of context-aware embeddings for each coordinate, accounting for neighbors at any potential distance within the batch. • While the spatial graph used still relies on pre-defined distance measure, the positional encoder embeds latitude and longitude values in a high-dimensional latent space. These high-dimensional coordinates are able to reflect spatial complexities much more flexibly and, added as node features, can communicate these throughout the learning process. • Batched PE-GNN training is not conducted on a single graph, but a new graph consisting of randomly sampled training points at each iteration. As such, at different iterations, focus is put on the relationships between different clusters of points. This helps our method to generalize better, rather than just memorizing neighborhood structures. • The differing training batches also help us to compute a "shuffled" version of the Moran's I metric, capturing autocorrelation at the same location for different (closer or more distant), random neighborhoods.

EXPERIMENTS 4.1 Data
We evaluate PE-GNN and baseline competitors on four real-world geographic datasets of different spatial resolutions (regional, continental and global): California Housing 1 : This dataset contains the prices of over 20, 000 California houses from the 1990 U.S. census and was introduced by [18]. The regression task at hand is to predict house prices using features x (e.g., house age, number of bedrooms) and location c. California housing is a standard dataset for assessment of spatial autocorrelation.
Election 2 : This dataset contains the election results of over 3, 000 counties in the United States and was proposed by [16]. The regression task here is to predict election outcomes using sociodemographic and economic features (e.g., median income, education) x and county locations c.
Air temperature 3 : The air temperature dataset [13] contains the coordinates of 3, 000 weather stations around the globe. For this regression task we seek to predict mean temperatures from a single node feature , mean precipitation, and location .
3d Road 4 : The 3d road dataset [17] provides 3-dimensional spatial co-ordinates (latitude, longitude, and altitude) of the road network in Jutland, Denmark. The dataset comprises over 430, 000 points and can be used for a regression task where altitude is predicted using latitude and longitude coordinates c. Note that this dataset does not contain any node features x, so it is only used for the spatial interpolation task.

Baselines
We compare PE-GNN with four different graph neural network backbones: The original GCN formulation, introduced by [20] and outlined in the Methods section, graph attention mechanisms (GAT) [31] and GraphSAGE [12]. Lastly, we use Kriging Convolutional Networks (KCN) [4]. KCN differs from GCN primarily in two ways: it transforms the distance-weighted adjacency matrix A using a Gaussian kernel and adds the outcome variable and features of neighboring points to the features of each node. To avoid a data leak, test set points can only access neighbors from the training set to extract these features. We compare the naive version of all these approaches to the same four backbone architectures augmented with our PE-GNN modules. Beyond GNN-based approaches, we also compare PE-GNN to the most popular method for modeling continuous spatial data: Gaussian processes. For all approaches, we compare a range of different training settings and hyperparameters, as discussed below.

Experimental setup
To allow for a fair comparison between the different GNN approaches, we equip all models with the same architecture, consisting of two GCN / GAT / GraphSAGE layers with ReLU activation and dropout, followed by a final linear layer as regression head. The KCN model also uses GCN layers, following the author specifications. We test four different auxiliary task weights = {0, 0.25.0.5, 0.75}, where = 0 implies no auxiliary task. Spatial graphs are constructed assuming = 5 nearest neighbors, following results from previous work [4,16] and our own rigorous testing. We include a sensitivity analysis of the parameter and different batch sizes in our results section. Training for the GNN models is conducted using PyTorch [27] and PyTorch Geometric [8]. We use the Adam algorithm to optimize our models [19] and the mean squared error (MSE) as loss function. Gaussian process models (exact and   Table 1: Spatial Interpolation: Test MSE and MAE scores from four different datasets, using four different GNN backbones with and without our proposed architecture. approximate) are trained using GPyTorch [11]. Due to the size of the dataset, we only provide an approximate GP result for 3d Road. All training is conducted on single CPUs. The code for PE-GNN can be accessed here: https://github.com/konstantinklemmer/pe-gnn.

Model
Cali  Table 2: Spatial Regression: Test MSE and MAE scores from three different datasets, using four different GNN backbones with and without our proposed architecture.

Results
We test our methods on two tasks: Spatial Interpolation, predicting outcomes from spatial coordinates alone, and Spatial Regression, where additional node features are available. The results of our experiments are shown in Table 1 and 2. We evaluate both models using the mean squared error (MSE) and mean absolute error (MAE) metrics.
For the spatial interpolation task, we observe that the PE-GNN approaches consistently and vastly improve performance for all four backbone architectures across the California Housing, Air Temperature and 3d Road datasets and by a small margin for the Election dataset.
For the spatial regression task, we observe that the PE-GNN approaches consistently and substantially improve performance for all four backbone architectures on the California Housing and Air Temperature datasets. Performance remains unchanged or decreases by very small margins in the Election dataset, except for the KCN backbone which benefits tremendously from the PE-GNN approach, particularly with auxiliary tasks.
Our results are interesting in several ways: Generally, PE-GNN substantially improves over baselines in regression and interpolation settings. Most of the improvement can be attributed to the positional encoder, however the auxiliary task learning also has substantial beneficial effects in some settings, especially for the KCN models. The best setting for the task weight hyperparameter seems to heavily depend on the data, which confirms findings by [22]. To our knowledge, PE-GNN is the first GNN-based learning approach that can compete with Gaussian Processes on simple spatial interpolation baselines, though especially exact GPs still sometimes have the edge. Nonetheless, PE-GNN is of course substantially more scalable than exact GPs, which rely on expensive pair-wise distance calculations across the full training dataset. Due to this problem, we do not run an exact GP baseline for the highdimensional 3d Road dataset. For KCN models, which include neighbor information as additional node features, we observe a proneness to overfitting. As the authors of KCN mention, this effect diminishes in large enough data domains [4]. For example, KCNs are the best performing method on the 3d Road dataset-by far our largest experimental dataset. We observe that in cases when KCN learns well, PE-KCN can still improve its performance, as highlighted by the 3d road experiments. The KCN experiments also highlight the strongest effects of the Moran's I auxiliary tasks: In cases when KCN overfits (Election, Cali. Housing datasets), PE-KCN without auxiliary task ( = 0) is not sufficient to overcome the problem. However, adding the auxiliary task ( = [0.25, 0.5, 0.75]) can mitigate most of the overfitting issue. This directly confirms a theory of [22] on the beneficial effects of auxiliary learning of spatial autocorrelation. Throughout our experiments, PE-GNN provides a powerful and flexible framework for learning with geogpraphic coordinates. The positional encoder allows for the learning of spatial context features to inform the downstream task, while the Moran's I auxiliary task helps to overcome potential overfitting issues.
We also assess the robustness of PE-GNN training cycles. Figure  3 highlights the confidence intervals of PE-GNN models with GCN, GAT and GraphSAGE backbones trained on the California Housing and 3d Road datasets, obtained from 10 different training cycles. We can see that training runs exhibit only little variability. These findings thus confirm that PE-GNN can consistently outperform naive GNN baselines. Lastly, Figure 4 highlights some results from our sensitivity analyses with the and ℎ (batch size) parameters. After rigorous testing, we opt for = 5-NN approach to create the spatial graph and compute the shuffled Moran's I across all models. We chose ℎ = 2048 for Cali. Housing and 3d Road datasets and ℎ = 1024 for the election dataset.

CONCLUSION
With PE-GNN, we introduce a flexible, modular new GNN-based learning framework for geographic data. PE-GNN leverages recent findings in embedding spatial context into neural networks to improve predictive models. Our empirical findings confirm a strong performance. This study highlights how domain expertise from applied academic disciplines like geography can help improve machine learning models for the particular data that is common in these disciplines. We hope to build on the foundations of PE-GNN to develop further methods for geospatial machine learning.