An introduction to geoGraph
geoGraph: walking through the geographic space using graphs.
This document describes the geoGraph package for the R software. geoGraph aims at implementing graph approaches for geographic data. In geoGraph, a given geographic area is modelled by a fine regular grid, where each vertice has a set of spatial coordinates and a set of attributes, which can be for instance habitat descriptors, or the presence/abundance of a given species. ‘Travelling’ within the geographic area can then be easily modelled as moving between connected vertices. The cost of moving from one vertex to another can be defined according to attribute values, which allows for instance to define friction routes based on habitat.
geoGraph harnesses the full power of graph algorithms implemented in R by the graph and RBGL (R Boost Graph Library) packages. In particular, RBGL is an interface between R and the comprehensive Boost Graph Library in C++, which provides fast and efficient implementations of a wide range of graph algorithms. Once we have defined frictions for an entire geographic area, we can easily, for instance, find the least costs path from one location to another, or find the most parsimonious way of connecting a set of locations.
Interfacing spatial data and graphs can be a complicated task. The
purpose of geoGraph is to provide tools to achieve and simplify
this ‘preliminary’ step. This is achieved by defining new classes of
objects which are essentially geo-referenced graphs with node attributes
(gGraph objects), and interfaced spatial data
(gData objects). In this vignette, we show how to install
geoGraph, construct and handle
gGraph/gData objects, and illustrate some
basic features of graph algorithms.
First steps
Installing the package
All the following instructions should be entered from a new R session to avoid errors due to installing attached packages.
devtools is also needed to install geoGraph:
Then, to install geoGraph, simply type:
Once installed, the package can be loaded using:
If you have an error regarding missing packages, you may need to install manually the packages graph and RBGL from Bioconductor:
And then attempt to reinstall geoGraph from GitHub.
Data representation
Data representation refers to the way a given type of data is handled
by a computer program. Two types of objects are used in
geoGraph: gGraph, and gData objects.
Both objects are defined as formal (S4) classes and often have methods
for similar generic function (e.g. getNodes is defined for
both objects). Essentially, gGraph objects contain
underlying layers of informations, including a spatial grid and possibly
node attributes, and covering the area of interest. gData
are sets of locations (like sampled sites, for instance) which have been
interfaced to a gGraph object, to allow further
manipulations such as finding paths on the grid between pairs of
locations.
gGraph objects
The definition of the formal class gGraph can be
obtained using:
## Class "gGraph" [package "geoGraph"]
##
## Slots:
##
## Name: coords nodes.attr meta graph
## Class: matrix data.frame list graphNELand a new empty object can be obtained using the constructor:
##
## === gGraph object ===
##
## @coords: spatial coordinates of 0 nodes
## lon lat
##
## @nodes.attr: 0 nodes attributes
## data frame with 0 columns and 0 rows
##
## @meta: list of meta information with 0 items
##
## @graph:
## A graphNEL graph with undirected edges
## Number of Nodes = 0
## Number of Edges = 0The documentation ?gGraph explains the basics about the
object’s content. In a nutshell, these objects are spatial grids with
nodes and segments connecting neighbouring nodes, and additional
informations on the nodes or on the graph itself. coords is
a matrix of longitudes and latitudes of the nodes.
nodes.attr is a data.frame storing attributes of the nodes,
such as habitat descriptors; each row corresponds to a node of the grid,
while each column corresponds to a variable. meta is a list
containing miscellanous informations about the graph itself. There is no
contraint applying to the components of the list, but some typical
components such as $costs or $colors will be
recognised by certain functions. For instance, you can specify plotting
rules for representing a given node attribute by a given color by
defining a component $colors. Similarly, you can associate
costs to a given node attribute by defining a component
$costs. An example of this can be found in already existing
gGraph objects. For instance, worldgraph.10k
is a graph of the world with approximately 10,000 nodes, and only
on-land connectivity
(\textit{i.e.no travelling on the seas). ```{r } worldgraph.10k worldgraph.10k@meta ``` Lastly, thegraphcomponent is agraphNEL`
object, which is the standard class for graphs in the graph and
RBGL packages. This object contains all information on the
connections between nodes, and the weights (costs) of these
connections.
Four main gGraph are provided with geoGraph:
rawgraph.10k, rawgraph.40k,
worldgraph.10k, and worldgraph.40k. These
datasets are available using the command data. The grid
used in these datasets are the best geometric approximation of a regular
grid for the surface of a sphere. One advantage of working with these
grids is that we do not have to use a projection for geographic
coordinates, which is a usual issue in regular GIS.
The difference between rawgraphs and worldgraphs is that the first
are entirely connected, while in the second connections occur only on
land. Numbers 10k' and40k’ indicate that the grids consist
of roughly 10,000 and 40,000 nodes. For illustrative purposes, we will
often use the 10k grids, since they are less heavy to handle. For most
large-scale applications, the 40k versions should provide sufficient
resolution. New gGraph can be constructed using the
constructor (new(...)), but this topic is not documented in
this vignette.
gData objects
gData are essentially sets of locations that are
interfaced with a gGraph object. During this operation,
each location is assigned to the closest node on the grid of the
gGraph, then allowing for travelling between locations
using the grid. Then, it is for instance possible to find the shortest
path between two locations through various types of habitats.
Like for gGraph, the content of the formal class
gData can be obtained using:
## Class "gData" [package "geoGraph"]
##
## Slots:
##
## Name: coords nodes.id data gGraph.name
## Class: matrix character ANY characterand a new empty object can be obtained using the constructor:
##
## === gData object ===
##
## @coords: spatial coordinates of 0 nodes
## lon lat
##
## @nodes.id: nodes identifiers
## character(0)
##
## @data: data
## NULL
##
## Associated gGraph:As before, the description of the content of these objects can be
found in the documentation (?gData). coords is
a matrix of xy (longitude/latitude) coordinates in which each row is a
location. nodes.id is vector of characters giving the name
of the vertices matching the locations; this is defined automatically
when creating a new gData, or using the function
closestNode. data is a slot storing data
associated to the locations; it can be any type of object, but a
data.frame should cover most requirements for storing data. Note that
this object should be subsettable (i.e. the [ operator
should be defined), so that data can be subsetted when subsetting the
gData object. Lastly, the slot gGraph.name
contains the name of the gGraph object to which the
gData has been interfaced.
Contrary to gGraph objects, gData objects
will frequently be constructed by the user. In the next sections, we
illustrate how we can build and use gData objects from a
set of locations.
Using geoGraph
An overview of the material implemented in the package is summarized the package’s manpage, accessible via:
The html version of this manpage may be preferred to browse more easily the content of geoGraph; it is accessible by typing:
To revert help back to text mode, simply type:
In the following, we go through various tasks that can be achieve using geoGraph.
Importing geographic data
GeoGraphic data consist of a set of locations, possibly accompanied
by additional information. For instance, one may want to study the
migrations amongst a set of biological populations with known geographic
coordinates. In geoGraph, geographic data are stored in
gData objects. These objects match locations to the closest
nodes on a grid (a gGraph object), and store additional
data if needed.
As a toy example, let us consider four locations: Bordeaux (France), London (UK), Malaga (Spain), and Zagreb (Croatia). Since we will be working with a crude grid (10,000 nodes), locations need not be exact. We enter the longitudes and latitudes (in this order, that is, xy coordinates) of these cities in decimal degrees, as well as approximate population sizes:
Bordeaux <- c(-1, 45)
London <- c(0, 51)
Malaga <- c(-4, 37)
Zagreb <- c(16, 46)
cities.dat <- rbind.data.frame(Bordeaux, London, Malaga, Zagreb)
colnames(cities.dat) <- c("lon", "lat")
cities.dat$pop <- c(1e6, 13e6, 5e5, 1.2e6)
row.names(cities.dat) <- c("Bordeaux", "London", "Malaga", "Zagreb")
cities.dat## lon lat pop
## Bordeaux -1 45 1.0e+06
## London 0 51 1.3e+07
## Malaga -4 37 5.0e+05
## Zagreb 16 46 1.2e+06We load a gGraph object which contains the grid that
will support the data:
##
## === gGraph object ===
##
## @coords: spatial coordinates of 10242 nodes
## lon lat
## 1 -180.00 90.00
## 2 144.00 -90.00
## 3 -33.78 27.19
## ...
##
## @nodes.attr: 1 nodes attributes
## habitat
## 1 sea
## 2 sea
## 3 sea
## ...
##
## @meta: list of meta information with 2 items
## [1] "$colors" "$costs"
##
## @graph:
## A graphNEL graph with undirected edges
## Number of Nodes = 10242
## Number of Edges = 6954
(we could use worldgraph.40k for a better resolution).
In this figure, each node is represented with a color depending on the
habitat type, either ‘sea’ (blue) or ‘land’ (green). We are going to
interface the cities data with this grid; to do so, we create a
gData object using new (see
?gData object):
cities <- new("gData", coords = cities.dat[, 1:2], data = cities.dat[, 3, drop = FALSE], gGraph.name = "worldgraph.10k")
cities##
## === gData object ===
##
## @coords: spatial coordinates of 4 nodes
## lon lat
## 1 -1 45
## 2 0 51
## 3 -4 37
## ...
##
## @nodes.id: nodes identifiers
## 1 2 3
## "5774" "6413" "4815"
## ...
##
## @data: 4 data
## pop
## Bordeaux 1.0e+06
## London 1.3e+07
## Malaga 5.0e+05
## ...
##
## Associated gGraph: worldgraph.10k
This figure illustrates the matching of original locations (black
crosses) to nodes of the grid (red circles). As we can see, an issue
occured for Bordeaux, which has been assigned to a node in the sea (in
blue). Locations can be re-assigned to nodes with restrictions for some
node attribute values using closestNode; for instance, here
we constrain matching nodes to have an habitat value
(defined as node attribute in worldgraph.10k) equalling
land (green points):
cities <- closestNode(cities, attr.name = "habitat", attr.value = "land")
plot(cities, type = "both", reset = TRUE)
plotEdges(worldgraph.10k)
Now, all cities have been assigned to a
land' node of the grid (again, better accuracy will be gained on 40k or finer grids - we use 10k for illustrative purposes only). Content ofcitiescan be accessed via various accessors (see?gData`).
For instance, we can retrieve original locations, assigned nodes, and
stored data using:
## lon lat
## 5775 -1 45
## 6413 0 51
## 4815 -4 37
## 7699 16 46## 5774 6413 4815 7699
## "5775" "6413" "4815" "7699"## pop
## Bordeaux 1.0e+06
## London 1.3e+07
## Malaga 5.0e+05
## Zagreb 1.2e+06We can also get the coordinates of the matching nodes (, red circle on previous figure) using:
## lon lat
## 5775 1.002e-05 43.73
## 6413 1.002e-05 51.38
## 4815 -3.788e+00 37.75
## 7699 1.548e+01 46.74More interestingly, we can now retrieve all the geographic
information contained in the underlying grid (, gGraph
object) as node attributes:
## habitat
## 5775 land
## 6413 land
## 4815 land
## 7699 landIn this example, the information stored in
worldgraph.10k is rather crude: habitat only
distinguishes the land from the sea. However, more complex habitat
information could be incorporated, for instance from GIS shapefiles (see
dedicated section below).
Visualizing data
An essential aspect of spatial analysis lies in visualizing the data.
In geoGraph, the spatial grids (gGraph) and
spatial data (gData) can be plotted and browsed using a
variety of functions.
Plotting gGraph objects
Displaying a gGraph object is done through
plot and points functions. The first opens a
new plotting region, while the second draws in the current plotting
region; functions have otherwise similar arguments (see
?plot.gGraph).
By default, plotting a gGraph displays the grid of nodes
overlaying a shapefile (by default, the landmasses). Edges can be
plotted at the same time (argument edges), or added
afterwards using plotEdges. If the gGraph
object possesses an adequately formed meta$colors
component, the colors of the nodes are chosen according to the node
attributes and the color scheme specified in meta$colors.
Alternatively, the color of the nodes can be specified via the
col argument in plot/points.
Here is an example using worldgraph.10k:
## habitat color
## 1 sea blue
## 2 land green
## 3 mountain brown
## 4 landbridge light green
## 5 oceanic crossing light blue
## 6 deselected land lightgray## habitat
## 1 sea
## 2 sea
## 3 sea
## 4 sea
## 5 sea
## 6 sea## habitat
## deselected land land sea
## 290 2632 7320
It may be worth noting that plotting gGraph objects
involves plotting a fairly large number of points and edges. On some
graphical devices, the resulting plotting can be slow. For instance, one
may want to disable cairo under linux: this graphical
device yields better graphics than Xlib, but at the expense
of increase computational time. To switch to Xlib,
type:
and to revert to cairo, type:
Zooming in and out, sliding, etc.
In practice, it is often useful to be able to peer at specific
regions, and more generally to navigate inside the graphical
representation of the data. For this, we can use the interactive
functions geo.zoomin, geo.zoomout,
geo.slide, geo.back,
geo.bookmark, and geo.goto. The zoom and slide
functions require to left-click on the graphics to zoom in, zoom out, or
slide to adjacent areas; in all cases, a right click ends the function.
Also note that geo.zoomin can accept an argument specifying
a rectangular region, which will be adapted by the function to fit best
a square area with similar position and centre, and zoom to this area
(see ?geo.zoomin). geo.bookmark and
geo.goto respectively set and go to a bookmark,
i.e. a tagged area. This is most useful when one has to switch
between distant areas repeatedly.
Here are some examples based on the previous plotting of
worldgraph.10k: Zooming in:
Zooming out:
Sliding to the east:
One important thing which makes plotting gGraph objects
different from most other plotting in R is that geoGraph
keeps the changes made to the plotting area in memory. This allows to
undo one or several moves using geo.back. Moreover, even if
the graphical device is killed, plotting again a gGraph
will use the old parameters by default. To disable this behavior, set
the argument reset=TRUE when calling upon
plot. Technically, this ‘plotting memory’ is implemented by
storing plotting information in an environment defined as the hidden
environment geoGraph:::.geoGraphEnv:
## [1] "bookmarks" "last.plot" "last.plot.param" "last.points"
## [5] "psize" "sticky.points" "usr" "zoom.log"You can inspect individual variables within this environment:
## $psize
## [1] 0.5
##
## $pch
## [1] 19However, it is recommended not to modify these objects directly,
unless you really know what you are doing. In any case, plotting a
gGraph object with argument reset=TRUE will
remove previous plotting history and undo possible wrong
manipulations.
Plotting gData objects
gData objects are by default plotted overlaying the
corresponding gGraph. For instance, using the
cities example from above:
Note the argument reset=TRUE, which tells the plotting
function to adapt the plotting area to the geographic extent of the
dataset.
To plot additional information, it can be useful to extract the
spatial coordinates from the data. This is achieved by
getCoords. This method takes an extra argument
original, which is TRUE if original spatial coordinates are
seeked, or FALSE for coordinates of the nodes on the grid. We can use
this to represent, for instance, the population sizes for the different
cities:
transp <- function(col, alpha = .5) {
res <- apply(col2rgb(col), 2, function(c) rgb(c[1] / 255, c[2] / 255, c[3] / 255, alpha))
return(res)
}
plot(cities, reset = TRUE)
par(xpd = TRUE)
text(getCoords(cities) + -.5, rownames(getData(cities)))
symbols(getCoords(cities)[, 1], getCoords(cities)[, 2], circ = sqrt(unlist(getData(cities))), inch = .2, bg = transp("red"), add = TRUE)
Editing gGraphs
Editing graphs is an essential task in geoGraph. While
available gGraph objects provide a basis to work with (see
?worldgraph.10k), one may want to adapt a graph to a
specific case. For instance, connectivity should be defined according to
biological knowledge of the organism under study. gGraph
can be modified in different ways: by changing the connectivity, the
costs of edges, or the attribute values.
Changing the global connectivity of a gGraph
There are two main ways of changing the connectivity of a
gGraph, which match two different objectives. The first
approach is to perform global and systematic changes of the connectivity
of the graph. Typically, one will want to remove all connections over a
given type of landscape, which is uncrossable by the organism under
study. Let’s assume we are interested in saltwater fishes. To model fish
dispersal, we have to define a graph which connects only nodes
overlaying the sea. We load the gGraph object
rawgraph.10k, and zoom in to a smaller area (Madagascar) to
illustrate changes in connectivity:
We shall set a bookmark for this area, in case we would want to get back to this place later on:
##
## Bookmark ' madagascar 'saved.What we now want to do is remove all but sea-sea connections. To do so, the easiest approach is to i) define costs for the edges based on habitat, with land being given large costs and ii) remove all edges with large costs.
Costs of a given node attribute (here,
habitat') are indicated in themeta$costs` slot:
## habitat cost
## 1 sea 100
## 2 land 1
## 3 mountain 10
## 4 landbridge 5
## 5 oceanic crossing 20
## 6 deselected land 100newGraph <- rawgraph.10k
newGraph@meta$costs[2:6, 2] <- 100
newGraph@meta$costs[1, 2] <- 1
newGraph@meta$costs## habitat cost
## 1 sea 1
## 2 land 100
## 3 mountain 100
## 4 landbridge 100
## 5 oceanic crossing 100
## 6 deselected land 100We have just changed the costs associated to habitat type, but this
change is not yet effective on edges between nodes. We use
setCosts to set the cost of an edge to the average of the
costs of its nodes:
On this new graph, we represent the edges with a width inversely
proportional to the associated cost; that is, bold lines for easy
travelling and light edges/dotted lines for more costly mouvement. This
is not enough yet, since travelling on land is still possible. However,
we can tell geoGraph to remove all edges associated to too
strong a cost, as defined by a given threshold (using
dropDeadEdges). Here, only sea-sea connections shall be
retained, that is, edges with cost 1.
Here we are: newGraph only contains connections in the
sea. Note that, although we restrained the plotting area to Madagascar,
this change is effective everywhere. For instance, travelling to the
nort-west Australian coasts:
##
## Bookmark ' australia 'saved.Changing local properties of a gGraph
A second approach to changing a gGraph is to refine the
graph by hand, adding or removing locally some connections, or altering
the attributes of some nodes. This can be necessary to connect
components such as islands to the main landmasses, or to correct
erroneous data. Adding and removing edges from the grid of a
gGraph can be achieved by geo.add.edges and
geo.remove.edges, respectively. These functions are
interactive, and require the user to select individual nodes or a
rectangular area in which edges are added or removed. See
?geo.add.edges for more information on these functions. For
instance, we can remove a few odd connections in the previous graph,
near the Australian coasts (note that we have to save the changes using
<-):
When adding connections within an area or in an entire graph, node
addition is based on another gGraph,i.e. only
connections existing in another gGraph serving as reference
can be added to the current gGraph. For graphs based on 10k
or 40k grids, the raw graphs provided in geoGraph should be
used, (rawgraph.10k, rawgraph.40k), since they
are fully connected.
In addition to changing grid connectivity, we may also want to modify
the attributes of specific nodes. This is again done interactively,
using the function geo.change.attr. For instance, here, we
define a new value shalowwater (plotted in light blue) for
the attribute habitat, selecting affected nodes using the
‘area’ mode first, and refining the changes using the ‘point’ mode:
plot(newGraph, edge = TRUE)
temp <- geo.change.attr(newGraph, mode = "area", attr.name = "habitat", attr.value = "shallowwater", newCol = "deepskyblue")
temp <- geo.change.attr(temp, attr.name = "habitat", attr.value = "shallowwater", newCol = "deepskyblue")
newGraph <- temp## habitat color
## 1 sea blue
## 2 land green
## 3 mountain brown
## 4 landbridge light green
## 5 oceanic crossing light blue
## 6 deselected land lightgray
## 7 shallowwater deepskyblue
Again, note that the changes made to the graph have to be save in an
object (using <-) to be effective.
Extracting information from GIS shapefiles
An important feature of geoGraph is serving as an interface
between geographic information system (GIS) layers and
geographic data. As currently implemented, geoGraph can extract
information from shapefiles with the Arc GIS (http://www.esri.com/software/arcgis/index.html) format,
using the function extractFromLayer. Here, we illustrate
this procedure using the shapefile world-countries.shp
provided with the package. geoGraph still relies on
SpatialPolygonDataFrame objects from the sp
package. A move to sf objects is planned, but for the
moment we have to use sf::st_read() to read a GIS shapefile
into an sf object, and then cast it to a
SpatialPolygonDataFrame with
sf::as_Spatial():
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUEworld.countries <- sf::st_read(system.file("files/shapefiles/world-countries.shp", package = "geoGraph"))## Reading layer `world-countries' from data source
## `/tmp/Rtmpx2Syg1/Rinste40764dc7be/geoGraph/files/shapefiles/world-countries.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 239 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -89.9 xmax: 181.8 ymax: 84.93
## CRS: NA## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"## Length Class Mode
## 239 SpatialPolygonsDataFrame S4The summary of world.countries shows the data
(attributes) stored in the layer. Let us assume that we are interested
in retrieving continent and country information for the
worldgraph.10k object. Note that
extractFromLayer can extract information to other types of
objects than gGraph (see
?extractFromLayer)
## habitat
## deselected land: 290
## land :2632
## sea :7320newGraph <- extractFromLayer(worldgraph.10k, layer = world.countries, attr = c("CONTINENT", "NAME"))
summary(getNodesAttr(newGraph))## habitat CONTINENT NAME
## deselected land: 290 Length:10242 Length:10242
## land :2632 Class :character Class :character
## sea :7320 Mode :character Mode :characterThe new object newGraph is a gGraph which
now includes, for each node of the grid, the corresponding continent and
country retrieved from the GIS layer. We can use the newly acquired
information for plotting newGraph, by defining new color
rules:
temp <- unique(getNodesAttr(newGraph)$"NAME")
col <- c("transparent", rainbow(length(temp) - 1))
colMat <- data.frame(NAME = temp, color = col)
head(colMat)## NAME color
## 1 <NA> transparent
## 2 Antartica #FF0000
## 3 Saudi Arabia #FF0B00
## 4 Yemen #FF1500
## 5 Somalia #FF2000
## 6 China #FF2B00## NAME color
## 140 Latvia #FF0040
## 141 Belarus #FF0035
## 142 Eritrea #FF002A
## 143 Djibouti #FF0020
## 144 East Timor #FF0015
## 145 Jordan #FF000B
This information could in turn be used to define costs for travelling on the grid. For instance, one could import habitat descriptors from a GIS, use these values to formulate a habitat model, and derive costs for dispersal on the grid.
As soon as a GIS layer has been extracted to a gGraph,
this information becomes also available for any gData
interfaced with this object. For instance, we can re-use the
cities example defined in a previous section, and interface
it with newGraph to retrieve continent and country
information for the cities of the dataset:
## lon lat pop
## Bordeaux -1 45 1.0e+06
## London 0 51 1.3e+07
## Malaga -4 37 5.0e+05
## Zagreb 16 46 1.2e+06cities <- new("gData", coords = cities.dat[, 1:2], data = cities.dat[, 3, drop = FALSE], gGraph.name = "newGraph")
cities <- closestNode(cities, attr.name = "habitat", attr.value = "land")
getData(cities)## pop
## Bordeaux 1.0e+06
## London 1.3e+07
## Malaga 5.0e+05
## Zagreb 1.2e+06## habitat CONTINENT NAME
## 5775 land Europe France, Metropolitan
## 6413 land Europe United Kingdom
## 4815 land Europe Spain
## 7699 land Europe AustriaFinding least-cost paths
One of the most useful applications of geoGraph is the
research of least-cost paths between couples of locations. This can be
achieved using the functions dijkstraFrom and
dijkstraBetween on a gData object which
contains all the locations of interest. These functions return
least-cost paths with the format gPath.
dijkstraFrom compute the paths from a given node of the
grid to all locations of the gData, while
dijkstraBetween computes the paths between pairs of
locations of the gData. Below, we detail the example of the
documentation of these functions, which uses the famous dataset of
native Human populations, HGDP:
##
## === gData object ===
##
## @coords: spatial coordinates of 52 nodes
## lon lat
## 1 -3 59
## 2 39 44
## 3 40 61
## ...
##
## @nodes.id: nodes identifiers
## 28179 11012 22532
## "26898" "11652" "22532"
## ...
##
## @data: 52 data
## Population Region Label n Latitude Longitude Genetic.Div
## 1 Orcadian EUROPE 1 15 59 -3 0.7259
## 2 Adygei EUROPE 2 17 44 39 0.7298
## 3 Russian EUROPE 3 25 61 40 0.7320
## ...
##
## Associated gGraph: worldgraph.40k
Populations of the dataset are shown by red circles, while the
underlying grid (worldgraph.40k) is represented with colors
depending on habitat (blue: sea; green: land; pink: coasts). Population
genetics predicts that genetic diversity within populations should decay
as populations are located further away from the geographic origin of
the species. Here, we verify this relationship for a theoretical origin
in Addis abeba, Ethiopia. We shall seek all paths through landmasses to
the HGDP populations.
First, we check that all populations are connected on the grid using
isConnected:
## [1] TRUENote that in practice, we may often want to assess graphically the
connectivity of the underlying grid, especially if not all locations of
the gData are connected. This can be done using
connectivityPlot, which has methods for both
gGraph and gData, and represents different
connected components using different colors. For instance, for
worldgraph.10k:
Since all locations in hgdp are connected, we can
proceed further. We have to set the costs of edges in the
gGraph grid. To do so, we can choose between i) strictly
uniform costs (using dropCosts) ii) distance-based costs –
roughly uniform – (using setDistCosts) or iii)
attribute-driven costs (using setCosts).
We shall first illustrate the strictly uniform costs. After setting a
gGraph with uniform costs, we use dijkstraFrom
to find the shortest paths between Addis abeba and the populations of
hgdp:
myGraph <- dropCosts(worldgraph.40k)
hgdp@gGraph.name <- "myGraph"
addis <- cbind(38, 9)
ori <- closestNode(myGraph, addis)
paths <- dijkstraFrom(hgdp, ori)## Loading required package: RBGLThe object paths contains the identified paths, which
are stored as a list with class gPath (see
?gPath). Paths can be plotted easily:
addis <- as.vector(addis)
plot(myGraph, col = NA, reset = TRUE)
plot(paths)
points(addis[1], addis[2], pch = "x", cex = 2)
text(addis[1] + 35, addis[2], "Addis abeba", cex = .8, font = 2)
points(hgdp, col.node = "black")
In this graph, each path is plotted with a different color, but
several paths overlap in several places. We can extract the distances
from the
origin' usinggPath2dist, and then examine the relationship between genetic diversity within populations (stored inhgdp`)
and the distance from the origin:
div <- getData(hgdp)$"Genetic.Div"
dgeo.unif <- gPath2dist(paths, res.type = "vector")
plot(div ~ dgeo.unif, xlab = "GeoGraphic distance (arbitrary units)", ylab = "Genetic diversity")
lm.unif <- lm(div ~ dgeo.unif)
abline(lm.unif, col = "red")
summary(lm.unif)##
## Call:
## lm(formula = div ~ dgeo.unif)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.07327 -0.00660 0.00074 0.01015 0.05449
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.70e-01 4.58e-03 168.2 <2e-16 ***
## dgeo.unif -8.39e-04 5.31e-05 -15.8 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0185 on 50 degrees of freedom
## Multiple R-squared: 0.833, Adjusted R-squared: 0.83
## F-statistic: 250 on 1 and 50 DF, p-value: <2e-16
Alternatively, we can use costs based on habitat. As a toy example, we will consider that coasts are four times more favourable for dispersal than the rest of the landmasses. We define these new costs, and then compute and plot the corresponding shortest paths:
## habitat cost
## 1 sea 100
## 2 land 1
## 3 mountain 10
## 4 landbridge 5
## 5 oceanic crossing 20
## 6 deselected land 100
## 7 coast 0.25plot(newGraph, col = NA, reset = TRUE)
plot(paths.2)
points(addis[1], addis[2], pch = "x", cex = 2)
text(addis[1] + 35, addis[2], "Addis abeba", cex = .8, font = 2)
points(hgdp, col.node = "black")
The new paths are slightly different from the previous ones. We can examine the new relationship with genetic distance:
dgeo.hab <- gPath2dist(paths.2, res.type = "vector")
plot(div ~ dgeo.hab, xlab = "GeoGraphic distance (arbitrary units)", ylab = "Genetic diversity")
lm.hab <- lm(div ~ dgeo.hab)
abline(lm.hab, col = "red")
summary(lm.hab)##
## Call:
## lm(formula = div ~ dgeo.hab)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.11183 -0.00976 0.00133 0.01216 0.06413
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.770137 0.007174 107.36 < 2e-16 ***
## dgeo.hab -0.001421 0.000145 -9.79 3.2e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0265 on 50 degrees of freedom
## Multiple R-squared: 0.657, Adjusted R-squared: 0.651
## F-statistic: 95.9 on 1 and 50 DF, p-value: 3.21e-13
Of course, the distinction between coasts and inner landmasses is a somewhat poor description of habitat. In practice, complex habitat models can be used as simply.