graph4lg
Build Graphs for Landscape Genetics Analysis

Package index
Search the graph4lg package
Vignettes
Functions
153
Source code
72
Man pages
80
Home
/
CRAN
/
graph4lg
/
3 - Landscape graph construction and analysis with Graphab and `graph4lg`

3 - Landscape graph construction and analysis with Graphab and `graph4lg`
In graph4lg: Build Graphs for Landscape Genetics Analysis 

knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)

library(graph4lg)
library(igraph)
library(ggplot2)

Introduction

The rationale of graph4lg package in R is to make the construction and analysis of genetic and landscape graphs easier for landscape genetic studies (hence the name graph4lg, meaning Graphs for Landscape Genetics). This package provides users with tools for:

Each one of the included vignettes focuses on one of these points. This third vignette will focus on landscape graph construction and analysis. It will describe the package functions allowing users to:

NOTA BENE: The package graph4lg integrates functions making possible the construction and analysis of landscape graphs but most of them are only wrappers that launch computations with Graphab software.

What is Graphab?

Graphab is an open source software tool dedicated to the modelling of landscape networks [@foltete2012software; @foltete2021graphab]. It was developed in the ThéMA laboratory (Besançon, France) to make possible the creation of landscape graphs and to integrate a complete set of connectivity analysis functions in a single application.

Integrating Graphab functionalities into graph4lg aims at facilitating the import of Graphab output into R for subsequent analyses, including the comparison of genetic graphs and landscape graphs. All the functions calling Graphab software are described in graph4lg help files. However, users are invited to read the complete documentation manuals as well as other resources on Graphab project website. Some computations are not possible when using graph4lg but are possible when using directly the Graphab software which can be downloaded for free on the same website. @foltete2012software and @foltete2021graphab papers also provide users with a synthetic description of the tool. In 2024, new functionalities were introduced in Graphab 3.0 to make it possible to include several types of habitats into the same habitat graph, creating so-called 'multiple habitat graphs'. This modelling option is described by @savary2024multiple and can now be used in graph4lg (version >= 2.0).

The function get_graphab from graph4lg package automatically downloads Graphab 3.0 in the user's machine, provided Java is installed on the machine. If it has already been downloaded, a message is displayed. Once Graphab has been downloaded in a machine, it is not necessary to download it again. Note that users familiar with the functions calling graphab-2.8.jar can still use the previous versions of graph4lg by installing them from the gitlab repository of the package. 

get_graphab()

In the table below, we present all the functions from graph4lg calling Graphab software:

| Function | Description | | ---- | -------------------------------------------------- | | get_graphab | Check if Graphab software (graphab-3.0.jar) is in the user's machine and downloads it if necessary, provided Java is installed. | | graphab_project | Creates a Graphab project from a raster file (.tif) | | graphab_habitat | Creates an habitat in the Graphab project from either the raster file or a vector spatial layer | | graphab_capacity | Modifies the capacity of habitat patches based on additional data or patch neighbourhood raster calculations | | graphab_link | Creates a set of links (Euclidean or least cost paths, planar or complete) | | graphab_graph | Creates a graph from a set of links | | graphab_merge_graph | Merges existing graphs from a set of links | | graphab_metric | Computes graph-theoretic metrics from a graph | | graphab_modul | Partition the graph into modules | | graphab_corridor | Creates corridor surrounding least-cost paths | | graphab_interpol | Interpolates Graphab metrics over continuous raster fields | | graphab_pointset | Imports a set of points in the project and relates them to habitat patches and corresponding metrics | | graphab_metapatch | Creates a new habitat type made of meta-patches, each corresponding to the connected components of a pruned/incomplete graph | | graphab_show | Lists the elements already included in the Graphab project | | graphab_project_desc | Describes all the project parameters | | get_graphab_metric | Imports a table with the metrics computed and the patch properties | | get_graphab_linkset | Imports a table with the link properties for a given link set | | get_graphab_linkset_cost | Gets the parameters of a cost-distance based link set | | get_graphab_raster_codes | Lists the categorical values of the raster used to create the project | | graphab_to_igraph | Imports a graph created with Graphab as an igraph graph object with node and link attributes |

NOTA BENE: This vignette only provides a basic introduction to all the functionalities offered by the Graphab software application. All the cases cannot be covered here and users interested in more options can find more information in the Graphab website and in publications of the Graphab group (in particular in @foltete2021graphab and @savary2024multiple for multiple habitat graphs).

