Circuitscape is a tool that is used to measure effective resistance between pairs of nodes. A more recent package in R called "gdistance" can also perform a similar measurement and offers other function to measure effective distance such as cost distance. In this script I compare all 3 and include euclidean distance.

CompareEffectiveDistances.R

# Compare CS, costDist and commuteDist

# Clear workspace
rm(list=ls())

# Load libraries
library(gdistance)
library(scales)
library(gstat)

# Setup CS
CS_exe <- 'C:/"Program Files"/Circuitscape/cs_run.exe' # Don't forget the "Program Files" problem

# Make a template CS .ini file
CS_ini <- c("[circuitscape options]",            
            "data_type = raster",
            "scenario = pairwise",
            paste(c("point_file =",
                    "habitat_file =",
                    "output_file ="),
                  paste(getwd(),c("CS_output/sites_rast.asc",
                                  "CS_output/landscape.asc",
                                  "CS_output/CS.out"),
                        sep="/")))

# Write the CS .ini it to your working directory
writeLines(CS_ini,"myini.ini")
rm(CS_ini)

# Make the CS run cmd
CS_run <- paste(CS_exe, paste(getwd(),"myini.ini",sep="/")) # Make the cmd
rm(CS_exe)

# Function that run different movement algorithms on a simulated landscape. 
# There is the option to simulate 1 pair of points or several pair of points
# It plots the landscape and points and runs the 3 different distance algorithms. 
# I included euclidean distance as a forth.
calculate_distances <- function(seed_num,dimension=100,multiples_pnts=F){
  # set seed
  set.seed(seed_num)
  
  # Print
  cat("\n")
  cat("Simulation # ",seed_num, ".................................\n")
  cat("\n")
  
  if(multiples_pnts){
    # make n sites = dimension
    sites <- sample(x = (dimension/10):(dimension-dimension/10),replace = F, size = 2*30)
  }else{
    # make 1 pair of points
    sites <- sample(x = (dimension/10):(dimension-dimension/10),replace = T, size = 4)
  }
  
  # Create sp points
  sites <- matrix(sites,ncol = 2)
  sites <- SpatialPoints(sites)
  
  # Simulate landscape
  xy <- expand.grid(1:dimension, 1:dimension)
  names(xy) <- c("x","y")
  rm(dimension)
  
  # define the gstat object (spatial model)
  g.dummy <- gstat(formula=z~x+y, locations=~x+y, dummy=T, beta=1, 
                   model=vgm(psill=0.002,model="Exp",range=2), nmax=20)
  
  # make a simulation based on the gstat object
  yy <- predict(g.dummy, newdata=xy, nsim=1)
  
  # show one realisation
  gridded(yy) = ~x+y
  landscape <- raster(yy)
  landscape[] <- round(rescale(landscape[],to = c(1,500)))
  
  # Remove junk
  rm(g.dummy,yy,xy)
  
  # Plot it
  plot(landscape)
  points(sites,pch=19,col=2)
  
  # Rasterize points using the cost extent
  sites_rast <- rasterize(x = sites,y = landscape) # So that it has same extent as landscape
  
  # Create CS directory
  dir.create("CS_output")
  
  # Write rasters to your working directory
  writeRaster(sites_rast,"CS_output/sites_rast.asc",overwrite=TRUE)
  writeRaster(landscape,"CS_output/landscape.asc",overwrite=TRUE)
  
  # Euclidean distance
  geo_dist <-as.dist(pointDistance(sites, lonlat = F))
  
  # Run the CS command
  system(CS_run)
  
  # Import the effective resistance
  cs_dist <- as.dist(read.csv("CS_output/CS_resistances.out",sep=" ",row.names=1,header=1))
  
  # Run gdistance
  trCost <- transition(x=1/landscape,transitionFunction=mean,directions=8) #Transition object from the raster
  cost_dist <- costDistance(x=trCost,fromCoords=coordinates(sites))
  
  # Get commute distance (~ resistance distance)
  comm_dist <- commuteDistance(x=trCost,coords=coordinates(sites))
  
  # Combine results
  results <- cbind(geo_dist=geo_dist,cost_dist=cost_dist,cs_dist=cs_dist,comm_dist=comm_dist)
  
  # Remove stuff
  rm(landscape,sites_rast,cs_dist,trCost,cost_dist,comm_dist,sites,geo_dist)
  
  # remove CS directory and everything in it
  unlink("CS_output", recursive = T, force = T)
  
  # Return results
  return(results)
}

### Part 1: Run 1 simulation with 100 pairs of points ###
part1_out <- calculate_distances(seed_num=1, multiples_pnts=T)
part1_out <- apply(part1_out,MARGIN = 2,scale) # Standardize output
colnames(part1_out)
pairs(part1_out,labels=c("Euclidean distance","Least-cost distance","Resistance distance","Commute distance"))
cor(part1_out)

### Part 2: Run 1000 landscape simulations with 1 pair of points each ####
part2_out <- lapply(X = 1:1000,FUN = calculate_distances)
file.remove("myini.ini") # Remove ini template
part2_out <- do.call("rbind",part2_out) # rbind the output to make a nice table
part2_out <- apply(part2_out,MARGIN = 2,scale) # Standardize output
colnames(part2_out)
pairs(part2_out,labels=c("Euclidean distance","Least-cost distance","Resistance distance","Commute distance"))
cor(part2_out)