Modeling spatial interactions is used to understand and quantify the level of interaction between different locations. It implies carrying out the three following steps:

**Conceptual formalization**of entities (e.g. individuals, spatial units) and relationships (e.g. work, power, goods).- Corresponent
**mathematical formalization**: \(T_{ij} = S_i D_j / d_{ij}^2\) the trade flow between units*i*and*j*(\(T_{ij}\)) is proportional to the supply of unit*i*(\(S_i\)) and the demand of unit*j*(\(D_j\)) and inversely proportional to the squared distance between both units (\(d_{ij}^2\)). **Informatic formalization**of the model: language, algorithms, etc.

The most ancient and common spatial interactions model is of the
gravity form (like below). It has roots in the late 19st century (Casey,
Ravenstein) and has been used in several fields (geography, economy,
demography) to model a high variety of flows (commuting, trade,
migrations). A brief presentation can be found in the Geography
of transport systems. For a detailed one see, among many others, the
book of Fotheringham and O’Kelly, *Spatial Interaction Models:
Formulations and Applications* (Kluwer Academic Publishers,
1989).

There are two main ways of modeling spatial interactions: the first
one focuses on links between places (flows), the second one focuses on
places and their influence at a distance. Following the physics
metaphor, the flow may be seen as the **gravitational
force** between two masses, the place influence as the
**gravitational potential**. The
`SpatialPosition`

package, as its name suggests, proposes an
implementation of place-based models: Stewart, Reilly and Huff
models.

The Stewart model [1] is also known in the literature, depending on
the discipline, as **potential access**,
**gravitational potential** or **gravitational
accessibility**. The concept was developped in the 1940s by
physicist John Q. Stewart from an analogy to the gravity model. In his
seminal work on the catchment areas of American universities, Stewart
computes **potentials of population**. This potential is
defined as a stock of population weighted by distance:

\[ A_i = \sum_{j=1}^n O_j f(d_{ij}) \]

The terms of this equation can be interpreted in a potential approach or an accessibility approach (in brackets):

- \(A_i\) is the potential at \(i\) (the accessibility)
- \(O_j\) is the stock of population at \(j\) (the number of opportunities)
- \(f(d_{ij})\) is a negative function of the distance between \(i\) and \(j\), mainly of the power or the exponential form.

The computation of potentials could be considered as a
**spatial interpolation method** such as inverse distance
weighted interpolation (IDW) or kernel density estimator. These models
aim to estimate unknown values of non-observed points from known values
given by measure points. Cartographically speaking, they are often used
to get a continuous surface from a set of discrete points. However, we
argue that the Stewart model is mainly a spatial interaction modeling
approach, with a possible secondary use for spatial interpolation.

The Reilly [2] and the Huff [3] models draw catchment areas. They have been widely used to study the location of economic activities and are part of the geomarketing toolbox. These models aim to quantify the attraction force on a location \(i\) generated by an opportunity \(j\):

\[ A_{ij} = O_j f(d_{ij}) \]

The attraction force between \(i\) and \(j\) (\(A_{ij}\)) is proportional to the mass of \(j\) (\(O_j\)) and inversely proportional to the distance between both locations (\(d_{ij}\)). Let \(i\) be the location of a consumer and \(j\) the locations of a set of shopping malls. We compute all the attraction forces sustained by \(i\) and we assign \(i\) to the mall with the strongest attraction force. Computed on a regular fine grid of points \(i\) the model produces a set of deterministic catchment areas.

The Huff model is a relative version of the Reilly model: for a given point \(i\) and a given attraction center \(j\), the attraction force of \(j\) is divided by the sum of all the possible attraction centers that affect \(i\). The result my be understood as the probability to choose \(j\) among the set of possible destinations. Computed on a regular fine grid of points \(i\) the model produces a raster representing probable catchment areas.

Modeling spatial interactions means quantifying the distance friction
or impedance. The role of the distance can be interpreted as a
disincentive to access desired destinations or opportunities (e.g. jobs,
shops). At the very place of the opportunity, the interaction function
equals 1, meaning that the potential access is 100%. Far away from the
opportunity, the interaction function tends to 0, meaning that the
potential access is 0 %. The **span** is defined as the
value where the interaction function falls to 0.5 (50%). From the
individuals’ point of vue, this function may be seen as a degree of
availability of a given opportunity. From the opportunity’s point of vue
(a store for example), the interaction function may be seen as a
decreasing catchment area: there is a maximal attraction close to the
opportunity and this attraction decreases progressively through
distance.

Tha example dataset consists in two spatial objects:

`paris`

(polygons): 1 spatial unit representing the Paris’ perimeter`hospital`

(points): 18 points representing the public hospitals with their capacity (number of beds)

