0

我从 Excel 文件上的一系列点生成了 Voronoi 图。我想创建一系列具有相同半径的圆，每个圆都以第二个坐标数据集为中心，并确定每个点的相应重叠区域。

我使用了https://gis.stackexchange.com/questions/190917/r-voronoi-tesselation-from-long-lat-data中定义的函数将 Voronoi 镶嵌创建为 SpatialPolygonsDataFrame。但是，我不确定如何从第二个数据集中生成圆圈。我st_buffer在定义的半径处使用，并尝试使用 转换 Voronoi 结果st_as_sf，但这似乎不起作用。在尝试使用st_intersection来确定区域时，我收到了错误st_crs(x) == st_crs(y)，这似乎来自未定义坐标的 Voronoi 集。

我不确定哪个函数最适合使用，以及是否sf有任何其他工具可用于输出 Voronoi 细分的每个部分与这些半径重叠的区域。

我从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，但我找不到使这些方法与 Voronoi 函数兼容的方法。

我一直在使用的示例测试数据集包含在此处：

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))
4

2 回答 2

1

示例数据：

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,]

创建 voronoi 镶嵌和圆圈

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
于 2018-02-07T18:00:58.190 回答
0

虽然我能够提出这个解决方案，但这可能不是最有效的解决方案。

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))
   }
}
于 2018-02-07T16:40:19.103 回答