Title: | Shortest Paths Between Points in Grids |
---|---|
Description: | Shortest paths between points in grids. Optional barriers and custom transition functions. Applications regarding planet Earth, as well as generally spheres and planes. Optimized for computational performance, customizability, and user friendliness. Graph-theoretical implementation tailored to gridded data. Currently focused on Dijkstra's (1959) <doi:10.1007/BF01386390> algorithm. Future updates broaden the scope to other least cost path algorithms and to centrality measures. |
Authors: | Christian Düben [aut, cre] |
Maintainer: | Christian Düben <[email protected]> |
License: | MIT + file LICENSE |
Version: | 1.1.3 |
Built: | 2024-11-09 05:25:23 UTC |
Source: | https://github.com/cdueben/spaths |
This function draws random unprojected (lonlat) locations.
rnd_locations( nobs, xmin = -180, xmax = 180, ymin = -90, ymax = 90, output_type = c("data.table", "data.frame", "SpatVector") )
rnd_locations( nobs, xmin = -180, xmax = 180, ymin = -90, ymax = 90, output_type = c("data.table", "data.frame", "SpatVector") )
nobs |
Number of observations |
xmin |
Minimum longitude |
xmax |
Maximum longitude |
ymin |
Minimum latitude |
ymax |
Maximum latitude |
output_type |
type of output object. Either |
By default, the function draws a global sample of random locations. You can restrict it to a certain region via specifying xmin
,
xmax
, ymin
, and ymax
. The function draws from a uniform spatial distribution that assumes the planet to be a perfect sphere. The
spherical assumption is common in GIS functions, but deviates from Earth's exact shape.
Returns a data.table, data.frame, or SpatVector object of unprojected (lonlat) points.
rnd_locations(1000)
rnd_locations(1000)
The shortest paths and/ or distance between locations in a grid according to Dijkstra's (1959) algorithm.
shortest_paths( rst, origins, destinations = NULL, output = c("distances", "lines", "both"), output_class = NULL, origin_names = NULL, destination_names = NULL, pairwise = FALSE, contiguity = c("queen", "rook"), spherical = TRUE, radius = 6371010, extent = NULL, dist_comp = c("spaths", "terra"), tr_fun = NULL, v_matrix = FALSE, tr_directed = TRUE, pre = TRUE, early_stopping = TRUE, bidirectional = FALSE, update_rst = NULL, touches = TRUE, ncores = NULL, par_lvl = c("update_rst", "points"), show_progress = FALSE, bar_limit = 150L, path_type = c("int", "unsigned short int"), distance_type = c("double", "float", "int", "unsigned short int") )
shortest_paths( rst, origins, destinations = NULL, output = c("distances", "lines", "both"), output_class = NULL, origin_names = NULL, destination_names = NULL, pairwise = FALSE, contiguity = c("queen", "rook"), spherical = TRUE, radius = 6371010, extent = NULL, dist_comp = c("spaths", "terra"), tr_fun = NULL, v_matrix = FALSE, tr_directed = TRUE, pre = TRUE, early_stopping = TRUE, bidirectional = FALSE, update_rst = NULL, touches = TRUE, ncores = NULL, par_lvl = c("update_rst", "points"), show_progress = FALSE, bar_limit = 150L, path_type = c("int", "unsigned short int"), distance_type = c("double", "float", "int", "unsigned short int") )
rst |
SpatRaster (terra), RasterLayer (raster), matrix, or list of matrices object. RasterLayers are converted to SpatRasters. Pixels with non-NA
values in all layers mark the cells through which the algorithm may pass. The values of |
origins |
Origin points at which the shortest paths start. If |
destinations |
Destination points to which the shortest paths are derived. It defaults to |
output |
|
output_class |
Class of the returned object. With |
origin_names |
Character specifying the name of the column in the |
destination_names |
Character specifying the name of the column in the |
pairwise |
Logical specifying whether to compute pairwise paths, if |
contiguity |
|
spherical |
Logical specifying whether coordinates are unprojected, i.e. lonlat, and refer to a sphere, e.g. a planet, if |
radius |
Radius of the object, e.g. planet, if |
extent |
Vector of length 4, specifying the extent of |
dist_comp |
Method to compute distances between adjacent cells. |
tr_fun |
The transition function based on which to compute edge weights, i.e. the travel cost between adjacent grid cells. Defaults to the
geographic distance between pixel centroids. Permitted function parameter names are |
v_matrix |
Logical specifying whether to pass values to |
tr_directed |
Logical specifying whether |
pre |
Logical specifying whether to compute the distances between neighboring cells before executing the shortest paths algorithm in C++.
|
early_stopping |
Logical specifying whether to stop the shortest path algorithm once the target cells are reached. It defaults to |
bidirectional |
Logical specifying whether to produce paths or distances in both directions, if destinations are not specified and no directed
transition function is given. In that case, the distance and the path from point A to point B is the same as the distance and path from point B to point
A. |
update_rst |
Object updating |
touches |
Logical specifying the |
ncores |
An integer specifying the number of CPU cores to use. It defaults to the number of cores installed on the machine. A value of 1 induces a single-threaded process. |
par_lvl |
|
show_progress |
Logical specifying whether the function prints messages on progress. It defaults to |
bar_limit |
Integer specifying until up to how many paths or list elements of |
path_type |
Data type with which C++ stores cell numbers. |
distance_type |
Data type with which C++ stores distances. |
This function computes shortest paths and/ or distance between locations in a grid using Dijkstra's algorithm. Examples are a ship navigating around land masses or a hiker traversing mountains.
Let us explore the ship example to illustrate how shortest_paths
works. To compute shortest paths between ports around the world, you start with
a global SpatRaster, in which all land pixels are set to NA and all ocean pixels are set a non-NA value, e.g. 1. A SpatVector marks port locations as
points on water pixels. Passing these two objects to the parameters rst
and origins
respectively derives the shortest paths from each port
to all other ports conditional on ships solely traversing water pixels and returns the distances, i.e. lengths of these paths. If you are not interested
in the distances, but in the spatial lines themselves, set output
to "lines"
. If you want to obtain both, set it to "both"
.
In a different application, you do not want to compute the paths between all ports, but only the paths between the ports on the northern hemisphere and
the ports on the southern hemisphere, but not the paths between ports within the same hemisphere. To assess this, you split the ports into two
SpatVectors and pass them to origins
and destinations
respectively. What if you do not want to connect all orgins to all destinations? Set
pairwise
to TRUE
to connect the first origin just to the first destination, the second origin to the second destination, etc.
By default, the distance or transition cost between adjacent cells of the input grid is the geographic distance between the cells' centroids. What if
the boat is a sailing vessel that minimizes travel time conditional on wind speed, wind direction, and ocean currents? Construct a SpatRaster with
three layers containing information on the three variables respectively. Define a transition function that combines the three layers into a travel time
measure and pass the SpatRaster to rst
and this function to tr_fun
. tr_fun
makes this package very versatile. With custom
transition functions, you can take this software out of the geo-spatial context and e.g. apply it to biomedical research.
Going back to the ship routing example, consider that there are hurricanes in the Caribbean. Ships traveling from India to Australia do not care, but
ships traveling from Mexico to the Netherlands have to go around the storm and must not take the shortest path through the hurricane. You have ten
SpatVector polygons delineating the extent of the hurricane on ten different days. You want to know what the shortest paths are given that ships must go
around the polygon on that day. Calling shortest_paths
ten times with ten different SpatRasters would be very inefficient. This would assemble
the graph ten times and recompute also paths unaffected by the hurricanes, such as the path between India and Australia, in each iteration. Instead,
pass the SpatVector polygons to update_rst
. shortest_paths
then produces the shortest paths for a hurricane-free route and all ten
hurricane days, only reestimating the paths that are affected by the hurricane polygon on a specific day.
Applications to Earth should always pass a SpatRaster or RasterLayer to rst
. The option to use a matrix or a list of matrices is meant for
applications to other planets, non-geo-spatial settings, and users who cannot install the terra package on their system.
The largest source of runtime inefficiency is the quantity of non-NA pixels in the rst
grid. Limit the rst
argument to the relevant area.
E.g. crop the grid to the North Atlantic when computing shipping routes between Canada and France. And set regions through which the shortest path does
certainly not pass to NA.
shortest_paths
is optimized for computational performance. Most of its functions are written in C++ and it does not use a general purpose graph
library, but comes with its custom graph-theoritical implementation tailored to gridded inputs.
The vignette provides further details on the package.
If 'output = "distances"', the output is by default returned as a data table. If you want the result to be a data frame only, not a data table, set 'output_class' to '"data.frame"'. If 'output' is '"lines"' or '"both"', the the function returns a SpatVector, if 'rst' is a SpatRaster or a RasterLayer, and a list, if 'rst' is a matrix or a list of matrices. Explicitly setting 'output_class' to '"list"' returns a list in any case. 'output_class = "SpatVector"', however, returns a SpatVector only if 'rst' is a SpatRaster or a RasterLayer.
If output = "distances"
or output = "both"
, the distances
variable marks which points are connected. Unconnected point pairs have
an Inf
distance, if distance_type = "double"
or distance_type = "float"
, and an NA distance, if distance_type = "int"
or
distance_type = "unsigned short int"
. If output = "lines"
, the connected
variable marks which points are connected. Points are
connected, when it is possible to travel between them via non-NA cells in rst
.
Dijkstra, E. W. 1959. "A note on two problems in connexion with graphs." Numerische Mathematik 1 (1): 269–71.
# Generate example data set.seed(2L) input_grid <- terra::rast(crs = "epsg:4326", resolution = 2, vals = sample(c(1L, NA_integer_), 16200L, TRUE, c(0.8, 0.2))) origin_pts <- rnd_locations(5L, output_type = "SpatVector") origin_pts$name <- sample(letters, 5) destination_pts <- rnd_locations(5L, output_type = "SpatVector") # Compute distances shortest_paths(input_grid, origin_pts) shortest_paths(input_grid, origin_pts, bidirectional = TRUE) shortest_paths(input_grid, origin_pts, destination_pts) shortest_paths(input_grid, origin_pts, origin_names = "name") shortest_paths(input_grid, origin_pts, destination_pts, pairwise = TRUE) # Compute lines shortest_paths(input_grid, origin_pts, output = "lines") # Compute distances and lines shortest_paths(input_grid, origin_pts, output = "both") # Use custom transition function input_grid[input_grid == 1L] <- stats::runif(terra::freq(input_grid, value = 1L)$count, max = 1000) shortest_paths(input_grid, origin_pts, tr_fun = function(d, v1, v2) sqrt(d^2 + abs(v2 - v1)^2), tr_directed = FALSE) # Compute distances with grid updating barrier <- terra::vect("POLYGON ((-179 1, 30 1, 30 0, -179 0, -179 1))", crs = "epsg:4326") shortest_paths(input_grid, origin_pts, update_rst = barrier) barriers <- list(barrier, terra::vect("POLYGON ((0 0, 0 89, 1 89, 1 0, 0 0))", crs = "epsg:4326")) shortest_paths(input_grid, origin_pts, update_rst = barriers)
# Generate example data set.seed(2L) input_grid <- terra::rast(crs = "epsg:4326", resolution = 2, vals = sample(c(1L, NA_integer_), 16200L, TRUE, c(0.8, 0.2))) origin_pts <- rnd_locations(5L, output_type = "SpatVector") origin_pts$name <- sample(letters, 5) destination_pts <- rnd_locations(5L, output_type = "SpatVector") # Compute distances shortest_paths(input_grid, origin_pts) shortest_paths(input_grid, origin_pts, bidirectional = TRUE) shortest_paths(input_grid, origin_pts, destination_pts) shortest_paths(input_grid, origin_pts, origin_names = "name") shortest_paths(input_grid, origin_pts, destination_pts, pairwise = TRUE) # Compute lines shortest_paths(input_grid, origin_pts, output = "lines") # Compute distances and lines shortest_paths(input_grid, origin_pts, output = "both") # Use custom transition function input_grid[input_grid == 1L] <- stats::runif(terra::freq(input_grid, value = 1L)$count, max = 1000) shortest_paths(input_grid, origin_pts, tr_fun = function(d, v1, v2) sqrt(d^2 + abs(v2 - v1)^2), tr_directed = FALSE) # Compute distances with grid updating barrier <- terra::vect("POLYGON ((-179 1, 30 1, 30 0, -179 0, -179 1))", crs = "epsg:4326") shortest_paths(input_grid, origin_pts, update_rst = barrier) barriers <- list(barrier, terra::vect("POLYGON ((0 0, 0 89, 1 89, 1 0, 0 0))", crs = "epsg:4326")) shortest_paths(input_grid, origin_pts, update_rst = barriers)