Traditional spatial modeling approaches assume that data are second-order stationary, which is rarely true over large geographical areas. A simple way to model nonstationary data is to partition the space and build models for each region in the partition. This has the side effect of creating discontinuities in the prediction surface at region borders. The regional border smoothing approach ensures continuous predictions by using a weighted average of predictions from regional models. The R package *remap* is an implementation of regional border smoothing that builds a collection of spatial models. Special consideration is given to distance calculations that make *remap* package scalable to large problems. Using the *remap* package, as opposed to global spatial models, results in improved prediction accuracy on test data. These accuracy improvements, coupled with their computational feasibility, illustrate the efficacy of the *remap* approach to modeling nonstationary data.

When observations exhibit spatial autocorrelation, geographic location can be leveraged to improve predictions of the response variable by considering responses in nearby observations. Typical spatial statistical models assume that the covariance between two observations can be modeled as a function of location difference, i.e., the relationship needs to be second-order stationary. For sufficiently large distances, the stationarity assumption often fails. Even when external drift or secondary variables are used, continental scale models may still not be second-order stationary.

The simplest way to use traditional spatial modeling with nonstationary data is to partition the global region into smaller sub-regions that are locally stationary and create separate models for each region. The naïve implementation of this approach leads to noncontinuous predictions at the borders of each region. Previous attempts to smooth out discontinuities at region boundaries often involved taking weighted averages of local regional model output where the weights of each local model prediction or covariance structure are a function of the distance between a new observation and the centers of each region (Fuentes 2001; Fuentes and Smith 2001; Gosoniu et al. 2006; Gosoniu et al. 2009; Konomi et al. 2014). While the center based approach may be appropriate for symmetrical regions, it is likely not appropriate for oddly shaped or disjoint regions not well represented by their centers. In some cases, a region may not even contain its center. The center-based approach also fails to respect major geographic features that can cause sharp changes in response variables over very short distances.

The regional border smoothing approach described in this article is a
novel method that uses a weighted average of predictions from local
models, but gives weight to regional models based on the nearest
distance to the *border* of each region rather than distance to the
center of each region. With this method, regions may be chosen that more
naturally reflect local climate and topography with smoothing only
occurring near the borders of each region. Figure 1 is
an example of the regional border smoothing approach applied to three
regions with different spatial models.

The R package *remap* is a
General Public License implementation of this regional border smoothing
method that scales well to large problems (Wagstaff 2022). Using *remap*,
regional border smoothing is applied to two different modeling problems.
The first problem is a national set of 50-year ground snow loads and the
second problem is modeling April 1st snow water content for Utah. Mapped
values from both projects show improvement in model accuracy using the
regional border smoothing approach over global models for a variety of
spatial modeling approaches.

The body of this article will proceed with a description of the regional
border smoothing approach. This will be followed by an illustration of
the available functions and tools in *remap* as well as two
demonstrations of the software on a state and national-level data set.
These examples show the utility of *remap* in producing smooth estimates
when applied to spatial modeling problems over large geographical areas
with irregularly shaped partitions.

There are many proposed methods to model nonstationary data using locally stationary models. Haas (1990b,a) describes a moving window approach where only data within a pre-specified bounding box are used to fit a local dependence structure and then make predictions. This process is computationally intensive and may not result in continuous predictions.

Fuentes (2001) and Fuentes and Smith (2001) propose a method for
nonstationary problems where a global covariance structure that changes
continuously as a function of location is used in a Gaussian process
model. The data are first partitioned into locally stationary regions
and a global covariance structure is calculated by taking a weighted
average of regional covariance structures. Weights are based on the
distance from the prediction location to a point in each region, usually
the center. Regions are either given *a priori* or by using subgrids
chosen using the Bayesian information criterion.

Applications of local partition modeling approaches include Kim et al. (2005), who describe a method to deal with sudden changes in spatial covariance structure that occur between layers of rock strata. The spatial domain is partitioned into independent regions using Voronoi Tessellations (Green and Sibson 1978), with each region fit using an independent Gaussian model. The resulting global model has sharp changes at the borders of each region, which was desirable given the context of the problem. Konomi et al. (2014) illustrate a decision tree based method for partitioning the spatial domain when modeling global Ozone levels. Heaton et al. (2017) use a hierarchical clustering method to partition the spatial domain for temperature data in Houston, TX. The hierarchical clustering method has the benefit of creating a partition that more naturally follows changes in the covariance structure rather than partitioning the space into symmetrical blocks or spheres.