Data used in this vignette

Here, we do not rely on genetic data sets and will only use Graphab projects already created for the vignettes.

load(file = paste0(system.file('extdata', package = 'graph4lg'), 
                   "/", "res_g.RDa"))

Landscape graph construction

Project creation

When using Graphab, the main input data consists of an 8-bit raster file (.tif) with discrete cell values. One (or several) of these values will correspond to habitat cells on the raster. Contiguous habitat cells will form habitat patches (8-neighbourhood criterion), as soon as they exceed the minimum patch area when specified.

The first step of habitat network modelling with Graphab is the creation of a Graphab project. This project is given a name (proj_name), which will be the name of a directory containing geographical information layers relative to the project as well as a .xml file ("proj_name.xml"), which is the main project file.

The function graphab_project creates the project directory called proj_name either in the current working directory (default) or in the directory whose path is given as a proj_path argument. The function takes as arguments:

The three last arguments are included in almost every function calling Graphab, and won't be mentioned anymore in the next sections.

As an example, we will create a landscape graph from a simulated landscape raster with 5 land use categories (see map below).

rast <- terra::rast(paste0(system.file('extdata', package = 'graph4lg'), 
                   "/", "rast_simul50.tif"))
# convert the raster into a df while keeping coords
r_df <- as.data.frame(rast,
                      xy = TRUE)
# add a categorical field with the raster codes
r_df$code <- as.factor(r_df$rast_simul50)

# Plot it with ggplot2
g <- ggplot(r_df, aes(x = x, y = y)) +
  geom_tile(aes(fill = code)) + coord_equal()+
  theme_bw()+
  #scale_fill_brewer(palette="Dark2")+
  scale_fill_manual(values = c("#396D35", "#FB9013", "#EDC951", "#80C342", "black", "#396D35"),
                    labels = c("0 - Forest", "1 - Shrublands", "3 - Crops",
                               "4 - Grasslands","5 - Artificial areas", "6 - Forest"),
                    name = "Land use type")+
  labs(x="Longitude",y="Latitude")
g

We can create a project to assess habitat connectivity in this landscape using the following command:

proj_name <- "graphab_example"

graphab_project(proj_name = proj_name,
                raster = "rast_simul50.tif")

proj_name <- "graphab_example"

Habitat creation

Once the project has been created, we can use the function graphab_habitat to create a new habitat type in the Graphab project (proj_name). As a result, a new repository is created in the project repository, including a geopackage layer with the corresponding habitat patches (polygons). The function takes as arguments:

Here, we create a first habitat type with habitat patches corresponding to contiguous patches of raster cells equal to 0 or 5, as soon as they are larger than 200 hectares. We assume they correspond to forests.

# Habitat creation
graphab_habitat(proj_name = proj_name, 
                name = "forest",
                type = "raster", 
                rast_codes = c(0, 5), 
                minarea = 200)

Since Graphab 3.0, it is possible to include several types of habitats in the same Graphab project. We illustrate that below by creating an habitat type corresponding to shrubland patches. 

# Bushland
graphab_habitat(proj_name = proj_name,
                name = "shrubland",
                type = "raster",
                rast_codes = 1,
                minarea = 20)

We can check that the two habitat types have been successfully created by listing the elements of the project, using the graphab_show function. Its only argument is the path to the .xml file storing all the project information:

graphab_show(proj_path = "graphab_example/graphab_example.xml")

graphab_show(proj_path = paste0(system.file('extdata', 
                                            package = 'graph4lg'), 
                   "/", "graphab_example/graphab_example.xml"))

Since Graphab 3.0, it is also possible to create habitat types by providing the corresponding patches as point or polygon entities of a geopackage or shapefile vectorial layer. In this case, a column of the attribute table must include the values used to set the capacity of the habitat patches. The name of this column is provided with the argument vec_capa_field. For example:

graphab_habitat(proj_name = proj_name,
                name = "grassland",
                type = "vector",
                vec_layer = "grassland_patches.gpkg",
                vec_capa_field = "capa")

Link set creation

The project and habitat types have been created and we will now create a link set between the habitat patches. Although the term "link set" is used at this stage, these links are not graph links but just spatial paths between habitat patches. Both the patches and these paths will eventually serve as the basis for the construction of habitat graphs having nodes connected by graph links.

