What is Graphab?

Graphab is an open source software tool dedicated to the modelling of landscape networks (Foltête et al. 2012, 2021). 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. Foltête et al. (2012) and Foltête et al. (2021) 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 Savary et al. (2024) 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:

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 Foltête et al. (2021) and Savary et al. (2024) 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.

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:

• proj_name: Name of the project (project and .xml file directory)
• raster: Path to the raster file on which habitat patches are identified
• nodata (optional): Cell value corresponding to no data values, thereafter ignored
• parallel.java: An integer indicating how many computer cores are used to run the .jar file. By default, parallel.java = NULL, and java sets it according to local settings.
• alloc_ram (optional): Integer or numeric value indicating RAM gigabytes allocated to the java process. Increasing this value can speed up the computations. Too large values may not be compatible with user’s machine settings
• proj_path(optional): Path to the directory that contains the project directory. It should be used when the project directory is not in the current working directory. Default is NULL. When 'proj_path = NULL', the project directory is equal to getwd().

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")

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:

• proj_name: Name of the project (project and .xml file directory).
• name: The name of the new habitat type (character string), arbitrarily defined by the user.
• type: A character string specifying whether the new habitat type is defined from (i) cells of the raster used to create the project ('raster') or (ii) an external vector spatial layer ('vector', can be points or polygons).
• rast_codes: Values of the raster cells to consider to create the new habitat type when type = 'raster'.
• vec_layer: Path to the vector spatial layer used to create the new habitat type when type = 'vector'.
• vec_capa_field: Column/field of the vector layer attribute table used to set the capacity of the habitat patches created when type = 'vector'.
• minarea (optional): Minimum habitat patch size in hectares.
• maxsize: (optional, default=NULL, only if ‘type=“raster”’). An integer or numeric value specifying the maximum side length of the rectangular full extent of each habitat patch in metric units. If this side length exceeds maxsize m, then several patches are created.
• con8: (optional, default=FALSE) A logical indicating whether a neighborhood of 8 pixels (TRUE) is used for patch definition. By default, con8=4, corresponding to 4 pixel neighborhood.
• parallel.java: An integer indicating how many computer cores are used to run the .jar file. By default, parallel.java = NULL, and java sets it according to local settings.
• alloc_ram (optional): Integer or numeric value indicating RAM gigabytes allocated to the java process. Increasing this value can speed up the computations. Too large values may not be compatible with user’s machine settings
• proj_path(optional): Path to the directory that contains the project directory. It should be used when the project directory is not in the current working directory. Default is NULL. When 'proj_path = NULL', the project directory is equal to getwd().

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")

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:

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)
#>   code cost
#> 1    0    1
#> 2    1    5
#> 3    2   60
#> 4    3   40
#> 5    4 1000
#> 6    5    1

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

• proj_name: Name of the project
• distance: Type of distance computed. By default, distance="cost", meaning that cost distances are computed between patches according to the cost values specified in the cost data.frame. Alternatively, distance="euclid", meaning that straight line Euclidean distances are computed.
• name: A character string indicating the name of link set.
• habitat: (optional, default=NULL) A character string indicating the name of the habitat types whose patches will be connected by the least-cost paths. If name='all', all existing habitats are considered. If a single habitat is indicated (e.g., habitat = 'forest'), then only this one is used. If several habitats are indicated, they are merged before the computation (e.g., habitat = c('forest', 'grass')). If habitat = NULL, each is separately considered.
• cost: A data.frame with the cost values (cost in our example).
• inter: (optional, default = FALSE) A logical indicating whether only links among different habitat types are computed (if TRUE, no “intra-type” links).
• topo: A character string indicating the topology of the created link set. It can be planar (topo='planar' (default)). It speeds up the computation but will prevent from creating complete graphs with graphab_graph, or it can be complete (topo='complete').

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, ])

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)

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

graphab_project_desc(proj_name = proj_name)

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: Name of the project.
• linkset: Name of the link set used to create the graph.
• name: A character string indicating the name of the graph created, if only one is created.
• thr: An integer or numeric value indicating the maximum distance associated with the links of the created graph. It allows users to create a pruned graph based on a distance threshold. Note that when the link set used has a planar topology, the graph is necessarily a pruned graph (not complete) and adding this threshold parameter can remove other links. When the link set has been created with cost-distances, the parameter is expressed in cost-distance units whereas when the link set is based upon Euclidean distances, the parameter is expressed in meters.
• cost_conv: Logical (TRUE or FALSE) indicating whether numeric thr values are converted from cost-distance into Euclidean distance using a log-log linear regression. See also convert_cd function presented in the first tutorial.

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.

