Definitions

• Raster - Raster is a digital encoding representation ofgeographic information as a grid of cells or pixels, where each cellcontains a value representing a specific property or attribute of alocation. Some details of a raster data include spatial resolution,dimensions, or geographic extent, which are useful in identification andprocessing. Here, we use raster to represent maps.

• TIFF(Tag Image File Format) is a type of file that storesgeographic information in the form of raster.(Keshav et al. 2023) has been tested using.tif and.nc files but it would be able to handle anyraster type of file supported by the R packageterra(Hijmans 2023) .

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

For estimating habitat connectivity for understanding the spread ofplant pathogens or pests, we use cropland density or host availabilityaround the world as a proxy for habitat availability (this assumptioncan be complemented with species-specific environmental suitability fora better proxy of habitat availability). We use publicly availablesources to obtain crop distribution in the form of raster layers.

• CROPGRIDS(Tang FHM 2024)

• EARTHSTAT(Monfreda, Ramankutty, and Foley2008)

• MapSpam(2010, 2017)(International Food PolicyResearch Institute 2019)

• CroplandCROS(Service 2025)

• GBIF(GBIF 2025)

geodata provides a set of APIs to access the EARTHSTAT andMapSPAM datasets directly in R (see code below).

The raster layers from the data sources listed above may need sometransformation before using them in(Keshav etal. 2023). A valid raster layer to be used is one that its entryvalues of each pixel range from 0 (no habitat available) to 1 (habitatis available in the entirety of the pixel). For example, data layers ofcropland fraction from EARTHSTAT can be directly provided to functionswithout the need to re-scale its entries. However, data layers ofharvested area from CROPGRIDS will need re-scaling the values of eachgrid cell in the raster (see an example of how to do this below).

For visualization and plots,geohabnet usesrnaturalearthrnaturalearth.

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/
• Garrett lab website -https://www.garrettlab.com/

Installation and loading

Installinggeohabnet will automatically install itsdependencies on other R packages. See the list of dependencies usingdesc::desc(package = "geohabnet")

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

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

Alternatively, you can install the development version ofgeohabnet from GitHub repository.

library("devtools")#> Loading required package: usethisif(!require("geohabnet")){install_github("GarrettLab/HabitatConnectivity", subdir="geohabnet")}#> Loading required package: geohabnet

Either way, now you can loadgeohabnet in R asfollows:

library(geohabnet)

At any point, you can access the help page using the following

?geohabnet# For description of the package# and?geohabnet::msean# For documentation of a function (in this case, the msean function)

Habitat connectivity analysis based on default parameters

The section below illustrate how the habitat connectivity analysiscan be conducted using the default set of values for the supportedparameters. For this example, we use the raster layer for avocadoharvested area from CROPGRIDS(Tang FHM2024) as a real-world case study for the habitat connectivity ofavocado-specific pathogens. This example would run the analysis onglobal geographical extent relatively quickly since the crop presencearound the world is relatively low.

For the purposes of demonstration, the code below sets the hostparameter programatically or directly in RStudio. Alternatively, theparameters can be fed viaparameters.yaml in an interactiveway using either the path or settingiwindow = TRUE for thewindows prompt inset_parameters()

Running your first habitat connectivity analysis

If you are also replicating this example, please first download theavocado layer from CROPGRIDS and put it in the folder you are currentlyworking. Note that we are downloading the raster layer for hostavailability (or habitat availability) manually from CROPGRIDS. However,an alternative option is getting the data layer directly in RStudiousing thegeodata package (depending on whether it is downor not). At the time, we created this vignettegeodata wasnot working on our local machine but the R code to do so isprovided.

#if (!require(geodata))#  install.packages("geodata")#avocado_sp <- geodata::crop_monfreda(crop = "avocado", var = "area_ha", path = tools::R_user_dir("geohabnet", which = "config"))library(terra)#> terra 1.8.50#avocado_rast_input <- system.file("CROPGRIDSv1.08_avocado.nc", package = "geohabnet", mustWork = TRUE)avocado_sp<-rast("CROPGRIDSv1.08_avocado.nc")cell.area<-(0.05*111111)*(0.05*111111)/10000#area in heactaresavocado_sp<-avocado_sp$harvarea/cell.area#area in hectaresvalues(avocado_sp)<-ifelse(values(avocado_sp)>0,values(avocado_sp),NaN)avocado_sp<-aggregate(avocado_sp, fact=2, fun="mean", na.rm=TRUE)#> |---------|---------|---------|---------|=========================================writeRaster(avocado_sp,"avocado_density.tif", overwrite=TRUE)