The link set can be built as the set of straight lines ("as the crow flies") connecting the patches, weighted by their corresponding geographical Euclidean distances. Yet, we can also compute least cost paths between habitat patches on the raster, which are more frequently used in connectivity analyses tailored to specific taxa whose dispersal behaviour is partially known. For that purpose, we need to create a data.frame specifying the cost values associated with every raster cell value, usually based on prior knowledge or analyses about dispersal of the study species in heterogeneous landscapes. In the following example, these values are:

| Code | Description | Cost value | | ---- | ------------------------ | ----------- | | 0 | Forests (habitat) | 1 | | 1 | Shrublands | 5 | | 2 | Crops | 60 | | 3 | Grasslands | 40 | | 4 | Artificial areas | 1000 | | 5 | Forest (habitat) | 1 |

They are stored into the cost object (data.frame):

cost <- data.frame(code = 0:5,
                   cost = c(1, 5, 60, 40, 1000, 1))

print(cost)

The link set is created based on these assumptions with the graphab_link function, which takes the main following arguments:

See more details in the help file (?graphab_link).

Importantly, users can control the topology of the created link set, which will constrain the topology of the graphs subsequently created, and can substantially fasten the computations. This can be achieved by specifying whether the link set has a complete topology (topo = complete) in which all patches are connected to every other patches ($n\times(n-1)/2$ links), or instead a minimum planar topology (topo = complete) in which only neighbour patches are connected (though indirect paths can be then computed from this subset of links). Similarly, users can set a maximum cost (maxcost), which is the accumulated cost of a least-cost path beyond which the computation stops.

Since Graphab 3.0, when there are multiple habitat types in the project, users have to specify which ones are considered in the link set computation. If habitat = 'all', then all existing habitat patches can be connected by links, regardless of their nature. Otherwise, the name of one or several habitat types has to be provided. Finally, when several habitat types are specified, it is now possible to compute only the links connecting a patch of one type to a patch of another type, with the following argument: inter = TRUE.

In our example, we first create a planar link set using cost-distances to connect forest patches:

graphab_link(proj_name = proj_name,
             distance = "cost",
             name = "forest_link_planar",
             habitat = "forest",
             cost = cost,
             topo = "planar")

We can create a similar link set for shrubland patches:

graphab_link(proj_name = proj_name,
             distance = "cost",
             name = "shrub_link_planar",
             habitat = "shrubland",
             cost = cost,
             topo = "planar")

To create a link set connecting only forest patches to shrubland patches (and reciprocally), we can select these two habitat types for the computation and add inter = TRUE, as follows:

graphab_link(proj_name = proj_name,
             distance = "cost",
             name = "inter_forest_shrub",
             habitat = c("forest", "shrubland"),
             inter = TRUE,
             cost = cost,
             topo = "planar")

Note that we would have obtained a link set combining the links of the three link sets created above by computing a single link set for the two habitat types (habitat = c("forest", "shrubland")), including links both among patches of each habitat type and among patches of different types (inter = FALSE). The planar topology would however have been computed differently in the latter case, leading to small differences in the results.

We can load the cost-distance values resulting from these computations in R using the get_graphab_linkset function. The result is returned as a data.frame.

link_forest <- get_graphab_linkset(proj_name = proj_name,
                                   linkset = "forest_link_planar")
print(link_forest[1:6, ])

link_forest <- get_graphab_linkset(proj_name = proj_name,
                                   linkset = "forest_link_planar",
                                   proj_path = system.file('extdata', 
                                            package = 'graph4lg'))
print(link_forest[1:6, ])

Similarly, we can recover the parameters used to create a given link set with the get_graphab_linkset_cost or graphab_project_desc functions.

link_param <- get_graphab_linkset_cost(proj_name = proj_name,
                                      linkset = "forest_link_planar")

print(link_param)

link_param <- get_graphab_linkset_cost(proj_name = proj_name,
                                       proj_path = system.file('extdata', 
                                            package = 'graph4lg'),
                    linkset = "forest_link_planar")

print(link_param)

The graphab_project_desc function provides information about all the elements of the project:

graphab_project_desc(proj_name = proj_name)

graphab_project_desc(proj_name = proj_name,
                     proj_path = system.file('extdata', 
                                            package = 'graph4lg'))

Graph construction

Once we have created at least one habitat and a link set among its patches, we can create a habitat graph (or landscape graph) with graphab_graph. This function takes as arguments:

proj_name is the only mandatory argument. Without any other argument, a non-thresholded graph is created for every link set present in the project. Names are created automatically.

Since a link set is necessarily associated with one or several habitat types, the nodes of the graph will be determined from the link set properties.