First we want to compute a global accessibility to hospitals. This computation could lead the residential choice of a highly hypocondriac individual searching for the optimal location with the best global access to hospitals. It could also lead an investment choice of the local authority searching for the optimal location to improve equity of access to health services.

The set of known points are the hospitals and the weighting variable
– i.e. the mass of the points – is their capacity
(`capacity`

). We want to create a raster of values and to
clip it with the perimeter of the study area so we provide clipping mask
(`paris`

). The function type, the span and the exponent
should reflect research hypothesis or known behavior of the individuals:
e.g. in a dense urban area such as the center of Paris, an average
consumer wouldn’t go farther than 500 meters to reach a small retail
shop and the friction of the distance is very high (exponential form
with a large beta coefficient). Here we set the span at 1000 meters, the
impedance function is exponential with beta = 3.

`library(SpatialPosition)`

```
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
## (status 2 uses the sf package in place of rgdal)
```

```
## Functions related to Stewart's potential are deprecated.
## Please use the `potential` package instead.
## https://riatelab.github.io/potential/
```

`library(sf)`

`## Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE`

```
data(hospital)
# Compute potentials (accessibility)
<- stewart(
potentials knownpts = hospital,
varname = "capacity",
typefct = "exponential",
span = 1000,
beta = 3,
resolution = 50,
mask = paris,
returnclass = "sf"
)
```

```
## Warning: stewart(), rasterStewart(), plotStewart(), quickStewart(), isopoly(),
## mcStewart() and smoothy() are deprecated. Use 'potential' package for all
## potential-related functions.
```

```
<- isopoly(x = potentials, mask = paris, returnclass = "sf")
isopotentials
<- paste0(round(isopotentials$min,0),' to ',
lab round(isopotentials$max,0))
par(mar = c(4,2,2,1))
plot(st_geometry(isopotentials), col = heat.colors(8))
legend(x = "topright", legend = lab, fill = heat.colors(8),
cex = 0.7, title = "Potentials")
mtext("Global Accessibility to Public Hospitals", side = 3,cex = 1.5)
mtext(text = "Potential nb. of beds
distance function: exponential, span = 1 km, beta = 3",
side = 1, line = 1)
```

We now want to draw the catchment areas corresponing to the public
hospitals. These areas could be used as a zoning tool, in order to set
the residential zones with its referent hospital. As for the global
accessibility, three steps are needed: computing the values
(`reilly`

), creating the raster (`rasterReilly`

),
plotting the raster (`plotReilly`

).

`row.names(hospital)`

```
## [1] "208" "207" "205" "204" "201" "200" "199" "198" "197" "195" "194" "193"
## [13] "192" "191" "190" "187" "186" "185"
```

```
<- reilly(knownpts = hospital, varname = "capacity",
catchReilly typefct = "exponential", span = 750, beta = 2,
resolution = 50, mask = paris, returnclass = "sf")
# Create a raster
<- rasterReilly(x = catchReilly, mask = paris)
rasterCatch
par(mar = c(4,2,2,1))
# Plot the raster and add the points
plotReilly(x = rasterCatch)
plot(st_geometry(hospital), pch = 20, add = TRUE)
mtext("Catchment Areas of Public Hospitals", side = 3,cex = 1.5)
mtext(text = "distance function: exponential, span = 0.75 km, beta = 2",
side = 1, line = 0)
```

We now want to draw the probabilistic catchment areas corresponing to
the public hospitals. These areas show the probability for a patient to
go to one hospital or to another. It is used in geomarketing to detect
the fuzzy areas where there are competing shopping malls and where the
brand should focus its advertising. The same three steps are needed:
computing the values (`huff`

), creating the raster
(`rasterHuff`

) and plotting the raster
(`plotHuff`

).

```
<- huff(knownpts = hospital, varname = "capacity",
catchHuff typefct = "exponential", span = 750, beta = 2,
resolution = 50, mask = paris, returnclass = "sf")
# Create a raster
<- rasterHuff(x = catchHuff, mask = paris)
rasterCatch
# Plot the raster and add the points
par(mar = c(4,2,2,1))
plotHuff(x = rasterCatch)
plot(st_geometry(hospital), pch = 20, col = "red", add = TRUE)
mtext("Probabilistic Catchment Areas \nof Public Hospitals",
side = 3,cex = 1.5, line=-1.5)
mtext(text = "distance function: exponential, span = 0.75 km, beta = 2",
side = 1, line = 0)
```

[1] STEWART J.Q. (1942) “Measure of the influence of a population at
a distance”, *Sociometry*, 5(1): 63-71.

[2] REILLY W.J. (1931) *The law of retail gravitation*, W. J.
Reilly, New York.

[3] HUFF D. (1964) “Defining and Estimating a Trading Area”, *Journal
of Marketing*, 28: 34-38.