Basic distance and connectivity functions in the Reconnect R-package
Demonstration of application of functions described in the CONN_functions.R and MPC_functions.R files.
#########################################################
### clean space
rm(list=ls(all=TRUE))
#########################################################
## set directories
## make sure you have loaded the "MSC_analysis.Rproj" correctly
if(!exists('sourcedir')){ sourcedir <- getwd() }; if(!exists('resdir')) { resdir <- getwd() }
##########################################################
## libraries
library(sf) # new spatial vector processing package
#> Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
library(raster) # old (but still useful) spatial raster processing
#> Lade nötiges Paket: sp
library(stars) # spatial processing (conversion between sf and raster possible, spatiotemporal datacube analyses)
#> Lade nötiges Paket: abind
library(terra)
#> terra 1.8.5
library(igraph)
#>
#> Attache Paket: 'igraph'
#> Die folgenden Objekte sind maskiert von 'package:terra':
#>
#> blocks, compare, union
#> Das folgende Objekt ist maskiert 'package:raster':
#>
#> union
#> Die folgenden Objekte sind maskiert von 'package:stats':
#>
#> decompose, spectrum
#> Das folgende Objekt ist maskiert 'package:base':
#>
#> union
library(Reconnect)
#>
#> Attache Paket: 'Reconnect'
#> Das folgende Objekt ist maskiert 'package:terra':
#>
#> add_legend
## some colors...
JOcols = c("cyan3","olivedrab3","olivedrab4")
JOcols1 = colorRampPalette(c("navy","lightblue","grey","mediumorchid","cyan3","cyan4","tan","olivedrab3","olivedrab2","olivedrab4"))Create a simple landscape
e.g. from the simple circle simulation in the NLMC_functions.R file
######################################################################
## set seed for reproducibility
set.seed(15)
### create a name of the simulation
## initialize parameters
nc1 = 100 ## define number of cells of landscape # 500 cells for dimx and dimy is reasonably good!!
hfr1 = c(0.15) ## define fraction of habitat in landscape, in combination with npatches, this defines patch area
npatches= c(3) ## define number of patches (maximum nr should not be more than 10 times the number of cells..) in combination with hfrs, this defines patch area
sdpa1 = 0 ## define sd of patch area
cf1 = 1 ## clumping factor: 1 = non-overlapping patches, >1: clumped & overlapping patches
remedge1= TRUE ## remove the edge of the landscape? - i.e. patches cannot go partially outside landscape
### prepare plot
par(mfrow=c(1,1))
### make a simple circle simulation
for(npatch1 in npatches){
res = simpcirc(dimx=nc1,dimy=nc1,hfr=hfr1,npatch=npatch1,sdpa=sdpa1,cf=cf1,form="circle",return="all",remedge=remedge1)
### unify new shape
sps1 = sf::st_union(res$sps,by_feature=FALSE) %>% st_cast("POLYGON") %>% st_sf # st_cast(sps1,"POLYGON")
## total new area
#tna = sum(res$rast[res$rast==1])
tna = sum(as.numeric(sf::st_area(sps1))) # new area
## total nr of new patches
tnpatch = length(sps1$geometry)
##plot result
plot(res$rast,main=sprintf("total area (ha) =%0.0f habitat fraction = %0.02f \n nr of patches = %0.0f clumping factor = %0.02f",tna,hfr1,npatch1,cf1),cex.main=0.8)
plot(sps1$geometry,add=TRUE,border="darkgreen")
}
#> [1] "creating patch areas with mean pa=500.0 and sdpa=0.0"
#> [1] "sample landscape for final number of patches = 3"
#> [1] "minimum distance = 25.2 units"
#> [1] 0
#> [1] "3 samples finished after 1 rounds"
#> [1] "finished strata 0 going to next..."
### from this simulation, you can derive habitat rasters and shapefiles
## habitat shapefile
lcshp = res$sps
plot(lcshp,main="habitat shapefile")
## for plotting later
ccoord = st_coordinates(st_centroid(lcshp))
#> Warning: st_centroid assumes attributes are constant over geometries
labels = lcshp$ID
## habitat raster
lcr = res$rast
#print(lcr)
plot(lcr,main="habitat landscape")
## resitance raster
lres = res$rast
## set all to 1 for the moment
lres[] = 1
## add a barrier whith resistance 10: cost distances should be 20 meters longer!
lres[5001:5099] = 20
#print(lres)
plot(lres,main="resistance landscape")
plot(lcshp$geometry,add=TRUE)
text(ccoord,labels=labels)
Apply some distance functions
e.g. by using functions in the CONN_functions.R file
### get euclidean distance among habitat patches!!
eucpm = eucdist_pm(x=lcshp,id="ID", dist_type = "centroid")
#> Warning: st_centroid assumes attributes are constant over geometries
print("euclidean distance, dist_type=centroid")
#> [1] "euclidean distance, dist_type=centroid"
print(eucpm)
#> $pa
#> [1] 499.9746 499.9746 499.9746
#>
#> $mdist
#> 1 2 3
#> 1 0.00000 50.44799 57.21888
#> 2 50.44799 0.00000 42.20190
#> 3 57.21888 42.20190 0.00000
eucpm = eucdist_pm(x=lcshp,id="ID", dist_type = "edge")
print("euclidean distance, dist_type=edge")
#> [1] "euclidean distance, dist_type=edge"
print(eucpm)
#> $pa
#> [1] 499.9746 499.9746 499.9746
#>
#> $mdist
#> 1 2 3
#> 1 0.00000 25.21750 31.98755
#> 2 25.21750 0.00000 16.97134
#> 3 31.98755 16.97134 0.00000
### get cost distance among habitat patches!!
costpm = costdist_pm(x=lcshp,id="ID", resi=lres, dist_type = "centroid")
print("cost distance, dist_type=centroid")
#> [1] "cost distance, dist_type=centroid"
print(costpm)
#> $pa
#> [1] 499.9746 499.9746 499.9746
#>
#> $mdist
#> 1 2 3
#> 1 0.00000 71.25483 59.07107
#> 2 71.25483 0.00000 63.35534
#> 3 59.07107 63.35534 0.00000
costpm = costdist_pm(x=lcshp,id="ID", resi=lres, dist_type = "edge")
print("cost distance, dist_type=edge")
#> [1] "cost distance, dist_type=edge"
print(costpm)
#> $pa
#> [1] 499.9746 499.9746 499.9746
#>
#> $mdist
#> 1 2 3
#> 1 0.00000 44.62742 32.00000
#> 2 44.62742 0.00000 36.14214
#> 3 32.00000 36.14214 0.00000
### get euclidean distance for each pixel in landscape!
euclm = eucdist_lm(x=lcr,exclude=NULL,dist_type = "edge")
## create a null raster to fill in values
lcr0 = lcr
lcr0[] = euclm[,1]
plot(lcr0, main="euclidean shortest \n distance to habitat")
plot(lcshp$geometry,add=TRUE)
text(ccoord,labels=labels)
### get euclidean distance from habitat centers to edge of habitat!
euchab = eucdist_hab(x=lcr,exclude=NULL)
lcr0 = lcr
lcr0[] = euchab[,1]
plot(lcr0, main="euclidean shortest \n distance to habitat edge")
plot(lcshp$geometry,add=TRUE)
text(ccoord,labels=labels)
### get cost distance at landscape scale!!
costlm = costdist_lm(x=lcr,resi=lres,scale=1,maxiter=100,dist_type="edge",cost_type ="least-cost")
lcr0 = lcr
lcr0[] = costlm[,1]
plot(lcr0, main="resistance shortest \n distance to habitat")
plot(lcshp$geometry,add=TRUE)
text(ccoord,labels=labels)
Apply some connectivity functions
e.g. by using functions in the CONN_functions.R & MPC_functions.R files
## graph centrality functions (betweenness centrality, closeness centrality and node degree from igraph)
centr = centr_igraph(pa=eucpm$pa,mdist=eucpm$mdist,alpha=40,dispfop="negex",weighted=TRUE,cutoffpr=0.01,MST=FALSE)
plot(lres,legend=FALSE,main="closeness centrality \n euclidean distance")
plot(lcshp$geometry,col=JOcols[order(centr$CC,decreasing=FALSE)],main="closeness centrality euclidean distance",add=TRUE)
text(ccoord,labels=labels)
legend("topleft",fill=JOcols,legend=round(sort(centr$CC,decreasing=FALSE),2),bty="n")
centr = centr_igraph(pa=costpm$pa,mdist=costpm$mdist,alpha=40,dispfop="negex",weighted=TRUE,cutoffpr=0.01,MST=FALSE)
plot(lres,legend=FALSE,main="closeness centrality \n cost distance")
plot(lcshp$geometry,col=JOcols[order(centr$CC,decreasing=FALSE)],main="closeness centrality cost distance",add=TRUE)
text(ccoord,labels=labels)
legend("topleft",fill=JOcols,legend=round(sort(centr$CC,decreasing=FALSE),2),bty="n")
#####
## Equivalent connected area, Probability of connectivity, ProtConn (i.e. ECA_AL) functions...
## numerator of probability of connectivity (needs to be divided by the area of the landscape to yield Probability of Connectivity, PC = pcnum/(AL)^2)
# pcnum = PCnum(pa=eucpm$pa,mdist=eucpm$mdist,alpha=40,dispfop="negex",savememory=FALSE)
# ECA_pc(pa=eucpm$pa,mdist=eucpm$mdist,alpha=40,dispfop="negex",savememory=FALSE)
# ECA_pc(pcnum=pcnum)
## wrapper function for all metrics at once:
#PCECA_fun(x=lcr,resi=lres,cost_type="least-cost")
## in order to speed up calculations, directly take distance matrices created above, but then you need to provide area of AL
AL = sf::st_area(st_as_sfc(st_bbox(lcr)))
pceca = PCECA_fun(pa=costpm$pa,mdist=costpm$mdist,AL=AL)
print(pceca)
#> $PCnum
#> [1] 2141225
#>
#> $PC
#> [1] 0.02141225
#>
#> $ECA
#> [1] 1463.292
#>
#> $ECA_AL
#> [1] 0.1463292
#>
#> $ECA_AP
#> [1] 0.9755778
#>
#> $AL
#> [1] 10000
#>
#> $AP
#> [1] 1499.924
######
## patch importance for function dI - there is a dI_fun and a dI_fun_alpha, where dispersal capacities can be indicated!
## this function can be slow if there are many patches (more than 10000...)
## define a connectivity function
confun = "PCnum"
## define an input this function can handle (it can also just be x = lcr)
x = list(pa=costpm$pa,mdist=costpm$mdist)
## importantly, the user needs to specify "k", i.e. how list objects should be subsetted;
expression = "pa=x[[1]][-k],mdist=x[[2]][-k,-k],alpha=alpha"
pid = colnames(costpm$mdist)
resnr = 1
alpha0 = 40
dIPCnum = dI_fun_alpha(fun=confun,expr=expression,x=x,pid=pid,resnr=resnr,alpha=alpha0)
## plot for overview
plot(lres,legend=FALSE,main="Patch importance for PCnum, \n based on cost distance")
plot(lcshp$geometry,col=JOcols[order(dIPCnum,decreasing=FALSE)],add=TRUE)
text(ccoord,labels=labels)
legend("topleft",fill=JOcols,legend=round(sort(dIPCnum,decreasing=FALSE),2),bty="n")
### try it with more difficult function - ECA and PC!
confun = "PCECA_fun"
x = list(pa=costpm$pa,mdist=costpm$mdist)
expression = "pa=x[[1]][-k],mdist=x[[2]][-k,-k],alpha=alpha,AL=10000"
pid = colnames(x[[2]])
resnr = c(2,3)
alpha0 = 40
dIPCECA = dI_fun_alpha(fun=confun,expr=expression,x=x,pid=pid,resnr=resnr,alpha=alpha0)
print(dIPCECA)
#> PC ECA
#> 1 47.61184 27.62033
#> 2 45.96385 26.49072
#> 3 50.49896 29.64303
dIECA = dIPCECA[,2]
## plot for overview
plot(lres,legend=FALSE,main="Patch importance for ECA, \n based on cost distance")
plot(lcshp$geometry,col=JOcols[order(dIECA,decreasing=FALSE)],add=TRUE)
text(ccoord,labels=labels)
legend("topleft",fill=JOcols,legend=round(sort(dIECA,decreasing=FALSE),2),bty="n")
####
## inverse cumulative resistance between x number of pairs in habitat network invCR
print("inverse cumulative resitance \n among all habitat patches")
#> [1] "inverse cumulative resitance \n among all habitat patches"
print(invCR_p(pa=costpm$pa,mdist=costpm$mdist,pnr=1))
#> [1] 2.660293
### patch importance is possible here as well!
confun = "invCR_p"
x = list(pa=costpm$pa,mdist=costpm$mdist)
expression = "pa=x[[1]][-k],mdist=x[[2]][-k,-k],pnr=1"
pid = colnames(x[[2]])
resnr = c(1)
dIinvCR_p = dI_fun(fun=confun,expr=expression,x=x,pid=pid,resnr=resnr)
## this has some negative value...
dIinvCR_p = dIinvCR_p + abs(min(dIinvCR_p)) + 1
## plot for overview
plot(lres,legend=FALSE,main="Patch importance for inverse \n cumulative resistance, based on cost distance")
plot(lcshp$geometry,col=JOcols[order(dIinvCR_p,decreasing=FALSE)],add=TRUE)
text(ccoord,labels=labels)
legend("topleft",fill=JOcols,legend=round(sort(dIinvCR_p,decreasing=FALSE),2),bty="n")
####
## Metapopulation metrics!
mpcres = MPC_fun(pa=costpm$pa,mdist=costpm$mdist,alpha=40)
## plot for overview
plot(lres,legend=FALSE,main="Patch importance for metapop. capacity, \n based on cost distance")
plot(lcshp$geometry,col=JOcols[order(mpcres$pimport*100,decreasing=FALSE)],add=TRUE)
text(ccoord,labels=labels)
legend("topleft",fill=JOcols,legend=round(sort(mpcres$pimport*100,decreasing=FALSE),2),bty="n")
####
## calculate a series of Metapopulation Capacity (MPC) - based metrics
## define region of interest
roi = st_sf(id=1,geometry=st_as_sfc(st_bbox(lcr)))
## define baseline MPC value (e.g. if all the landscape would be habitat)
MPCbasl = MPC_fun(x=roi)$mpc
###
mpcres = MPCser(pa=costpm$pa,mdist=costpm$mdist,alpha=40,roi=roi,MPCbasl=MPCbasl)
## plot for overview
plot(lres,legend=FALSE,main="Patch importance for metapop. capacity, \n based on cost distance")
plot(lcshp$geometry,col=JOcols[order(mpcres$MPCev*100,decreasing=FALSE)],add=TRUE)
text(ccoord,labels=labels)
legend("topleft",fill=JOcols,legend=round(sort(mpcres$MPCev*100,decreasing=FALSE),2),bty="n")
## some further connectivity metrics..
print("some further connectivity metrics...")
#> [1] "some further connectivity metrics..."
## inverse resistance / distance of landscape...
print(invCR_l(x=lcr,resi=lres,cost_type = "least-cost",dist_type="edge"))
#> [1] 0.07044483
## some general Landscape statistics metrics (patch area, number of patches, nearest distances...)
print(LandStats(pa=costpm$pa,mdist=costpm$mdist))
#> tpa mnpa sdpa tnp tnnd mnnnd
#> 1.499924e+03 4.999746e+02 6.785590e-13 3.000000e+00 1.001421e+02 3.338071e+01
#> sdnnd
#> 2.391463e+00
## some general Fragstats metrics (clumpiness index, cohesion index, mesh size, percent of like adjacencies...)
print(FragStats(x=lcr))
#> Warning: Please use 'check_landscape()' to ensure the input data is valid.
#> ai clumpy cohesion mesh nlsi np pd
#> 97.52066 NA 98.05375 0.04970 1.00000 3.00000 2012.07243
#> pladj
#> 94.96982