Note that, apart from getting the raster in R, we are also scalingthe entry values of each grid cell in the raster layer, from harvestedarea in hectares to avocado cropland fraction. We are aggregating thisraster layer to reduce computational cost when running the analysis (forillustrative purposes only, you may not need to do this for your owndata set). We are also saving the new avocado area fraction as a newraster layer in our local machines.

avocado_sp#> class       : SpatRaster#> dimensions  : 1800, 3600, 1  (nrow, ncol, nlyr)#> resolution  : 0.1, 0.1  (x, y)#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)#> coord. ref. :#> source(s)   : memory#> name        :     harvarea#> min value   : 3.268422e-06#> max value   : 1.700697e-01gplot(avocado_sp)

geo_net<-msean(avocado_sp)#> |---------|---------|---------|---------|=====#> Warning: [+] CRS do not match

#> |---------|---------|---------|---------|======#> Warning: [+] CRS do not match#> |---------|---------|---------|---------|==========#> Warning: [+] CRS do not match

#> |---------|---------|---------|---------|====#> Warning: [+] CRS do not match

Using configuration

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

Get the initial configuration file. By default this function willsave the file in temporary directorytempdir(), however werecommend saving to path where program will have write permissions.Usingiwindow = TRUE will prompt a selection window to saveconfig file.

config_file<-get_parameters(out_path=tempdir())config_file#> [1] "/var/folders/r5/zggvft9d3yn5kh51wqp78rd00000gn/T//RtmphRPeWc/parameters.yaml"

The file should look something like this -

Initial parameters.yaml

The values must be modified without modifying the structure. Theorder don’t matter for the program. The new values in the configurationmust be fed to the workflow. For the purposes of demonstration, the codebelow sets the host parameter programmitically. The parameters can befed viaparameters.yaml in an interactive way as well usingeither the path or settingiwindow = TRUE for the windowsprompt inset_parameters()

if(!require("geodata")){install.packages("geodata")}#> Loading required package: geodatapath_to_avo<-"avocado_density.tif"#terra::sources(avocado_sp)params_file<-geohabnet::get_parameters()if(!require(yaml))install.packages("yaml")#> Loading required package: yamlparams_yaml<-yaml::yaml.load_file(params_file)params_yaml$default$`CCRI parameters`$Host<-path_to_avoyaml::write_yaml(params_yaml,params_file)geohabnet::set_parameters(params_file)#> [1] TRUE#using iwindow = true will prompt a selection window to choose config file.

get_parameters() was only to fetch the initialparameters. While you can, it is not required to re-fetch if theparameters has not been modified in the configuration. Modify the valueand feed it to workflow usingset_parameters().

Hosts

The library supports a spatial raster file (e.g. TIFF). The casestudies were done on crop data from Monfreda(Monfreda, Ramankutty, and Foley 2008) andMapspam(International Food Policy ResearchInstitute 2019) which is available viageodatapackage(geodata2025?)

results<-sensitivity_analysis()#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> |---------|---------|---------|---------|====#> Warning: [+] CRS do not match

#> |---------|---------|---------|---------|======#> Warning: [+] CRS do not match#> |---------|---------|---------|---------|==========#> Warning: [+] CRS do not match

#> |---------|---------|---------|---------|====#> Warning: [+] CRS do not match

Thresholds

Thresholds are used to select subset of values from theSpatRaster on which the operations are applied. It directlyeffects the connectivity and gives a sense of sensitivity in thenetwork. The intermediate goal is to produce a adjacency graph whichessentially determines the connectivity. Cells which doesn’t meet thethreshold are removed from the consideration by dispersal models.

risk_indexes<-msean(avocado_sp, global=FALSE, hd_threshold=0.0001, link_threshold=0.001)#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> |---------|---------|---------|---------|==#> Warning: [+] CRS do not match

#> |---------|---------|---------|---------|==#> Warning: [+] CRS do not match

#> |---------|---------|---------|---------|==#> Warning: [+] CRS do not match

Aggregation

Aggregation strategy refers to the function used to create a new mapof host density with a lower resolution (larger cells). Reducing thespatial resolution helps to reduce the computational power needed to runthe analysis.

• If AggregationStrategy: [sum], then the sum of the croplandproportion of all initially small grids within a large grid is dividedby the total number of initially small grids within that largegrid.

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

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

Distance methods

For each pair of locations in the host map with values greater thanthe host density threshold, the sensitivity_analysis() function willcalculate the physical distances and use them to calculate the relativelikelihood of pathogen movement between locations based on theirpairwise geographical proximity.

