{'fraction_non_alphanumeric': 0.031185243930517536, 'fraction_numerical': 0.009184169459438717, 'mean_word_length': 5.065166795965865, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Auxiliary-task learning for geographic data with autoregressive embeddings

Machine learning is gaining popularity in a broad range of areas working with geographic data, such as ecology or atmospheric sciences. Here, data often exhibit spatial effects, which can be difficult to learn for neural networks. In this study, we propose SXL, a method for embedding information on the autoregressive nature of spatial data directly into the learning process using auxiliary tasks. We utilize the local Moran's I, a popular measure of local spatial autocorrelation, to"nudge"the model to learn the direction and magnitude of local spatial effects, complementing the learning of the primary task. We further introduce a novel expansion of Moran's I to multiple resolutions, thus capturing spatial interactions over longer and shorter distances simultaneously. The novel multi-resolution Moran's I can be constructed easily and as a multi-dimensional tensor offers seamless integration into existing machine learning frameworks. Throughout a range of experiments using real-world data, we highlight how our method consistently improves the training of neural networks in unsupervised and supervised learning tasks. In generative spatial modeling experiments, we propose a novel loss for auxiliary task GANs utilizing task uncertainty weights. Our proposed method outperforms domain-specific spatial interpolation benchmarks, highlighting its potential for downstream applications. This study bridges expertise from geographic information science and machine learning, showing how this integration of disciplines can help to address domain-specific challenges. The code for our experiments is available on Github: https://github.com/konstantinklemmer/sxl.


INTRODUCTION
When monitoring the physical environment, from crop yields to air pollution, the gathered data are often geographic in nature and follow some spatial process: data values depend on their spatial locations. This violates a key assumption of many statistical learning frameworks, that data are independent and identically distributed (iid). Even models which are explicitly designed to account for these dependencies, such as kernel methods, may struggle with over-and under-smoothing of spatial patterns or with distinguishing stationary from non-stationary spatial effects. The complexities of geographic data concern researchers in many academic fields and are a key focus of the geographic information sciences (GIS). The GIS community has a long tradition of analyzing spatial phenomena, developing metrics to measure spatial effects and deploying models to account for spatial dependencies [1,12,29]. With the growing popularity of scalable machine learning methods, particularly deep neural networks, applications of these approaches to geospatial data domains have become more and more common. Nevertheless, these applications have only rarely inspired methodological innovation [22] in neural networks. Recently, [39] specifically pointed to a lack of deep learning methods tailored to geospatial and spatiotemporal data in the context of earth system science. In light of this call-to-action, we propose SXL, a novel neural network method for geospatial data domains, that is inspired by domain expertise from GIS and explicitly learns spatial dependencies contained within the data. This is facilitated through a multi-task learning process, where a spatial embedding capturing local autoregressive effects at each data point is learned as an auxiliary task, sharing the model parameters with the main task and hence nudging the model to learn both the original data and its embedding in parallel. We integrate one of the most prominent and widely used metrics in GIS into the model: the Moran's I measure of local spatial autocorrelation [2]. To also account for longer-distance spatial relations, we propose a novel multi-resolution local Moran's I by gradually coarsening the input data. By providing a learner with prior knowledge on the autoregressive nature of the data we seek to improve performance on a broad range of spatial modeling tasks.
Our main contributions can be summarized as follows: (1) We propose a novel, flexible multi-resolution expansion of the Moran's I measure of local spatial autocorrelation. 2) We use the traditional arXiv:2006.10461v3 [cs.LG] 19 Aug 2021 and multi-resolution local Moran's I as embeddings to be learned in an auxiliary learning procedure to capture both short-and longdistance spatial effects. (3) We provide a practical framework to adapt our method to arbitrary neural network architectures and different supervised (predictive) and unsupervised (generative) spatial learning tasks. (4) For the purpose of balancing the losses of multiple tasks in a generative modeling setting, we develop a novel auxiliary task GAN loss based on uncertainty weights [10]. (5) We evaluate SXL on both generative and predictive spatial modeling tasks, providing empirical evidence for consistent and robust performance gains across multiple synthetic and real-world experimental settings.

RELATED WORK 2.1 Machine learning with geospatial data
Geography has a long tradition of modeling and empirical analysis. Ideas from GIS and spatial statistics have inspired popular approaches in modern machine learning, from Gaussian Processes (GPs), pioneered by the development of Kriging in the 1960s [31], to spatial scan statistics [26], widely used for event and pattern detection. Nevertheless, since the emergence of the era of deep neural networks, the relationship between the GIS and machine learning communities have been mostly defined through applications of existing neural network models to geographic data. Only rarely have concepts and ideas from GIS and spatial statistics motivated methodological advancements in neural networks-whereas application areas such as vision, bioinformatics, and computational linguistics have strongly influenced the deep learning state-of-theart we know today. Nevertheless, some recent advances in machine learning have been exceptionally useful to the GIS community. Kernel methods such as GPs have seen huge progress towards overcoming computational bottlenecks stemming from the complexities of working with pair-wise distance matrices. Kernel interpolation techniques and GPU acceleration have allowed GPs to easily scale up to a million data points [15,46]. In neural networks, the emergence of graph-based networks and particularly graph convolutions have allowed for the modeling of asymmetric and non-Euclidean spatial relationships [28], and the emergence of physics-informed deep learning has reinforced the need for neural networks to model complex spatio-temporal patterns [47].
Nevertheless, as geographic context has proven to be of relevance in many machine learning applications [3,5,9,51], core concepts from GIS and geography have gradually attracted more attention in the ML community. While methodologies for geospatial data and problems have previously utilized traditional machine learning approaches like tree-based models [24] or ensemble learners [23], focus has recently switched more towards neural networks: for example, [30] and [53] propose a vector embedding to capture spatial context, inspired by word embeddings such as Word2Vec [33]. [54] propose to model complex spatial covariance patterns through injective warping functions, learned by a deep net. While there exist approaches for capturing spatial autocorrelation in neural networks [56], the Moran's I metric has previously only been used once in an explicit machine learning setting: [25] propose to use it as a heuristic for early stopping in generative adversarial network (GAN) training, but not explicitly for learning. More generally, geographic metrics describing spatial phenomena have barely been used in neural network models-which we believe is a missed opportunity.

Learning with auxiliary tasks
In this paper, we use auxiliary learning, an approach using multitask learning to improve performance on a primary task, originally conceptualized by [43]. The authors propose to give learners "hints" related to the original task throughout training in order to improve training speed and model performance. This can be understood as forcing the learner (e.g., a neural network) to focus its attention on certain patterns in the data, highlighted by the auxiliary objective. It also implies that the auxiliary task has to provide some meaningful embedding of the primary task. Auxiliary learning has been particularly successful in deep reinforcement learning and is widely utilized [14]: for example, learning to steer a wheel can be improved using auxiliary tasks related to image segmentation and optical flow estimation [20]. These auxiliary tasks are already related to spatial perception, which is common when working with image data: recent work has highlighted the applicability of pixel control tasks [21] or depth estimation [34]. Auxiliary tasks have also been successfully deployed for generative modeling with GANs: Auxiliary Classifier GAN (AC-GAN) [36] extends the original GAN loss function by an auxiliary classifier to improve the fidelity of generated images. However, this approach comes with complications. As recent work [16] notes, AC-GANs can lead to perfect separability-where the GAN discriminator is easily winning the two-player game against the generator, preventing efficient and balanced learning. Another use of spatial semantics in an auxiliary task GAN setting is proposed by [48]: the authors propose to generate surface normal maps, a 3D representation of the 2D image, and use these in an auxiliary task. In an explicit geospatial setting, auxiliary tasks have been used in semantic visual localization [40] and semantic segmentation [7], in addition to dedicated spatial and spatio-temporal multi-task frameworks [50]. Nevertheless, to the best of our knowledge, measures of spatial autocorrelation such as Moran's I have never previously been used in any kind of multi-task learning setting, be it for generative or predictive spatial modeling.

METHODOLOGY 3.1 Multi-resolution local Moran's I
As outlined above, working with geospatial data requires a careful assessment of and accounting for potential autoregressive effectsan intuition which neural networks traditionally do not provide. One of the most prominent measures, capturing spatial autocorrelation at the point-level, is the local Moran's I metric [2]. Local Moran's I measures the direction and extent of similarity between each observation and its local spatial neighbourhood. As such, it provides an indication for both local spatial clusters and spatial outliers.
Formally, let be a 2-spatial matrix (image)  define the spatial neighbourhood of observation to be N = { ∈ : , > 0}. Here, , corresponds to a binary spatial weight matrix, indicating whether any observation is a neighbor of . For continuous data, the creation of this matrix requires computing a pair-wise distance matrix or kd-tree, however in our case using discrete × matrices, this problem is trivial. Throughout this study, we utilize queen contiguity (i.e. all adjacent grid cells, including diagonals, are neighbors), but the approach generalizes to arbitrary neighborhood definitions. We can compute the local Moran's I statistic, , of observation with the mean over all observationsā s: Combining all local Moran's I values gives the matrix × , of the same size as × . can take positive or negative values: a positive value suggests that a data point is similar to its neighbors, which could indicate latent cluster structure. A negative value suggests that the data point is distinctly different from neighboring data points, which could indicate a changepoint or edge. While the Local Moran's I statistic is closely correlated to its input, their relationship can take different forms, depending on the complexity of the inputs' spatial structure. Figure 1 provides examples of different geospatial data and their respective ( ) representations. One of the main limitations of the local Moran's I metric is its restriction to represent local spatial dependence only at the scale provided by the immediate neighbors as defined by the neighborhood matrix . Thus longer-range spatial dependencies can be lost. This scale sensitivity of the local Moran's I is known as a common challenge in applications [13,32,55]. Nevertheless, to our knowledge, no alternative metric accounting for this issue exists. Here, we propose a novel, multi-resolution representation of the local Moran's I by increasingly coarsening the input data for the Moran's I computation and then upsampling the output back to the original data shape. The coarsening step here is analogous to a 2-average pooling operation, so that our coarsened input is given as: where gives the kernel size and the subscript indicates downsampling. The coarsened spatial matrix corresponds to the vector x ( ) of length ( ) . The coarsened local Moran's I, (x ( ) ), is then computed according to Equation (2), using the spatial weight matrix ( ) , corresponding to the new size / × / of the downsampled input. In the last step, the coarsened Moran's I is upsampled again to the original input size × using nearest-neighbor interpolation. This whole process can be repeated several times to compute the local Moran's I at increasingly coarse resolutions. The local Moran's I values at different resolutions can then be stacked on top of one another, much like a multi-channel image (e.g. RGB image). As such, tensors provide an ideal data structure for our metric. We illustrate this with an example in Figure 2.

Auxiliary learning of spatial autoregressive structures
Auxiliary task learning shares the benefits of multi-task learning: auxiliary tasks hint at specific patterns in the data for the model to focus attention on. Further, they introduce a representation bias, whereas the model prefers latent representations of the data that work for both primary and auxiliary tasks, thus helping with generalization. Lastly, auxiliary tasks can work as regularizers by introducing inductive bias and decreasing the risk of overfitting the model. Here, we want to use the local Moran's I embedding and our newly introduced multi-resolution Moran's I as auxiliary tasks. The main motivation for any auxiliary tasks is "relatedness" to the primary task: spatial theory characterizes a spatial pattern as a reflection of underlying spatial processes. Accordingly, [8] concludes that "[...] the capability of generalizing and quantifying spatial patterns is a prerequisite to understanding the complicated processes governing the distribution of spatial phenomena. "-explicitly mentioning the power of the Moran's I metric to capture these effects. This statement can be translated directly into a learning algorithm, where the learning of a spatial pattern is constrained by a simultaneous learning of the underlying spatial process. Together with the well documented success of spatial auxiliary tasks in computer vision (see Related Work), this makes auxiliary tasks based on the local Moran's I well motivated by both spatial theory and machine learning research. Recent research further highlights the importance of learning at multiple resolutions to support a comprehensive understanding of spatial processes [37,44,49]. Lastly, the local Moran's I (and the multi-resolution variant) can be constructed for any numerical input, and can thus be seamlessly integrated in a broad range of learning settings and with arbitrary neural network architectures.
With our experiments, we focus on two distinct settings: generative spatial modeling using GANs [17], and predictive spatial  modeling in the form of spatial interpolation. To outline the application of our proposed auxiliary task approach, we introduce the GAN example in detail-the predictive modeling formulation follows from this straightforward. GANs are a family of generative models comprised of two neural networks, a Generator G that produces fake data and a Discriminator D that seeks to distinguish between real and fake data. These two networks are agents in a two-player game, where G learns to produce synthetic data samples that are faithful to the true data generating process, and D learns to separate real from fake samples, thus pushing G to produce increasingly realistic synthetic data. The standard GAN loss function is thus given as: consisting of the Discriminator and Generator losses Our auxiliary task approach augments the Discriminator with a loss based on the Moran's I embeddings of the real and the fake data: so that the composite loss for auxiliary tasks (single-or multiresolution) is given as: Both loss functions use a customary weight hyper-parameter for the auxiliary losses. Alternatively, we could fit separate weights for each auxiliary task. The approaches are further illustrated in Figure 3. Integrating the auxiliary tasks into predictive models for spatial interpolation is more straightforward: We simply let a regressor predict the (multi-resolution) local Moran's I of the output, ( ) ∼ ( ( )), simultaneously with the main task ∼ ( ).

Loss balancing using task uncertainty
In addition to the hard loss weight outlined above, we also test a setting for automatically learning the loss weight of main and auxiliary tasks. For this, we follow the approach conceptualized by [10] and utilize each task's homoskedastic uncertainty to inform the loss weight. To formalize this, [10] define a probabilistic multi-task regression problem with tasks and likelihood where y 1 , . . . , y are the main and auxiliary model outputs and x is the model input. Using the maximum likelihood method, the minimization objective of the multi-task regression is min L ( 1 , . . . , ): where 1 , . . . , are the model noise parameters. Minimizing this objective can be interpreted as learning the relative weight of each task's contribution to the composite loss. The noise is kept from decreasing infinitely by the last term of the loss, which serves as a regularizer. This loss constitutes the objective we use for our predictive modeling task. However, no approach for uncertainty weighted auxiliary task GANs exists. We hence propose an extension to the auxiliary GAN loss outlined in the previous section: Instead of using a fixed auxiliary task weight , we augment both main and auxiliary discriminator losses with uncertainty weights, so that This constitutes the first adaption of the uncertainty task balancing principles to the multi-task GAN family. In the following experiments, we report results from both fixed weight parameters and uncertainty weights. Throughout rigorous experiments, we find weights of = [0.1, 0.01] to be particularly helpful and hence utilize these throughout our experiments. We use auxiliary learning settings with hard-parameter sharing, where the top layers of the respective models are task-specific. For the GAN experiment, just the last layer is task-specific, and for the interpolation experiment the last two layers are task specific.
Having outlined our approach for GANs in detail, the adaptation of SXL to predictive spatial models is trivial: Just as with the GAN discriminator, a second prediction head is added to the predictive model (e.g. a neural network regressor for spatial interpolation), aiming to predict the (multi-resolution) local Moran's I embedding of the output parallely to the main task. The losses of both tasks are then combined the same way as for the GAN discriminator, using either a hard weighting parameter or uncertainty weights, as outlined in Equation 10.

EXPERIMENTS 4.1 Experiment 1: Generative spatial modeling
Selection & Evaluation: In our first experiment, we want to examine whether our proposed method can improve the learning of spatial data generating processes. We generate synthetic data from several types of GANs, with and without including the Moran's I auxiliary tasks, and compare how faithful the generated samples are compared to the true distribution of samples. To assess model quality, we use the Maximum Mean Discrepancy (MMD) metric [6], a distance measure between distributions based on mean embeddings of the features. For data distributions and , the MMD is defined as ( , ) = || − || R . The empirical MMD for random variables and of length is given as where : X × X represents a positive-definite kernel-in our case a radial basis function (RBF) kernel. The more similar the data distributions and are, the closer the MMD metric gets to 0. A lower MMD score between samples of real and synthetic data thus indicates higher quality of the synthetic samples. As GAN training is notoriously difficult and prone to mode collapse and other issues, we opt for the following training and selection procedure to ensure our findings are robust: For each architecture and training strategy combination (e.g. EDGAN MRES MAT UW) we train ten cycles of GANs. For each cycle, we save the generator which optimizes the MMD metric on a held-out validation set (separate from the training data and from the test data used for evaluation), rather than choosing the final model after all training iterations. We then choose the one of these ten generators that optimizes the MMD metric on the validation set, and finally evaluate the MMD score of that generator on a separate, held-out test set. In short, this training process allows each architecture and training strategy combination ten cycles to train the best possible generator (as measured on the validation set), which is then evaluated on the test set and compared to all other combinations of architecture and training strategy.
Data: We select four datasets for our experiments: (1) A toy dataset of 7000 32 × 32 tiles with a Gaussian peak and a Gaussian dip, where the position of the dip mirrors the position of the peak.
(2) The PetrelGrid seabed relief dataset [27], processed into a grid of 195 32 × 32 tiles. (3) Digital Elevation Model (DEM) data of the area surrounding Lake Sunapee (NH, USA), as found in the elevatr R package 1 , processed into a grid of 1156 32 × 32 tiles. (4) Tree canopy data of the University of Maryland's "Global Tree Change" project [18], processed into a grid of 1800 64 × 64 tiles. These datasets are chosen to represent a range of different geospatial patterns occurring in real-world physical environments and relate to important modeling challenges in the earth sciences, ecology, or geography. For more information on the data used, including summary statistics, please refer to Appendix A.
Benchmark Models: The modularity of our proposed auxiliary task learning method allows us to test it on a range of different GAN models. We chose the original GAN implementation, denoted here as GAN [17], the DCGAN [38] and lastly an Encoder-Decoder GAN (EDGAN ) architecture recently proposed by [57] and explicitly designed for geospatial applications. All models are optimized using the same, traditional GAN objective lined out in the previous section. We test all benchmark models with the single-and multi-resolution Moran's I auxiliary task (MAT / MRES-MAT) as well as with fixed ( ) and uncertainty (UW) based task weights. For more details on    Table 2 from the hard loss weights . We can see that for both strategies the auxiliary task settings improve performance for most experiments, agnostic of the underlying GAN architecture, by usually 3-10%. We believe the auxiliary tasks support the learning process in two ways: (1) GAN Discriminators are known to exhibit spare capacity (i.e., they can be too powerful), which can cause them to over-specialize, leading to worse generalization performance [19]-thus adding a second, closely correlated task should not pose a problem but might help prevent over-specialization in the bottom Discriminator layers. (2) The losses stemming from the auxiliary tasks have a regularizing effect throughout training, further preventing Discriminator over-specialization.
This leads to several beneficial effects, some even visually apparent, as highlighted in Figure 4. Figure 4a shows how generators trained with MRES MAT appear to be better at smoothing (but not over=smoothing) spatial artifacts-residual, noisy spatial patterns and hard edges introduced by the model-compared to the same GAN backbone trained without our auxiliary task. Figure 4c shows two EDGANs, one trained with and one trained without MRES-MAT. We observe that the model trained without the auxiliary task exhibits "mode collapse", a phenomenon common with GANs where a Generator always produces the same image or some parts of the image are always the same, while the Generator including the auxiliary task does not exhibit this behavior. It is important to note that these examples are not cherry-picked but represent a pattern  Table 2: Test MMD scores (lower is better) of different GAN architectures. We compare synthetic samples from these generators to held-out test data to compute the scores. Shown are models trained with hard auxiliary task weights .
we can observe throughout all our experiments: the auxiliary tasks consistently improved model performance.
However, the optimal setting for applying the auxiliary task appears to vary, both in terms of the task weighting strategy and using the simple MAT versus the MRES-MAT. For example, while for the Vanilla GAN architecture, single-resolution MAT with hard task weights seems superior, both DCGAN and EDGAN appear to benefit especially from MRES-MAT with uncertainty weights. We also observe cases where adding the auxiliary tasks with a particular weighting strategy massively increases the MMD score (e.g. Toy dataset, GAN-MAT UW). This can happen when the auxiliary task "overpowers" the main task, causing the Generator to produce synthetic Moran's I embeddings. This further justifies the generator selection strategy employed throughout our experiments. Task weighting: Tables 1 and 2 show that models employing our MAT / MRES MAT auxiliary tasks can reliably produce generators that outperform naively trained models. Uncertainty weight models with MRES MAT represent the "winning" Generator in 9 of 12 cases. For the hard loss weights, we are always able to find a combination of and MAT / MRES MAT that outperforms naive training. We thus conclude that while training strategies appear to be highly data and model dependent, one can find a performanceincreasing setting in almost all cases. Here, the MRES MAT UW strategy seems to be the safest bet, as it does not require further, manual weight parameter tweaking. Figure 5 shows the losses and learned task uncertainties in an example training cycle of an EDGAN MRES MAT UW model using the DEM dataset. More details on all our data, model architectures and training settings can be found in the Appendix.

Experiment 2: Predictive spatial modeling
Evaluation: In the second experiment, we tackle spatial interpolation, that is, obtaining high-resolution spatial data from a lowresolution input. Spatial interpolation is widely used in real world applications, for example with meteorological measurements [45], air quality assessment [52] or mobile sensing [35]. It is a regression task and can be evaluated using the residual mean squared error (RMSE) between real and predicted high-resolution output. As such, this task is comparable to image super-resolution, a popular task in computer vision. Nevertheless, spatial interpolation is particularly focused on reconstructing the spatial patterns of the output. We again train 10 models per strategy and compare their performance when no model selection is used (final model used for prediction on test set) and when model selection on a validation set is applied, saving the 10 best models, one from each run.
Data: As a common use case in geography and ecology, we use hillshades of DEM data from the National Ecological Observatory Network (NEOS) 2 for the interpolation task. The data is pre-processed into a grid of 1674 64 × 64 tiles (output) and a low-res 32×32 version (input) by removing every second row and column of the image matrix. For more information on the data used, including summary statistics, please refer to Appendix A. Benchmark Models: Spatial statistics provides a range of tools to tackle interpolation problems. Commonly used methods we focus on here are: (1) Bicubic interpolation (BicInt), commonly used for interpolating 2-regular grids, (2) Inverse Distance Weighting (IDW ) [41], a weighted rolling-average approach, (3) Ordinary Kriging (OK) [11], spatial interpolation closely related to Gaussian Process regression and (4) Universal Kriging (UK) [42], a generalization of OK assuming a polynomial trend model. We compare these established methods to a simple CNN implementation with two hidden layers (5). Again, the modularity of SXL allows us to simply plug-in our MAT by having the CNN interpolate the spatial pattern and its Moran's I embedding from low to high resolution, making the last layers task-specific. In this setting, we do not use the MRES MAT, as our experiments show that further coarsening the already-reduced image is counterproductive. The CNN model main tasks are optimised using MSE loss, while the auxiliary tasks use ℓ1 loss. For more details on the experimental setup, please refer to Appendix B.II.
Findings: The results of our experiments are presented in Table  3 and Figure 6. We can again see a positive effect of the MAT on the performance of the CNN model-outperforming all other benchmarks. If no model selection is deployed, both hard loss weights and uncertainty weights produce models that outperform the naively trained CNN. MAT UW models provide the best average performance increase, of around 5%. If model selection is utilized, the MAT UW strategy outperforms the naive CNN by about 1.5%. Both of these performance increases are statistically significant, according to a paired t-test of the mean prediction scores.
Task weighting: Finding the optimal task weighting strategy appears much less tricky as for generative modeling, as the MAT UW strategy prevails in all interpolation experiments, whether model selection is applied or not. We can thus conclude that the MAT UW both consistently and significantly improves training. More details on all our data, model architectures and training settings can be found in the Appendix.

CONCLUSION
In general, our experiments give some insight into the way the auxiliary learning mechanism works, allowing us to compare the method to related ideas in machine learning: First, as mentioned before, we believe the auxiliary tasks to have a regularizing effect on the learning process, preventing models from overfitting on the primary task by forcing them to follow "spatial rules". Second, we believe the MRES MAT shares the intuition of moment matching, as we seek to simultaneously minimize the loss of one function at several coarsened resolutions. Third, the MRES MAT also shares the same goal as recent developments in visual self-attention: moving beyond the short distance spatial learning of convolutional layers and accounting for longer distance spatial effects. We also make some empirical observations with respect to the most important design choices of our models: (1) The use of MAT vs. MRES MAT appears to depend on both the input data and the model architecture used, with MRES MAT prevailing in most cases. (2) The optimal auxiliary task weighting strategy varies across generative and predictive modeling experiments, but in most cases uncertainty weights appear to have the edge. In summary, with SXL we propose the use of single-and multiresolution measures of local spatial autocorrelation for improving the learning of geospatial processes. We introduce a novel, flexible multi-resolution version of the local Moran's I statistic using coarsened inputs. We demonstrate its integration as an auxiliary task into generative and predictive neural network models, using both hard (static) and task uncertainty (automatically learned) loss weights. We empirically show robust, consistent and significant performance gains of up to 10% for generative spatial modeling and up to 5% for predictive spatial modeling when using this strategy. We comment on the importance of the exact configuration of the auxiliary tasks, especially choosing single-versus multi-resolution auxiliary tasks and the weighting strategy for auxiliary losses. We see this study further evidence of the importance of integrating domain expertise from GIS into neural network methods for geospatial data, opening a broad range of further research directions: In future work, we plan to expand this idea beyond discrete, regularly gridded data to continuous spatial processes and to assess the applicability of other common spatial metrics for a more comprehensive learning of complex geospatial patterns.

APPENDIX A: Data Description
Due to a lack of geospatial benchmark datasets within the machine learning community, we run our experiments using one toy dataset and three datasets from real-world geospatial applications. All data is chosen to represent different spatial patterns and to be closely related to important applications in fields such as climate science or geology.
Toy: The Toy dataset consists of 32 × 32 matrices with values coming from a function creating a Gaussian peak at a random location, which is mirrored diagonally by a Gaussian dip. This function is given as: where c X and c Y are the spatial coordinates mapping the values to the 32 × 32 matrix (so in our cases, integers in the range [0, 31], is a positive constant determining the size of the Gaussian peak and dip (we use = 7), and are random draws from integers in the range [0, 10], determining the location of the Gaussian peak and = 10 − and = 10 − are the location of the Gaussian dip, mirroring the peak diagonally.
PetrelGrid: The PetrelGrid dataset [27] is composed of geo-referenced seabed relief data. It can be accessed via R here: https://rdrr.io/cran/ spm/man/petrel.grid.html DEM (Generative Modeling): We use two different digital elevation model (DEM) based datasets, one for the generative modeling experiments and one for the spatial interpolation experiments. The DEM for generative modeling is chosen as it is rather small, enabling us to assess how our proposed method deals with data scarcity. An applicably small DEM dataset, providing a DEM of the area surrounding Lake Sunapee, NH, USA can be found as part of the elevatr R package; accessible via: https://rdrr.io/cran/elevatr/ man/lake.html.
TreeCanopy: This dataset contains data on global forest coverage. We use tree canopy, which describes canopy closure for all vegetation taller than 5m in height. The data comes from the University of Maryland's "Global Forest Change" project [18], documenting the global loss of forests in the light of climate change and forest exploitation. Specifically, we use data within the geographic area 50-60N / 100-110W; an area lying in continental Canada and representing a broad range of forest coverage types. The data can be accessed via: http://earthenginepartners.appspot.com/science-2013-global-forest/download_v1.6.html.

DEM (Spatial Interpolation):
The second, larger DEM dataset used for the spatial interpolation experiments is part of a LiDAR data collection conducted by the National Ecological Observatory Network (NEOS). Specifically, we use DEM hillshades from the NEOS training exercise outlined here: https://www.neonscience. org/da-viz-neon-lidar-co13flood-R. Hillshades are used to visualize terrain as shaded reliefs, where shades depend on a (synthetic) light source (e.g. the sun shining at a modelled angle).
All our data is processed into regular grids of either size 32 × 32 or 64 × 64. For the exact processing of each of the datasets, please refer to our code. B: Experimental Setting, Model Architectures and Compute B.I: Generative Spatial Modeling Setup: Our main experimental findings, the MMD scores displayed in Table 1 and Table 2 (main paper), are obtained from training generative models on 60% of the data, holding out 20% of the data for validation and model selection, and 20% for computing the displayed test scores. This setting is used for all four experimental datasets. For each dataset, model architecture and auxiliary task setting, we train 10 GANs with different random initializations, in each cycle saving the best generator according to tests on validation data. We then choose the best out of the 10 trained generators (again according to the validation score) to compute test scores.
Model Architecture and Optimization: Here we briefly describe the model architectures of the different generative models used in the experiments working with 32 × 32 inputs (the models for the 64 × 64 input are adapted to fit the larger input). For the implementation of these models, please refer to our code.
The Vanilla GAN architecture used consists of a Generator with four hidden linear layers, supported by Leaky ReLU and 1 Batch-Norm layers. The Discriminator has two hidden linear layers supported by Leaky ReLU layers and one linear task-specific layer.
The DCGAN architecture used consists of a Generator with a linear initialization layer, followed by three hidden (de-)convolutional layers, supported by ReLU and 2 BatchNorm layers. The Discriminator contains two convolutinal layers supported by Leaky ReLU and 2 BatchNorm layers, followed by one task-specific convolutional layer with a final linear transformation. For more information on DCGAN, please refer to the original publication [38].
The EDGAN architecture used consists of an Encoder-Decoder Generator, where the Encoder contains three convolutional layers, supported by Leaky ReLU and 2 BatchNorm layers and the Decoder contains three (de-)convolutional layers supported by ReLU layers. The Discriminator has five hidden convolutional layers suppoerted by Leaky ReLU and 2 BatchNorm layers, followed by a last, taskspecific convolutional layer. For more information on the EDGAN architecture, please refer to [57], the study which motivated the use of this benchmark.
Model Training: All models are trained using the binary cross entropy criterion to compute losses. Optimization through backpropagation is conducted using the Adam algorithm with a learning rate of 0.001 and values of [0.5, 0.999]. Experiments with the Toy dataset run for 40 epochs, with the PetrelGrid dataset for 500 epochs, with the DEM dataset for 100 epochs and with the TreeCanopy dataset for 100 epochs. All training is conducted on GPUs provided via Google Colab, which includes Tesla K80, Tesla T4 and Tesla P100 GPUs. The model training times do not exceed 30 minutes at the longest.
Evaluation: To evaluate our models, we generate synthetic data from the different types of GANs and compare how faithful the generated samples are compared to true samples. To assess model quality, we use the Maximum Mean Discrepancy (MMD) metric