MCP95_from_polygons.R

# Clear workspace
rm(list=ls())

# load libraries
library(rgeos)
library(maptools)

# Set seed
set.seed(1234)

# The number of spatial units in our example
nber = 100

# Range of the coordinates
range_coords = c(1,1000)

# Mean Size of the spatial units
mean_size = 0.001*(range_coords[2] - range_coords[1])^2 # 10% of the area
radius = sqrt(mean_size/pi)

# Standard deviation of the radius of this size
sd_size <- 1.5*radius # want it to vary alot

# Make some random polygons spatial units, in my case there were traplines and township polygons
centroids <- SpatialPointsDataFrame(coords = cbind("x" = sample(x = range_coords[1]:range_coords[2], size = nber),
                                                   "y" = sample(x = range_coords[1]:range_coords[2], size = nber)), 
                                    data = data.frame("radius" = rnorm(n = nber, 
                                                                       mean = radius, 
                                                                       sd = sd_size)))

polygons <- gBuffer(spgeom = centroids, byid = T, width = centroids$radius)

# Plot, these are my spatial units
plot(polygons)
points(centroids, pch = 19, col = 2)

# Now we want to calculate some distance between these polygons 
# to be able to make a 95% MCP of these polygons

# Edge distance
poly_dist <- gDistance(polygons, byid = T)

# Now we sum these distance by polygon ID
sum_poly_dist <- colSums(poly_dist)

# Now lets find the 95% distance cutoff for both distance using quantiles
qcut95_poly <- quantile(sum_poly_dist, 0.95)

# Now we select polygons for each distance type using this 95% cutoff
polys95 <- polygons[which(sum_poly_dist < qcut95_poly),]

# How many polygons were selected, there were 100 polygons, we should have lost 5
length(polys95)

# Now lets make the 95% mcp
mcp95_bypolys <- gConvexHull(polys95)

# Plot them
plot(polygons, col = 2)
plot(polys95, col = 3, add = T)
plot(mcp95_bypolys, border = "blue", add = T)

# Lets compare this to a mcp created from just the centroid points
centroid_dist <- as.matrix(dist(coordinates(centroids)))
sum_centroid_dist <- colSums(centroid_dist)
qcut95_centroid <- quantile(sum_centroid_dist, 0.95)
centroids95 <- centroids[which(sum_centroid_dist < qcut95_centroid),]
mcp95_bypnts <- gConvexHull(centroids95)

# Plot them
plot(centroids, pch = 19, col = 2)
plot(centroids95, pch = 19, col = 3, add = T)
plot(mcp95_bypnts, border = "blue", add = T)

# Compare both 95% MCPs
plot(polygons, col = 2)
points(centroids, col = 1, pch = 19)
plot(mcp95_bypolys, border = "blue", add = T)
plot(mcp95_bypnts, border = "green", add = T)