Gosoniu et al. (2006; Gosoniu et al. 2009) provide an additional
application of local partitioning models mapping malaria risk using *a
priori* partitions of West Africa. The 2006 study uses three large
rectangular regions and the 2009 study uses agro-ecological regions to
partition the spatial domain. In these studies, spatial random effects
are modeled as a weighted sum of regional stationary effects based on
the distance to region centroids. The authors note problems with sudden
changes at region borders as a result of using region centroids for the
weighted sum of effects.

Most of the methods discussed so far create a global covariance structure from local covariance structures of Gaussian process models (Fuentes 2001; Fuentes and Smith 2001; Kim et al. 2005; Konomi et al. 2014; Heaton et al. 2017). Each method has a different way of smoothing regional transitions that are specifically tied to their methodology. Many of the methods use a Bayesian framework that requires computationally expensive Markov chain Monte Carlo simulations (Kim et al. 2005; Gosoniu et al. 2006; Gosoniu et al. 2009; Konomi et al. 2014; Heaton et al. 2017). Many of the techniques are only applied to purely spatial data rather than multivariate data (Kim et al. 2005; Konomi et al. 2014; Heaton et al. 2017). There is a lack of methodology that allow for smooth transitions between partitions and works for multiple modeling techniques.

The regional border smoothing method described in this article can be thought of as stitching together images to form a larger image with no sharp changes. The approach is similar to those used to combine black and white images from microscopes into a larger image (Thévenaz and Unser 2007). Individual microscope images are aligned and the overlapping regions are smoothed by taking a weighted average of the overlapping pixel brightness. The weights are based on the distance from the pixel to the outer edge of each image with more weight being applied to pixels closer to the center of an image.

The process described in this article provides a simple way of combining regional model predictions to form a continuous global prediction surface. The method can be applied to problems that are not strictly spatial. For example, Osborne and Suárez-Seoane (2002) show that building models for partitioned space can improve the accuracy of large scale species distribution models. The regional border smoothing method works for any model that produces continuous predictions.