In our example, we can create a planar graph with the following command:

graphab_graph(proj_name = proj_name,
              linkset = "forest_link_planar",
              name = "graph_forest")

If we select a linkset linking multiple habitat types, the graph will have nodes corresponding to multiple types of habitats. In the present case, since the only linkset connecting multiple habitat types only includes inter-links among patches of different types, if we want to include both intra- and inter-links in the graph, we can create 1 graph for each of the linksets (2 intra and 1 inter) and then merge the 3 graphs into one using the graphab_merge_graph function (created for that purpose). We would proceed as follows:

# Forest graph
graphab_graph(proj_name = proj_name,
              linkset = "forest_link_planar",
              name = "graph_forest_1")
# Shrub graph
graphab_graph(proj_name = proj_name,
              linkset = "shrub_link_planar",
              name = "graph_shrub_2")
# Inter graph
graphab_graph(proj_name = proj_name,
              linkset = "inter_forest_shrub",
              name = "graph_inter_3")
# Merge the 3 graphs
graphab_merge_graph(proj_name = proj_name,
              name = "graph_multi",
              graphs = c("graph_forest_1",
                         "graph_shrub_2",
                         "graph_inter_3"))

A simpler option consists in creating a linkset connecting multiple habitat types (inter=FALSE, see above) and then a graph using this linkset.

Landscape graph analysis

Once the graph has been created, it can be analysed by computing graph-theoretic connectivity metrics or by partitioning its nodes, among other possibilities.

Metric calculation

Many connectivity metrics can be computed from a landscape graph (see @baranyi2011contribution and @rayfield2011connectivity among other reviews on the subject). The Graphab software includes a large range of connectivity metrics and its manual provides users with a comprehensive description of every one of them (see Graphab 3.0 manual).

The graphab_metric function computes these metrics. It takes as arguments:

Metrics fall into different categories. PC and IIC are global metrics and take only one value for the entire graph, which is returned in R environment when return_val=TRUE.

The other metrics are computed at the node level. When return_val=TRUE, a data.frame is returned in R environment specifying the value of the metric for every graph node.

Importantly, the multihab argument now (since Graphab 3.0) allows users to decompose the value of the metrics into several components when the graph considered for the computation includes several types of habitats. If multihab='all', all the pairwise combinations are considered, including the within-habitat case. If multihab='inter', only the inter-habitat types combinations are considered. Note that this argument is only required if the graph on which you compute the metrics is based on several habitat types. Please find details in @savary2024multiple.

The description of every metric is beyond the scope of this tutorial. We again invite users to read the help file (?graphab_metric) as well as the Graphab 3.0 user manual.

We will compute metrics on the graph "graph_forest" in our example. First, we will compute the probability of connectivity at the global level. We set the dist and prob parameter such that: $p(10km)=e^{-\alpha \times d_{ij}}=0.05$. We convert 10 km in cost-distance units for the computation.

# Global metric: PC
pc <- graphab_metric(proj_name = proj_name,
               graph = "graph_forest",
               metric = "PC",
               dist = 10000,
               prob = 0.05,
               beta = 1,
               cost_conv = TRUE)
pc

res_g[["PC"]]

Note that we could have specified with the argument resfile the name of the .txt file storing the result of the global metric value.

We obtain the value in an object of class list.

Now, using the same parameters, we compute the local metric Flux.

f <- graphab_metric(proj_name = proj_name,
                    graph = "graph_forest",
                    metric = "F",
                    dist = 10000,
                    prob = 0.05,
                    beta = 1,
                    cost_conv = FALSE)

print(res_g[["F"]][1:6, ])

Every time a local metric is computed, it can be returned in R environment (here stored in f) but is also stored in the habitat patch attribute table, such that at the end this table contains every metric computed. Because of this, we can load the metric values in R later on, using the get_graphab_metric function. For instance:

metric_table <- get_graphab_metric(proj_name = proj_name,
                   graph = "graph_forest")

Metrics can be filtered based on their specific name (one metric), the graph on which they were computed or the habitat types considered. See ?get_graphab_metric for more details.

In the specific case of multiple habitat graphs, we can compute different components of the metrics, each corresponding to the contribution of a specific habitat type. We add the multihab argument for that purpose. For instance with the EC metric computed on the multiple habitat graph graph_multi:

