Function to calculate geodesic distances

Usage

getGeoDist(xygrid, polys, datalocations = NULL, plot.transition = TRUE)

Arguments

xygrid: Regular grid of x and y coordinates, must have an extent greater than that of datalocations if specified.
polys: Either a single polygon (defined by a dataframe with x and y) or a list object containing multiple polygons that define exclusion zones.
datalocations: data frame of x and y coordinates of locations to calculate distances between. Default = NULL, in which case distances are returned between locations in xygrid
plot.transition: (default=TRUE). Logical stating whether to plot the transition matrix. Useful to see if the boundaries are being obeyed.

Details

Using the raster, spancs and gdistance packages to calculate the geodesic distance between pairs of points. The polygons define the areas that an animal/object cannot be. E.g. Land if the animal is a marine mammal.

The object returned is a matrix of distances between pairs of points and the x,y coordinates of each pair of points.

Examples


Nx=100
Ny=100

xygrid<-expand.grid(x=1:Nx, y=1:Ny)

# make exclusion polygon,
bnd<-data.frame(x=c(78  ,79, 52, 53 , 82 , 82, 100, 103 ,104 , 70), 
y= c( 99  ,65 , 66 , 43 , 42 , 16, 19, 16 ,102 ,105))
bnd2<-data.frame(x= c(37 ,54 ,55, 37), y=c(35 ,35, 17, 18))

geodistsoutput<-getGeoDist(xygrid=xygrid, polys=list(bnd, bnd2))
#> Loading required package: raster
#> Loading required package: sp
#> Loading required package: gdistance
#> Warning: package 'gdistance' was built under R version 4.2.3
#> Loading required package: igraph
#> Warning: package 'igraph' was built under R version 4.2.3
#> 
#> Attaching package: 'igraph'
#> The following object is masked from 'package:raster':
#> 
#>     union
#> The following objects are masked from 'package:stats':
#> 
#>     decompose, spectrum
#> The following object is masked from 'package:base':
#> 
#>     union
#> 
#> Attaching package: 'gdistance'
#> The following object is masked from 'package:igraph':
#> 
#>     normalize
#> Error in as.double(y): cannot coerce type 'S4' to vector of type 'double'

# show on plot
i=2150
a.dist.mesh = data.frame(geodistsoutput$xydata, value=geodistsoutput$distance[i,])
#> Error in data.frame(geodistsoutput$xydata, value = geodistsoutput$distance[i,     ]): object 'geodistsoutput' not found
iLCdsitance = rasterFromXYZ(a.dist.mesh)
#> Error in rasterFromXYZ(a.dist.mesh): object 'a.dist.mesh' not found
plot(iLCdsitance, col=topo.colors(100))
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'iLCdsitance' not found
points(geodistsoutput$xydata[i,1],geodistsoutput$xydata[i,2])
#> Error in points(geodistsoutput$xydata[i, 1], geodistsoutput$xydata[i,     2]): object 'geodistsoutput' not found
polymap(bnd, add=T)
#> Error in polymap(bnd, add = T): could not find function "polymap"
polymap(bnd2, add=T)
#> Error in polymap(bnd2, add = T): could not find function "polymap"