The *remap* package is
available on the Comprehensive R Archive Network (see
https://cran.r-project.org/web/packages/remap/index.html) with the
most current version available at
https://github.com/jadonwagstaff/remap. The code and data used to
generate the results in this article are available as supplementary
materials.

The figures and tables in this article are made using the R programming
language. Figures are created with the
*tidyverse*
(Wickham et al. 2019),
*gridExtra*
(Auguie 2017), *cowplot*
(Wilke 2020), and *maps*
(Becker et al. 2021) packages. The
*sf* (Pebesma 2018),
*nngeo* (Dorman 2022), and
*raster* (Hijmans 2022)
packages are used to manipulate spatial data. Kriging models are built
with *automap*
(Hiemstra et al. 2009) and *gstat*
(Pebesma 2004; Gräler et al. 2016) and generalized additive models are built
with *mgcv* (Wood 2011). A
docker container with these packages installed and all code is available
at https://hub.docker.com/r/jadonwagstaff/remap_manuscript_code.

This research was funded by the American Society of Civil Engineers and the Structural Engineering Institute (award number 202827). The authors have no affiliation with any organization with a direct or indirect financial interest in the subject matter discussed in the manuscript. This article is based on the first author’s master’s thesis (Wagstaff 2021).

Regional border smoothing is the process of using regional models to make predictions that are globally continuous. Regional border smoothing may be used on any spatially referenced data (\(\pmb{X}\)) in conjunction with a modeling approach and a finite set of regions \(\mathcal{R}_p\) where each \(\mathcal{R}_{p_i} \in \mathcal{R}_p\), \(i = 1 \dots m\), is a closed set of points contained in the region of interest \(\mathcal{R}\). The intention is that \(\mathcal{R}_p\) is a set of non-overlapping polygons with shared borders and \(\bigcup\limits_{i=1}^m \mathcal{R}_{p_i} = \mathcal{R}\), but these are not necessary conditions. While \(\mathcal{R}_p\) does not meet the strict definition of a partition, \(\mathcal{R}_p\) will be referred to as a partition throughout this article.

The border smoothing approach described in this article was originally designed for regression-based models, but can be used with any modeling approach that produces a continuous response. This technically includes classification techniques that make continuous probability predictions prior to classifying based on probability threshold, such as logistic regression. Presumably, the predictions from a chosen modeling approach will result in continuous predictions as a function of location.

Modeling regions are defined by the borders of \(\mathcal{R}_p\). The data contained by a modeling region \(\mathcal{R}_{p_i}\) are used to inform a distinct regression model (\(f_{p_i}(\pmb{X})\)) for that region. In some cases it may be desirable to include observations near each \(\mathcal{R}_{p_i}\) when building regional regression models. In these cases, data within \(\mathcal{R}_{p_i}\) and within a buffer zone around each modeling region are used to build each \(f_{p_i}(\pmb{X})\) (Figure 2). Using points within a buffer zone that extends beyond the region boundaries avoids edge extrapolation when using a regional model for interpolation in the smoothing zones of neighboring regions.

Simply making predictions within each \(\mathcal{R}_{p_i}\) using each corresponding \(f_{p_i}(\pmb{X})\) results in noncontinuous predictions at region boundaries. Regional border smoothing results in continuous predictions for the entire space by taking a weighted average of the predictions provided by each \(f_{p_i}(\pmb{X})\). The smoothed global prediction surface is continuous, but not necessarily differentiable.

Let \(\hat{y}_{1}, \hat{y}_{2}, ..., \hat{y}_{m}\) represent predictions from regional models \(f_{p_1}(\pmb{x}^*), f_{p_2}(\pmb{x}^*), ..., f_{p_m}(\pmb{x}^*)\) for location of interest \(\pmb{x}^*\) where \(\pmb{x}^*\) represents both spatial and non-spatial information. A final prediction \(\hat{y}'\) is calculated by the weighted average of each \(\hat{y}\), i.e.,

\[\label{eq:smooth} \hat{y}' = \frac{\sum_{i=1}^m w(d_i|S)\hat{y}_{i}}{\sum_{i=1}^m w(d_i|S)}. \tag{1}\]

The weights are calculated based on the smallest great-circle distances (\(d_{1}, d_{2}, ..., d_{m}\)) between the location of \(\pmb{x}^*\) and the boundaries of \(\mathcal{R}_{p_1}, \mathcal{R}_{p_2}, ..., \mathcal{R}_{p_m}\). If \(\pmb{x}^*\) is located within a region, the distance between \(\pmb{x}^*\) and that region is \(0\). The weight \(w(d_i|S)\) given to \(\hat{y}_{i}\) is non-zero when \(d_{i}\) is within some threshold \(S\), i.e.,

\[\label{eq:weight} w(d_i | S) = \begin{cases} \left( \frac{S - d_{i}}{S} \right)^2 & d_{i} \leq S\\ 0 & d_{i} > S. \end{cases} \tag{2}\]

Consider now the special case when a prediction is made in region \(j\) and \(d_{i\ne j} > S\), then \(w(d_{i \ne j}|S) = 0\) and Equation (1) reduces to \(\hat{y}' = \hat{y}_{j}\). As prediction location in region \(j\) approach \(\mathcal{R}_{p_k}\), then the weight of \(\hat{y}_{k}\) increases gradually and \(\hat{y}' = \frac{ \hat{y}_{j} + \left(\left(S - d_{k}\right)/S\right)^2 \hat{y}_{k}}{1 + \left(\left(S - d_{k}\right)/s\right)^2}\). At the border of regions \(j\) and \(k\), \(\hat{y} = \left(\hat{y}_{j} + \hat{y}_{k}\right) / 2\). Finally, as prediction locations progress into \(\mathcal{R}_{p_k}\), weights for \(y_{i_j}\) decrease gradually to zero and \(\hat{y}' = \frac{\left(\left(S - d_{j}\right)/S\right)^2 \hat{y}_{j} + \hat{y}_{k}}{\left(\left(S - d_{j}\right)/S\right)^2 + 1}\). All locations within \(S\) of a region are referred to as the smoothing zone for that region.

The smoothing approach described in this paper represents a spatially
weighted ensemble of model predictions along region boundaries, with
each model (possibly) estimating prediction standard error (SE) or
variance. There is no consensus in the literature on how SE should be
calculated for ensemble model predictions. The methods for SE
calculation that do exist tend to be specific to model type (e.g.
(Wager et al. 2014)) or averaging approach (e.g.
(Hoeting et al. 1999)). For this reason, the method for estimating SE
in *remap* relies on the general properties of variance for the
summation of random variables.

Suppose that each \(\hat{y}_i\) is an unbiased estimate from of \(f_{p_i}(\pmb{x}^*)\) with the SE of \(\hat{y}_i\) represented as \(\hat{\sigma}_i\). Then the combined model SE for \(\hat{y}'\) can be represented as \[\hat{\sigma}' = \sqrt{\frac{\sum_{i=1}^m w(d_i|S)^2\hat{\sigma}_{i}^2 + \sum_{i=1}^m\sum_{j<i} 2w(d_i|S)w(d_j|S)\rho_{ij}\hat{\sigma}_{i}\hat{\sigma}_{j}}{\left(\sum_i w(d_i|S)\right)^2}},\] where \(\rho_{ij}\) represents the correlation between predictions for model \(i\) and \(j\). One primary difficulty with estimating \(\hat{\sigma}'\) in the smoothing zones is a lack of information regarding \(\rho_{ij}\). It is highly likely that SE estimates are correlated for models in adjacent regions. In the absence of a computationally efficient and theoretically robust estimate for \(\rho_{ij}\) the Cauchy–Schwarz inequality can be used to provide an upper bound on \(\hat{\sigma}'\) given as \[\label{eq:cs_in} \hat{\sigma}' \le \sqrt{\frac{\sum_{i=1}^m w(d_i|S)^2\hat{\sigma}_{i}^2 + \sum_{i=1}^m\sum_{j<i} 2w(d_i|S)w(d_j|S)\hat{\sigma}_{i}\hat{\sigma}_{j}}{\left(\sum_i w(d_i|S)\right)^2}}, \tag{3}\] where \(\rho_{ij}\) has been replaced with a constant value of one.

The *remap* package provides a method for estimating the upper bound of
the prediction SE using the inequality specified in ((3)).
It is important to note that it may not always be appropriate to combine
model SE estimates in this way. This is particularly true for models
that ignore spatial autocorrelation. It is also important to remember
that the primary focus of *remap* is to provide a practical approach for
smoothing model predictions, and not to preserve the theoretical
integrity of the SE estimates. It is up to the end user to determine
that the SE calculations are valid and appropriate for their modeling
purposes.

The “smoothing” described throughout this article refers to smoothing in the colloquial sense. A continuous prediction surface is created with a steady transition between regions. The prediction surface is not always differentiable. If a region is not convex, the rate of change in the distance to a region can shift suddenly at locations that are equidistant to different parts of the region. Figure 3 shows an example of a prediction surface near a non-convex region where the prediction surface is not differentiably smooth.

Since each region is a closed set of points, \(d_i\) is continuous. It therefore follows that the weight function \(w(d_i|S)\) described in Equation (2) is continuous as \(\lim_{x\to0^{+}} w(x|S) = w(0|S) = 1\) and \(\lim_{x\to S^{-}} w(x|S) = w(S|S) = 0\) for \(x \in [0, \infty)\) and \(S > 0\). Since the right side of Equation (1) is the multiplication and addition of continuous functions, it is guaranteed to have \(\max\left(w(d_i|S)\right) > 0\) as long as \(\pmb{x}^*_s\) is within \(S\) units of any \(\mathcal{R}_p\) and it follows that Equation \(\hat{y}'\) is also continuous. Given continuous predictions as a function of location for each region, the regional border smoothing method is guaranteed to be continuous for any location within the smoothing zone of at least one region. Once all \(d_i > S\), then the denominator of \(y'\) is equal to zero, which means \(y'\) is no longer well-defined. As long as all of the data \(\pmb{X}\) are located within \(\mathcal{R}_p\), any prediction location outside of \(\mathcal{R}_p\) is spatial extrapolation that is generally discouraged.

The R package *remap* is an implementation of the regional border
smoothing approach to spatial modeling. The function `remap`

creates a
set of regionalized models given:

- A set of observations as spatially projected or geographic points.
- A set of regions as spatially projected or geographic polygons.
- A desired buffer zone distance.
- A modeling function to apply to observations in each region.

Predictions can be made on new observations given a regionalized model
and the smoothing parameter `smooth`

used for weighted averages in
smoothing zones (variable \(S\) in Equation (2)). Detailed
descriptions of function parameters can be found in the package
documentation. Some working examples using the *remap* package are
provided via the vignette which accompanies the package. The development
version of the code for *remap* can be found at
https://github.com/jadonwagstaff/remap.

The weights for regional predictions require the distances from all
observations to the boundary of each region. The modeling process of the
`remap`

function also requires the distances from all observations to
each region to assign the correct observations to each regional model.
The process of fine-tuning a model can result in recalculating these
distances many times. To avoid recalculating distances, the function
`redist`

is included in the *remap* package to pre-compute the distances
from a given set of observations to a given set of regions. These
pre-computed distances can be supplied as a parameter to the `remap`

function to reduce computational time while fine-tuning models.

Calculating distances in the *remap* package takes advantage of tools
already available for spatial analysis in the R package *sf*
(Pebesma 2018). The function `sf::st_distance`

is used to find either
Euclidean or great-circle distances depending on the spatial projection
of the locations and regions. The `sf::st_distance`

function uses the
*S2* geometry library (Google 2022) when calculating great-circle
distances. The *S2* geometry library is designed specifically to
efficiently compute distances between nearby objects given large sets of
geographic data. Regardless, calculating distances can still take a lot
of computational time if the regions are complex and/or there are a lot
of observations.

The naïve approach to regional border smoothing is to find all distances between each location and each region. This approach may be necessary during the modeling process if any region does not contain a required minimum number of observations (detailed in the implementation considerations section). If predictions are being made using new observations, then distances do not need to be calculated between every observation and every region.

Because the weight for predictions in different regions is zero when the
distance to those regions is greater than the smoothing parameter
(Equation (2)), distances do not need to be calculated
between *every* region and *every* prediction location. To determine
which observations require distance calculations, an approximate polygon
is constructed that encompasses the original region plus the smoothing
zone around the region. The function `sf::st_within`

is used to
determine which observations are within the new polygon and
equivalently, which observations are within the smoothing zone of the
original region. Observations in the approximate polygon become
candidates for precise distance calculations.

An approximate polygon that contains a region and the region’s smoothing
zone is created using the `sf::st_buffer`

function. For geographic
coordinates, `sf::st_buffer`

requires a buffer value in degrees. Since
the great-circle distance between a unit degree of longitude changes
with latitude, Equation (4) is used to find \(c\): the
shortest distance (measured in kilometers) required to move one degree
longitude at the observation in the data set that is nearest to a
geographical pole. The set \(\pmb{\lambda}\) is the set of Latitude values
of the observations and 6380 km is the radius of the earth at a pole
(rounded down). A sufficient buffer value in degrees for `sf::st_buffer`

is calculated by dividing the smoothing zone length in km by \(c\). The
resulting buffered polygon contains all of the observations within the
smoothing zone of the original polygon,

\[\label{eq:km_per_degree} c = \frac{\pi * 6350}{180} * cos\left(\frac{\pi * \max\left(|\pmb{\lambda}|\right)}{180}\right). \tag{4}\]

Assuming each region contains a sufficient number of observations to fit
each regional model, the same process used for reducing the number of
distance calculations for model predictions can be used to reduce the
number of distance calculations required for model training. The
function `redist`

is able to restrict distance calculations to only
observations within a certain buffer zone or smoothing zone of each
region using the `max_dist`

parameter. This eliminates the need to
calculate distances to points with a known weight of zero when
smoothing.

Reliable regression results depend upon sufficient sample sizes, which
differ based on the variability of the response and dimensionality of
the inputs. A strict (and obvious) minimum sample size for the *remap*
function is one observation per region, though this threshold will
rarely, if ever, ensure reasonable results. An additional parameter
`min_n`

is included in the `remap`

function to specify a minimum number
of observations be used to build each model. If a region and the buffer
around that region do not contain the minimum number of observations
specified, the `min_n`

observations closest to the boundaries of the
region are used to train that region’s model.

The *remap* package has the ability to make predictions outside of all
modeling regions. If the prediction location is within the smoothing
zone of a region, then Equation (1) applies and the nearest
region will have the most weight. If the prediction location is outside
the smoothing zone of all regions, then the model from the closest
region is used to make a prediction. This may result in non-continuous
transitions in predictions when the closest region changes across
geographic space, but this is only possible at locations outside of the
smoothing zone of all regions. As a general rule, extrapolation of
predictions to a location beyond those represented in the input data is
not recommended. See Figure 4 for a visual depiction of how
predictions behave outside of all modeling regions.

The ability to predict outside of all regions means that the *remap*
package can make smooth predictions across small gaps at polygon
borders, provided the gaps are no larger than two times the smoothing
zone distance. This encourages the use of `rgeos::gSimplify`

(Bivand and Rundel 2021) or `sf::st_simplify`

(Pebesma 2018) to simplify complex
polygons and ease the computational burden associated with distance
calculations, even if those simplifications slightly compromise the
topology of the original geometry.

Snow loads are obtained from weather stations throughout the United
States in the National Oceanic and Atmospheric Administration’s (NOAA)
Global Historical Climatology Network (Menne et al. 2012). Snow loads are
either measured directly, or estimated from measured snow depth using a
depth to load conversion model. Engineers have historically used
estimates of 50-year ground snow loads when designing structures. The
50-year snow load is traditionally obtained by fitting a probability
distribution to yearly maximum snow load measurements and extracting the
98^{th} percentile. A recent effort by the American Society of Civil
Engineers has resulted in a new set of 50-year ground snow loads at 7964
measurement locations (Bean et al. 2021).

Building design requirements call for continuous maps that estimate design loads between measurement locations. This has historically been accomplished using various mapping techniques (Tobiasson et al. 2002; Liel et al. 2017; Bean et al. 2019). The problem is that the relationship between predictor variables and snow load can change drastically on a continental scale. For example, the typical loads at 1500 meter elevation in the Rocky Mountains of Montana are much lower than loads at 1500 meters in the Cascade Mountains of Washington State (Figure 5). Commonly used geospatial models for mapping are not well suited for the nonstationary nature of this problem.

Snow loads are typically assumed to share a log-linear relationship with
elevation. This relationship can be modeled directly using ordinary
least squares (OLS), which ignores any potential spatial dependencies
among the observations. A generalized additive model (GAM) built with
the *mgcv* R package (Wood 2011) characterizes the log of 50-year loads
as a function of elevation and a spatial smoother called splines on the
sphere (Wood 2003). Universal kriging interpolates values using a
Gaussian process model after accounting for the log-linear trend in
elevation. Variograms are individually fit within each region using the
*automap* R package (Hiemstra et al. 2009).

Geographic regions defined by the US Environmental Protection Agency
(EPA) define regions with similar ecology and climate called eco-regions
(Commission for Environmental Cooperation 1997). The eco-regions provide a natural partition of the
conterminous United States and give no regard to political boundaries.
Snow load can be modeled using observations within each eco-region where
the relationship between predictor variables and the response are more
consistent on a local level. The *remap* package facilitates modeling in
separate eco-regions and creates a smooth model on the national scale.
There are 86 eco-regions that fall within the conterminous United
States, so 86 separate models are built using the `remap`

function with
a buffer zone of 50 km and a `min_n`

of 150 observations. The regional
models are smoothed to a single continuous model using a `smooth`

parameter of 25 km.

The following code demonstrates how `remap`

is used to build these
models. The script and data for the following examples are provided as
supplementary materials with this article. The example data include
three `"sf"`

objects: a spatial points data frame with 50-year snow
loads at locations within the US and Canada (`loads`

), a spatial
polygons data frame of eco-regions (`eco3`

), a spatial polygon of the
conterminous United States (`cont_us`

), and a `"raster"`

object with
elevations comprising a 0.8 km grid over the conterminous United States
(`grd`

).

```
library(tidyverse)
library(sf)
library(raster)
library(mgcv)
library(automap)
library(gstat)
library(remap)
load("wagstaff-bean.RData")
```

The eco-regions have more complex borders than is necessary for this problem, so they can be simplified to the same extent as shown in Figure 6.

```
<- eco3 %>%
eco3_simp ::st_simplify(dTolerance = 10000) %>%
sf::filter(!sf::st_is_empty(.)) %>%
dplyr::st_cast("MULTIPOLYGON") sf
```