ec_multi <- graphab_metric(proj_name = proj_name,
                     graph = "graph_multi",
                     multihab = "all",
                     metric = "EC",
                     dist = 10000,
                     prob = 0.05,
                     beta = 1,
                     cost_conv = TRUE)
ec_multi

res_g[["EC_multi"]]

The X0.0 column indicates the contribution of the forest habitat nodes (habitat code 0) and the connections among them to the EC component (intra), whereas the X0.1 column indicates the contribution of the connections among forest and shrubland patches (inter), and X1.1 that of the shrubland patches and the connections among them (intra). Please see details about the formulas used for the decomposition in @savary2024multiple.

Similarly, in the case of local metrics, the contribution of each habitat type to the local connectivity metric computed for each node can be assessed. For the Flux metric for instance, we can know how much a node is connected to the nodes of each other habitat type:

f_multi <- graphab_metric(proj_name = proj_name,
                    graph = "graph_multi",
                    multihab = "all",
                    metric = "F",
                    dist = 10000,
                    prob = 0.05,
                    beta = 1,
                    cost_conv = FALSE)
f_multi

print(res_g[["F_multi"]][1:6, ])

Graph partitioning

Another way to analyse landscape graph connectivity and topology is to make a partitioning. The principle is the same as that of modularity analyses of genetic graphs described in the second tutorial.

The graphab_modul function carries out such an analysis, but for the moment, it only creates a shapefile polygon layer with the Voronoi polygons corresponding to the groups of patches forming modules.

It takes as arguments:

For example:

graphab_modul(proj_name = proj_name,
              graph = "graph_forest",
              dist = 10000,
              prob = 0.05,
              beta = 1)

NOTA BENE: All these analyses can be performed for several graphs created in the same project, provided their name is correctly indicated in the function arguments.

Using landscape graphs in landscape genetic analyses

Landscape graphs have been created and analysed by computing connectivity metrics and/or by partitioning them. The output of these analyses can be used in landscape genetic analyses, e.g., to assess the relationship between genetic variables measured in populations inhabiting the study landscape and the connectivity of its habitat patches. To do that, we can import Graphab output in R environment to use it in subsequent (statistical) analyses.

As seen before, users can import in the R environment the link set properties in a data.frame with an edge list format with the function get_graphab_linkset which takes as arguments the project name and the link set name.

get_graphab_linkset(proj_name = proj_name,
                    linkset = "forest_link_planar")

print(res_g[["LK"]][1:6, ])

Similarly, users can import the table in which are stored all the metrics computed in the project as well as the node properties (including their areas). The function get_graphab_metric takes as argument the project name. For instance for the metrics relative to the forest habitat: 

get_graphab_metric(proj_name = proj_name,
                   habitat = "forest")

print(res_g[["MET"]][1:6, ])

Finally, users can also create a graph from a Graphab link set and a given habitat type, and convert it into a graph object of class igraph. The imported graph has weighted links. Habitat patch attributes present in the Graphab project are included. A figure can be displayed automatically representing the graph on a map. The graph is given the topology of the selected link set.

The function graphab_to_igraph takes as arguments:

land_graph <- graphab_to_igraph(proj_name = proj_name,
                                linkset = "forest_link_planar",
                                habitat = "forest",
                                weight = "cost",
                                fig = FALSE,
                                crds = TRUE)

crds_patches <- land_graph[[2]]
land_graph <- land_graph[[1]]

crds_patches <- res_g[["CRDS"]]
land_graph <- res_g[["LGRAPH"]]

The function returns a list of two objects:

This graph can be plotted on a map with plot_graph_lg with node sizes proportional to habitat patch area and link width inversely proportional to cost-distances:

plot_graph_lg(land_graph,
              crds = crds_patches,
              mode = "spatial",
              node_size = "area")

Computational performance

Some computations involved in the construction and analysis of landscape graphs can be heavy and take a lot of time and computational resources on your machines. Two main arguments are now included in most of the functions calling Graphab to adjust the RAM allocated to each computation, and the number of cores on which the computation can be run in parallel. These arguments are the following:

Conclusion

We have seen in this third vignette how to construct, analyse, and import landscape graphs with Graphab and graph4lg. In the next and last vignette, we will see how to compare genetic graphs and landscape graphs.

Note that other functions of graph4lg call Graphab and can be used to create corridors along least-cost paths, interpolate metrics, add pointsets to a project, among other examples.

References



Try the graph4lg package in your browser

Any scripts or data that you put into this service are public.

graph4lg documentation built on Sept. 24, 2026, 5:12 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close