require(tidyr) require(dplyr) require(raster) require(geoRcb) require(sp) require(spatstat) require(ggplot2) require(viridis) require(cbK) theme_set(theme_bw())
## 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()
res <- cbK( loc = indigenous_settlements, cs = cambaceres_cs, directions = 16 )
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.
plot(res$K) plot(res$env_K)
Or, its stabilized transformation: the $L$ function
plot(res$L) plot(res$env_L)
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.