Metric calculation

Many connectivity metrics can be computed from a landscape graph (see Baranyi et al. (2011) and Rayfield et al. (2011) 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:

• proj_name: Name of the project

• graph: Name of the graph on which the metric is computed

• metric: Name of the metric among:

• Probability of Connectivity ("PC"),

• Integral Index of Connectivity ("IIC"),

• Flux ("F"),

• Interaction Flux ("IF"),

• Degree ("Dg"),

• Clustering Coefficient ("CCe"),

• Current Flow ("CF")

• delta Probability of Connectivity ("dPC").

• dist: A numeric or integer value specifying the distance at which dispersal probability is equal to prob. This argument is mandatory for weighted metrics (PC, F, IF, BC, dPC, CCe, CF) but not used for others. It is used to set \(\alpha\) for computing dispersal probabilities associated with all inter-patch distances such that dispersal probability between patches \(i\) and \(j\) is \(p_{ij}= e^{-\alpha d_{ij}}\).

• prob: A numeric or integer value specifying the dispersal probability at distance dist. By default, code=0.05. It is used to set \(\alpha\) (see argument dist above).

• beta: A numeric or integer value between 0 and 1 specifying the exponent associated with patch areas in the computation of metrics weighted by patch area. By default, beta=1. When beta=0, patch areas do not have any influence in the computation.

• multihab: (optional, default = NULL) A character string indicating whether the metric value should be decomposed across the different habitat types and their pairwise combinations.

• cost_conv: Logical (TRUE or FALSE) indicating whether numeric dist values are converted from cost-distance into Euclidean distance using a log-log linear regression.

• return_val: Logical (default = TRUE) indicating whether metric values are returned in R (TRUE) or only stored in the patch attribute layer (FALSE)

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 Savary et al. (2024).

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

#> $`Metric name`
#> [1] "pc_d1479.8358851873484_p0.05_graph_forest"
#> 
#> $`Metric value table`
#>          Graph        d    p beta           PC
#> 1 graph_forest 1479.836 0.05    1 0.0007436003

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)

#>   habitat id_habitat id_patch    area perim capacity
#> 1  forest          0        1 2610000 13000  2610000
#> 2  forest          0        2 8800000 23000  8800000
#> 3  forest          0        3 3270000 10800  3270000
#> 4  forest          0        4 5330000 18800  5330000
#> 5  forest          0        5 6580000 18200  6580000
#> 6  forest          0        6 2620000  8800  2620000
#>   f_d1479.8358851873484_p0.05_beta1_graph_forest
#> 1                                     22659462.5
#> 2                                     12467998.8
#> 3                                      2264071.3
#> 4                                       340302.5
#> 5                                     13018161.5
#> 6                                     16655121.0

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

#> $`Metric name`
#> [1] "EC_graph_multi"
#> 
#> $`Metric value table`
#>         Graph       d    p beta     X0.0     X0.1     X1.1
#> 1 graph_multi 813.815 0.05    1 82578203 66570524 29719307

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 Savary et al. (2024).

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

#>   habitat id_habitat id_patch    area perim capacity
#> 1  forest          0        1 2610000 13000  2610000
#> 2  forest          0        2 8800000 23000  8800000
#> 3  forest          0        3 3270000 10800  3270000
#> 4  forest          0        4 5330000 18800  5330000
#> 5  forest          0        5 6580000 18200  6580000
#> 6  forest          0        6 2620000  8800  2620000
#>   f_d10000_p0.05_beta1.0_graph_multi f_d10000_p0.05_beta1.1_graph_multi
#> 1                          117669295                           49122576
#> 2                          105967219                           46842866
#> 3                           64066487                           35918701
#> 4                           42748208                           29043326
#> 5                          119935793                           51041965
#> 6                          123382351                           50838307

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:

• proj_name: Name of the project
• graph: Name of the graph on which the partition is performed
• dist: A numeric or integer value specifying the distance at which dispersal probability is equal to prob. This argument is used to weight the links when modularity is computed.
• prob: A numeric or integer value specifying the dispersal probability at distance dist. By default, code=0.05. It is used to set \(\alpha\) (see argument dist above).
• beta: A numeric or integer value between 0 and 1 specifying the exponent associated with patch areas in the computation of the modularity. By default, beta=1. When beta=0, patch areas do not have any influence in the computation.
• nb: Optional argument indicating the number of modules to create. By default, the number of modules maximising the modularity index is used.

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.