There are two different options to calculate the distance betweenlocations.

· Vincenty ellipsoid distance

This option is highly accurate but more computationallyexpensive.

· Geodesic distance

This option is less computationally expensive and less accurate thanthe option above.

You can set the distance option either asDistanceStrategy:“vincentyEllipsoid” orDistanceStrategy: “geodesic”. Oneof these options should be used as input to run the analysis. Check forsupported methods in analysis by runningdist_methods() inthe console.

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

Resolution

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

If not provided, the defaulted value is selected fromreso()

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 distancematrix with specified weights. The weights must be specified in % andsum of all the weights should be equal to 100. We use functions from(Csardi and Nepusz 2006) to calculatemetrics for each dispersal model. The 2 dispersal models that areapplied to parametersinverse power law andnegativeexponential. More formally, metrics are way to determineconnectivity among nodes in a network.

In a graph functions of(Csardi and Nepusz2006), the links are interpreted as distances. However, in thecontext of habitat connectivity, the network is adapted to interpretlinks as weights which means that the likelihood of pathogen spread islower 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 undirectedgraph. W’ is the transformed weight vector for calculating networkcentrality.

Geographical Extent

Geographical extent is a subset of world map defined by coordinatereference system. The corresponding parameter to set the area insean() andsensitivity_analyis() isgeoscale andGeoExtent respectively. Defaultsetting isglobal = TRUE which will ignore the value ofgeoscale. This will consider taking the world map intoaccount using values fromglobal_scales() . For non-globalanalysis, either setglobal = FALSE with or without thesettinggeoscale. By default,geoscale will beextracted from the extents from the input SpatRaster. We recommend usingEPSG:4326 as coordinate reference system because the functions have beentested on it.

1. Using function

   results<-msean(avocado_sp, global=FALSE, geoscale=c(-115,-75,5,32))#> Warning: [+] CRS do not match

   

   #> Warning: [+] CRS do not match

   

   #> Warning: [+] CRS do not match

   

2. Using config

   SetGlobal = FALSE andCustomExt = [-115, -75, 5, 32] . The initialparameters.yaml already contains this value which would runin combination with other parameters.

When provided withgeoscale, program will take thesubset of provided raster (data-set of a crop). The workflow will applygraph operations and network connectivity only to the subset.

final<-msean(avocado_sp, link_threshold=0.0001, hd_threshold=0.0025)#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> Warning in max(wv, na.rm = TRUE): no non-missing arguments to max; returning#> -Inf#> |---------|---------|---------|---------|====#> Warning: [+] CRS do not match

# checkout the type of an objectclass(final)#> [1] "GeoNetwork"#> attr(,"package")#> [1] "geohabnet"

final#> class            :  GeoNetwork#> habitat density     :  98 , 267  (nrow, ncol) SpatRaster#> mean             :  /var/folders/r5/zggvft9d3yn5kh51wqp78rd00000gn/T//RtmphRPeWc/plots/mean_20251101192717060378.tif#> mean raster      :  1176 , 3204  (nrow, ncol)#> variance          :  /var/folders/r5/zggvft9d3yn5kh51wqp78rd00000gn/T//RtmphRPeWc/plots/variance_20251101192718031884.tif#> variance raster   :  1176 , 3204  (nrow, ncol)#> difference        :  /var/folders/r5/zggvft9d3yn5kh51wqp78rd00000gn/T//RtmphRPeWc/plots/difference_2025110119271871457.tif#> difference raster :  1176 , 3204  (nrow, ncol)

# checkout the resultsfinal@rasters#> class   :  GeoRasters#> global  :  TRUE#> globals :  1# global is TRUE because we original set the global analysis# thus, we will have set of 2 risk indices, eastern and wetern hemispherefinal@rasters$global_rast#> [[1]]#> class :  GlobalRast#> east  :  14#> west  :  14# Number of elements from above determines the the number of parameter values provided# To access the adjacency matrix,final@rasters$global_rast[[1]]$east[[1]]# this is also s4 class 'GeoModel'#> class            :  GeoModel#> adjacency matrix :  1 ,  1  (nrow, ncol)#> risk index       :  98 ,  170  (nrow, ncol)

# replace the indexing with any arbitary index,# uncomment line below to see the results.#final@rasters$global_rast[[1]]$east[[1]]$amatrix

#avocado <- geohabnet::cropharvest_rast("avocado", "monfreda")# see ?future::plan for detailsfuture::plan(future::multicore())msean(avocado_sp)future::plan(future::multisession())msean(avocado_sp)