require(tidyr)
require(dplyr)
require(raster)
require(geoRcb)
require(sp)
require(spatstat)
require(ggplot2)
require(viridis)
require(cbK)
theme_set(theme_bw())

Exploration

## The locations are LazyLoaded with the package
## If you had yourself a shapefile, for example, you 
## would do something like
# basename <- "indigenous_settlements"
# pts <- rgdal::readOGR(dsn = path, layer = basename)

## The cost surface is provided with the package as a
## separate raster file
dn <- system.file("extdata", package="cbK")
cambaceres_cs <- raster(file.path(dn, "cambaceres_cs"))

There seem to be some pixels with outlying cost values (an order of magnitude higher). These have been filtered in the surface plot.

boxplot(getValues(cambaceres_cs), horizontal = TRUE)
data.frame(
  coordinates(cambaceres_cs),
  cost = getValues(cambaceres_cs)
) %>% 
  ggplot(aes(x, y)) +
  geom_raster(aes(fill = cost)) +
  geom_point(
    data = as.data.frame(indigenous_settlements),
    aes(coords.x1, coords.x2)) +
  scale_fill_viridis() +
  coord_fixed()

Cost-based clustering calculations

res <- cbK(
  loc = indigenous_settlements,
  cs  = cambaceres_cs,
  directions = 16
)

Multi-Dimensional Scaling

In this case, the 2D MDS representation captures the r round(res$mds$GOF[1]*100, 1) % of the variability in the cost-based distances.

Ripley's $K$ function

plot(res$K)
plot(res$env_K)

Or, its stabilized transformation: the $L$ function

plot(res$L)
plot(res$env_L)

Comparison with Euclidean distances

The cost-based distances are an order of magnitude higher than the Euclidean distances.

data.frame(
  eucl = unclass(res$eucl$dm),
  cb = res$dm[lower.tri(res$dm)]
) %>% 
  ggplot(aes(eucl, cb)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, col = 'darkgray')

For reference, we have also computed $K$ and $L$ functions with euclidean coordinates.

plot(res$eucl$K)
plot(res$eucl$env_K)

plot(res$eucl$L)
plot(res$eucl$env_L)

Conclusions from Cabaceres case study



famuvie/cbK documentation built on May 16, 2019, 10:04 a.m.