Definitions

• Raster - Raster is a digital encoding of 2D digital image with underlying pixels as underlying. It also encapsulates other details like resolution, dimensions etc which is useful in identification and processing. Here, we use raster to represent maps.

• TIFF(Tag Image File Format) is a file that stores raster and information.

External materials

• Download and install R, R Studio - https://rstudio-education.github.io/hopr/starting.html/

• Raster - https://www.adobe.com/creativecloud/file-types/image/raster.html

• TIFF - https://www.adobe.com/creativecloud/file-types/image/raster/tiff-file.html/

Data sets

We use publicly available sources to obtain crop information -

• Monfreda (Monfreda, Ramankutty, and Foley 2008)

• MapSpam (2010, 2017) (International Food Policy Research Institute 2019)

Geodata provides set of APIs to access these data-sets. For visualization and plots, we use rnaturalearth.

Meta information

• Package information - https://garrettlab.github.io/HabitatConnectivity/
• CRAN site - https://CRAN.R-project.org/package=geohabnet
• Source code - https://github.com/GarrettLab/HabitatConnectivity/tree/main/geohabnet/
• Report issues - https://github.com/GarrettLab/CroplandConnectivity/issues/
• Lab website - https://www.garrettlab.com/

Installation and loading

Installing geohabnet will also install its dependencies. Please see the list of dependencies using desc::desc(package = "geohabnet")

if (!require("devtools")) {
  install.packages("devtools")
}

if (!require("geohabnet")) {
  utils::install.packages("geohabnet")
}

Alternatively, install from github source.

library("devtools")
#> Loading required package: usethis

if (!require("geohabnet")) {
  install_github("GarrettLab/HabitatConnectivity", subdir = "geohabnet")
}

library(geohabnet)

At any point, access the help page using following -

?geohabnet
# and
?geohabnet::msean # any function

Run analysis on default configuration

This is to run the analysis on default set of values for the supported parameters. Initially, the values are set based on the Xing et al(2021) (Xing et al. 2020) and crop is Avocado. This would run the workflow on global geographical extent quickly since the crop presence is relatively low.

res <- sensitivity_analysis()
#> |---------|---------|---------|---------|==========                                          |---------|---------|---------|---------|==============                                          |---------|---------|---------|---------|==========                                          
#> Warning in value[[3L]](cond): beep() could not play the sound due to the following error:
#> Error in play.default(x, rate, ...): no audio drivers are available

The results are the side effects and the similar information is captured in the object returned from it. Its based on the values set for each of the parameters which we will see in details in the further sections.

Running new analysis

We provide 2 methods or entry points to run analysis by setting new values for the supported parameters.

avocado <- cropharvest_rast("avocado", "monfreda")

# verify the raster object
avocado
#> class       : SpatRaster 
#> dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
#> resolution  : 0.08333333, 0.08333333  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : avocado_HarvestedAreaFraction.tif 
#> name        : avocado_HarvestedAreaFraction
gplot(avocado)

geo_net <- msean(avocado)
#> |---------|---------|---------|---------|=======                                          |---------|---------|---------|---------|==========                                          |---------|---------|---------|---------|=======                                          

Using configuration

More parameters are available under configuration and thus more control over the analysis. The configuration file name is parameters.yaml, currently supporting up to 10 parameters. The intuition behind this methodology is to provide a basic interface for setting new values. The snippet below describes the basic usage of configuration.

Get the initial configuration file. By default this function will save the file in temporary directory tempdir(), however we recommend saving to path where program will have write permissions. Using iwindow = TRUE will prompt a selection window to save config file.

config_file <- get_parameters(out_path = tempdir())
config_file
#> [1] "/tmp/RtmpH57uOF/parameters.yaml"

The file should look something like this -

Initial parameters.yaml

The values must be modified without modifying the structure. The order don’t matter for the program. The new values in the configuration must be fed to the workflow using -

set_parameters(new_params = config_file)
#> [1] TRUE
#using iwindow = true will prompt a selction window to choose config file.

get_parameters() was only to fetch the initial parameters. While you can, it is not required to re-fetch if the parameters has not been modified in the configuration. Modify the value and feed it to workflow using set_parameters().

Hosts

We have some helper functions to make the analysis easier. Available sources can be seen using -

geohabnet::supported_sources()
#> [1] "monfreda"          "mapspam2010"       "mapspam2017Africa"

More specifically, we support 2 data set sources - Monfreda and Mapspam using the APIs provided by Geodata. These sources have crop presence in the geographical areas in the form of raster(tiff). To start with, we can check if a certain crop is available in the source -

search_crop("avocado") #1
#> [1] "monfreda"
search_crop("bean") #2
#> [1] "monfreda" "spam"

The #1 result suggests that we have crop information available under monfreda . Now, #2 suggests that banana is present in 2 sources. We will use the results from here to obtain its raster data from the one of the specified sources.

# get avocado data
rast_avocado <- crops_rast(list(monfreda = "avocado"))

old_to <- options("timeout")
options(timeout = 6000)

#get data of banana crop
rast_ban <- crops_rast(list(mapspam2010 = "banana"))

The returned is actually if type SpatRaster. In above, we will have separate raster information for each of the crop. Crops can be combined and be returned as one raster.

# Back to back downloads causes connection to fail, this is a workaround
rast_avo_ban <- crops_rast(list(monfreda = c("avocado", "banana"), mapspam2010 = c("banana")))
old_to <- options(old_to)

