Области перекрывающихся кругов на диаграмме Вороного

Я создал диаграмму Вороного из ряда точек в файле Excel. Я хотел бы создать серию кругов одного и того же радиуса, каждый из которых сосредоточен вокруг второго набора данных координат, и определить соответствующие перекрывающиеся области для каждой точки.

Я использовал функцию, определенную в https://gis.stackexchange.com/questions/190917/r-voronoi-tesselation-from-long-lat-data для создания тесселяции Вороного в виде SpatialPolygonsDataFrame. Однако я не уверен, как генерировать круги из второго набора данных. Я использовал st_buffer для определенного радиуса и попытался преобразовать результаты Вороного, используя st_as_sf, но это, похоже, не сработало. При попытке использовать st_intersection для определения областей я получил ошибку st_crs(x) == st_crs(y), которая, по-видимому, исходит из набора Вороного, не имеющего определенных координат.

Я не уверен, какую функцию лучше всего использовать и есть ли у sf какие-либо дополнительные инструменты, которые можно использовать для вывода области, в которой каждый сегмент тесселяции Вороного перекрывается этими радиусами.

Я нашел похожие проблемы из https://gis.stackexchange.com/questions/229453/create-a-circle-of-defined-radius-around-a-point-and-then-find-the-overlapping-a?noredirect=1&lq=1 и https://gis.stackexchange.com/questions/140504/extracting-intersection-areas-in-r, но я не могу найти способ сделать эти методы совместимыми с функцией Вороного.

Пример набора тестовых данных, который я использовал, приведен здесь:

Hubs <- cbind('X Coord' = c(52.37639999,52.36989975,52.86299896,52.01011658,51.67409897,50.84980011,51.88669968,52.1048943,52.0746994), 
          'Y Coord' = c(4.894589901,4.876679897,6.493730068,4.703330517,4.548630238,5.687580109,4.491449833,5.0528496,4.310130119),
          'Info 1'  = c(13,15,62,24,9,46,73,97,69))

References <- cbind('X Coord' = c(51.88529968,52.3360136,52.37440109,51.92269897,51.9192276,51.43019867,52.0780896,51.90299988,52.04100037,50.98810196),
                'Y Coord' = c(4.493549824,4.8750552,4.895800114,4.47382021,4.481627464,5.451980114,5.1470653,4.458670139,4.318089962,5.78059721))

r voronoi areas
person AlwaysInTheDark    schedule 29.01.2018    source источник

Ответы (2)


arrow_upward
1
arrow_downward

Пример данных:

library(dismo)
p <- cbind(x=c(17, 42, 85, 70, 19, 53, 26, 84, 84, 46, 48, 85, 4, 95, 48, 54, 66, 74, 50, 48), 
      y=c(28, 73, 38, 56, 43, 29, 63, 22, 46, 45, 7, 60, 46, 34, 14, 51, 70, 31, 39, 26))

p1 <- p[11:20,]
p2 <- p[1:10,]

Создание мозаики Вороного и кругов

vr <- voronoi(p1)
names(vr) <- "vID"
cr <- circles(p2, 5, lonlat=FALSE, dissolve=FALSE)
cr <- SpatialPolygonsDataFrame(polygons(cr), data.frame(cID=1:10))

Пересечение и вычисление площади

x <- intersect(vr, cr)
crs(x) <- "+proj=utm +zone=1 +datum=WGS84" 
x$area <- area(x)
data.frame(x)

Это использовало Spatial объектов. Теперь с sf после AlwaysInTheDark, но без цикла

library(sf)
sf2 <- st_as_sf(data.frame(p2), coords = c("x", "y"), crs = 3035)
circles <- st_buffer(sf2, dist = 5)
circles$cID <- 1:nrow(circles)

sf1 <- st_as_sf(data.frame(p1), coords = c("x", "y"), crs = 3035)
voroi <- st_collection_extract(st_voronoi(st_combine(sf1)))

d <- st_intersection(circles, voroi)
d$area <- st_area(d)

d
person Robert Hijmans    schedule 07.02.2018

arrow_upward
0
arrow_downward

Хотя я смог придумать это решение, возможно, это не самое эффективное решение.

library(raster)
library(tidyverse)
library(sf)
library(sp)
library(rgdal)
library(rgeos)
library(maptools)
#> Linking to GEOS 3.5.0, GDAL 2.1.1, proj.4 4.9.3

voronoipolygons = function(layer) {
require(deldir)
crds = layer@coords
z = deldir(crds[,1], crds[,2])
w = tile.list(z)
polys = vector(mode='list', length=length(w))
require(sp)
for (i in seq(along=polys)) {
    pcrds = cbind(w[[i]]$x, w[[i]]$y)
    pcrds = rbind(pcrds, pcrds[1,])
    polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
}
SP = SpatialPolygons(polys)
voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(dummy = 
seq(length(SP)), row.names=sapply(slot(SP, 'polygons'), 

function(x) slot(x, 'ID'))))
}

Hubs = read.csv("Test_Hub_Data.csv", header=TRUE)
coordinates(Hubs) <- c("X.Coord", "Y.Coord")

vp <- voronoipolygons(Hubs)

#Voronoi Setup
crs(vp) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
vp_test <- sf::st_as_sf(vp) %>% st_transform(3035)

References = read.csv("Test_Reference_Points.csv", header = TRUE)

# Convert to sf, set the crs to EPSG:4326 (lat/long), 
# and transform to EPSG:3035
Ref_sf <- st_as_sf(References, coords = c("X.Coord", "Y.Coord"), crs = 4326) %>% st_transform(3035)

#Initialise matrix of areas x: Voronoi Segments, y: Nodes
areas <- matrix(0,length(References[,1]),length(Hubs[,1]))

for (i in 1:length(References[,1])){
   #New circles need to be generated for each point
   # Buffer circles by 100m
   Ref_circle <- st_buffer(Ref_sf[i,1], dist = 100)

   for(j in 1:length(Hubs[,1])){
   # Intersect the circles with the polygons
        vp_int_circle <- st_intersection(vp_test[j,1], Ref_circle)

        #Append array
        areas[i,j] = sum(st_area(vp_int_circle))
   }
}
person AlwaysInTheDark    schedule 07.02.2018