Я рассчитал площадь поверхности распределения видов (слияние полигонов из шейп-файлов), но, поскольку эта область может состоять из довольно удаленных полигонов, я хотел бы рассчитать некоторую меру дисперсии. До сих пор я извлекал центроиды каждого многоугольника, вычислял расстояние между ними и использовал их для вычисления коэффициента вариации, как в фиктивном примере ниже;
require(sp)
require(ggplot2)
require(mapdata)
require(gridExtra)
require(scales)
require(rgeos)
require(spatstat)
# Create the coordinates for 3 squares
ls.coords <- list()
ls.coords <- list()
ls.coords[[1]] <- c(15.7, 42.3, # a list of coordinates
16.7, 42.3,
16.7, 41.6,
15.7, 41.6,
15.7, 42.3)
ls.coords[[2]] <- ls.coords[[1]]+0.5 # use simple offset
ls.coords[[3]] <- c(13.8, 45.4, # a list of coordinates
15.6, 45.4,
15.6, 43.7,
13.8, 43.7,
13.8, 45.4)
# Prepare lists to receive the sp objects and data frames
ls.polys <- list()
ls.sp.polys <- list()
for (ii in seq_along(ls.coords)) {
crs.args <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
my.rows <- length(ls.coords[[ii]])/2
# create matrix of pairs
my.coords <- matrix(ls.coords[[ii]],nrow = my.rows,ncol = 2,byrow = TRUE)
# now build sp objects from scratch...
poly = Polygon(my.coords)
# layer by layer...
polys = Polygons(list(poly),1)
spolys = SpatialPolygons(list(polys))
# projection is important
proj4string(spolys) <- crs.args
# Now save sp objects for later use
ls.sp.polys[[ii]] <- spolys
# Then create data frames for ggplot()
poly.df <- fortify(spolys)
poly.df$id <- ii
ls.polys[[ii]] <- poly.df
}
# Convert the list of polygons to a list of owins
w <- lapply(ls.sp.polys, as.owin)
# Calculate the centroids and get the output to a matrix
centroid <- lapply(w, centroid.owin)
centroid <- lapply(centroid, rbind)
centroid <- lapply(centroid, function(x) rbind(unlist(x)))
centroid <- do.call('rbind', centroid)
# Create a new df and use fortify for ggplot
centroid_df <- fortify(as.data.frame(centroid))
# Add a group column
centroid_df$V3 <- rownames(centroid_df)
ggplot(data = italy, aes(x = long, y = lat, group = group)) +
geom_polygon(fill = "grey50") +
# Constrain the scale to 'zoom in'
coord_cartesian(xlim = c(13, 19), ylim = c(41, 46)) +
geom_polygon(data = ls.polys[[1]], aes(x = long, y = lat, group = group), fill = alpha("red", 0.3)) +
geom_polygon(data = ls.polys[[2]], aes(x = long, y = lat, group = group), fill = alpha("green", 0.3)) +
geom_polygon(data = ls.polys[[3]], aes(x = long, y = lat, group = group), fill = alpha("lightblue", 0.8)) +
coord_equal() +
# Plot the centroids
geom_point(data=centroid_points, aes(x = V1, y = V2, group = V3))
# Calculate the centroid distances using spDists {sp}
centroid_dists <- spDists(x=centroid, y=centroid, longlat=TRUE)
centroid_dists
[,1] [,2] [,3]
[1,] 0.00000 69.16756 313.2383
[2,] 69.16756 0.00000 283.7120
[3,] 313.23834 283.71202 0.0000
# Calculate the coefficient of variation as a measure of polygon dispersion
cv <- sd(centroid_dist)/mean(centroid_dist)
[1] 0.9835782
Участок трех полигонов и их центроидов
Я не уверен, что этот подход очень полезен, поскольку во многих случаях некоторые из полигонов (как синий в приведенном выше примере) довольно велики по сравнению с остальными, что увеличивает расстояние еще больше. Например, центр тяжести Австралии имеет почти такое же расстояние до своих западных границ, как и до Папау.
То, что я хотел бы получить, является некоторой информацией об альтернативных подходах. Например, как или с помощью какой функции можно рассчитать расстояние между полигонами?
Я протестировал преобразование указанного выше фрейма данных SpatialPolygon в PointPatterns (ppp), {spatstat}
чтобы можно было бежать nndist() {spatstat}
для расчета расстояния между всеми точками. Но так как я имею дело с довольно большими областями (многоугольниками и большими), матрица становится огромной, и я не уверен, как продолжать добираться до минимального расстояния между многоугольниками .
Я также посмотрел на функцию gDistance {rgeos}
, но я думаю, что она работает только с проецируемыми данными, что может быть проблемой для меня, так как мои области могут пересекать несколько EPSG areas
. Та же проблема возникнет для функции crossdist {spatstat}
.
postgres/postgis
в дополнение кR
? Я использовал рабочий процесс, в котором я выполняю большую часть своей работыR
, но храню данные в базе данных, к которой я обращаюсьsqldf
. Это позволяет вам использовать всеpostgis
функции (из которых расстояние между полигонами прямое)postgres
создавать базу данных, но остановился, когда не знал (не смотрел), как связать рабочий процесс / геостаты между базой данных иR
...Ответы:
Вы можете сделать этот анализ в пакете "spdep". В соответствующих соседних функциях, если вы используете "longlat = TRUE", функция вычисляет расстояние по большому кругу и возвращает километры в качестве единицы расстояния. В приведенном ниже примере вы можете привести получившийся объект списка расстояний («dist.list») к матрице или data.frame, однако, это довольно эффективный способ расчета сводной статистики с использованием lapply.
источник
spdeb
пакета. Просто для пояснения, этот подход дает тот же результат, что и в моем примере, верно?