In this case for banana, information from all the sources will be first combined using mean and then all the crop information i.e. banana and avocado both are again combined into one raster. This is done to produce a single map which contains information of all the requested crops. We have one more method to obtain raster, using tiff - tiff_torast() . Finally, we have get_rasters() method, which is an abstraction of . To see the difference in data that was generated in last 3 methods -

gplot(c(rast_avocado, rast_ban, rast_avo_ban))

The result obtained from one of these methods can be used as an argument to sean() . Value of results is list of risk indices.

results <- lapply(list(rast_avocado, rast_avo_ban), msean)

Internally, cropharvest_rast() is used to fetch the crop data-set in it’s supported form.

So far, we have now run the sensitivity analysis workflow for avocado and banana on global scale and default set of parameter values. Alternatively, we can just specify the crops under parameters.yaml and run the workflow -

results <- sensitivity_analysis()
#> |---------|---------|---------|---------|=====                                          |---------|---------|---------|---------|=======                                          |---------|---------|---------|---------|==============                                          |---------|---------|---------|---------|=====                                          
#> Warning in value[[3L]](cond): beep() could not play the sound due to the following error:
#> Error in play.default(x, rate, ...): no audio drivers are available

Thresholds

Thresholds are used to select subset of values from the SpatRaster on which the operations are applied. It directly effects the connectivity and gives a sense of sensitivity in the network. The intermediate goal is to produce a adjacency graph which essentially determines the connectivity. Cells which doesn’t meet the threshold are removed from the consideration by dispersal models.

risk_indexes <- msean(avocado, global = FALSE, hd_threshold = 0.00001, link_threshold = 0.00001)
#> |---------|---------|---------|---------|==                                          |---------|---------|---------|---------|==                                          |---------|---------|---------|---------|=                                          

Aggregation

Aggregation strategy refers to the function used to create a new map of host density with a lower resolution (larger cells). Reducing the spatial resolution helps to reduce the computational power needed to run the analysis.

• If AggregationStrategy: [sum], then the sum of the cropland proportion of all initially small grids within a large grid is divided by the total number of initially small grids within that large grid.

• If AggregationStrategy: [mean], then the sum of the cropland proportion of all initially small grids within a large grid is divided by the total number of initially small grids containing only land (where small grids with water are excluded) within that large grid.

By default, analysis is run on both but can be opted out from one. If only one method is used, then the difference map is skipped from the outcome.

Distance methods

For each pair of locations in the host map with values greater than the host density threshold, the sensitivity_analysis() function will calculate the physical distances and use them to calculate the relative likelihood of pathogen movement between locations based on their pairwise geographical proximity.

There are two different options to calculate the distance between locations.

· Vincenty ellipsoid distance

This option is highly accurate but more computationally expensive.

· Geodesic distance

This option is less computationally expensive and less accurate than the option above.

You can set the distance option either as DistanceStrategy: “vincentyEllipsoid” or DistanceStrategy: “geodesic”. One of these options should be used as input to run the analysis. Check for supported methods in analysis by running dist_methods() in the console.

dist_methods()
#> [1] "geodesic"          "vincentyellipsoid"

Resolution

The aggregation factor or granularity is the number of small grid cells that are aggregated into larger grid cells in each direction (horizontally and vertically). The finest value is 1 which can require analysis to run up to hours because of large number of cells in SpatRaster . The resolution is also used in calculation of variance while dis-aggregating the risk indices into coarser resolution for producing maps.

If not provided, the defaulted value is selected from reso()

Metrics

See available metrics using

supported_metrics()
#> [1] "betweeness"               "node_strength"           
#> [3] "sum_of_nearest_neighbors" "eigenvector_centrality"  
#> [5] "closeness"                "degree"                  
#> [7] "page_rank"

Metrics corresponding to dispersal models are applied to distance matrix with specified weights. The weights must be specified in % and sum of all the weights should be equal to 100. We use functions from (Csardi and Nepusz 2006) to calculate metrics for each dispersal model. The 2 dispersal models that are applied to parameters inverse power law and negative exponential. More formally, metrics are way to determine connectivity among nodes in a network.

In a graph functions of (Csardi and Nepusz 2006), the links are interpreted as distances. However, in the context of habitat connectivity, the network is adapted to interpret links as weights which means that the likelihood of pathogen spread is lower if the distance is larger.$$ L = \frac{1}{W}, \\ W' = \sum_{i=1}^{N} \max(W - W_i) $$

L is link weights and W is the original weights in an undirected graph. W’ is the transformed weight vector for calculating network centrality.

Geographical Extent

Geographical extent is a subset of world map defined by coordinate reference system. The corresponding parameter to set the area in sean() and sensitivity_analyis() is geoscale and GeoExtent respectively. Default setting is global = TRUE which will ignore the value of geoscale. This will consider taking the world map into account using values from global_scales() . For non-global analysis, either set global = FALSE with or without the setting geoscale. By default, geoscale will be extracted from the extents from the input SpatRaster. We recommend using EPSG:4326 as coordinate reference system because the functions have been tested on it.

1. Using function

   results <- msean(avocado, global = FALSE, geoscale = c(-115, -75, 5, 32))

2. Using config

   Set Global = FALSE and CustomExt = [-115, -75, 5, 32] . The initial parameters.yaml already contains this value which would run in combination with other parameters.

When provided with geoscale, program will take the subset of provided raster (data-set of a crop). The workflow will apply graph operations and network connectivity only to the subset.