我试图'剪辑'沃罗诺多边形(创建与sf包)到一个边界框。。。但它给了我一个我无法定义的错误。我对R的空间世界不是很有经验。
感谢所有的帮助。
样本数据
stations <- structure(list(STN = c(209L, 210L, 215L, 225L, 235L, 240L, 242L,
248L, 249L, 251L, 257L, 258L, 260L, 267L, 269L, 270L, 273L, 275L,
277L, 278L, 279L, 280L, 283L, 285L, 286L, 290L, 308L, 310L, 311L,
312L, 313L, 315L, 316L, 319L, 323L, 324L, 330L, 331L, 340L, 343L,
344L, 348L, 350L, 356L, 370L, 375L, 377L, 380L, 391L), geometry = structure(list(
structure(c(4.518, 52.465), class = c("XY", "POINT", "sfg"
)), structure(c(4.43, 52.171), class = c("XY", "POINT", "sfg"
)), structure(c(4.437, 52.141), class = c("XY", "POINT",
"sfg")), structure(c(4.555, 52.463), class = c("XY", "POINT",
"sfg")), structure(c(4.781, 52.928), class = c("XY", "POINT",
"sfg")), structure(c(4.79, 52.318), class = c("XY", "POINT",
"sfg")), structure(c(4.921, 53.241), class = c("XY", "POINT",
"sfg")), structure(c(5.174, 52.634), class = c("XY", "POINT",
"sfg")), structure(c(4.979, 52.644), class = c("XY", "POINT",
"sfg")), structure(c(5.346, 53.392), class = c("XY", "POINT",
"sfg")), structure(c(4.603, 52.506), class = c("XY", "POINT",
"sfg")), structure(c(5.401, 52.649), class = c("XY", "POINT",
"sfg")), structure(c(5.18, 52.1), class = c("XY", "POINT",
"sfg")), structure(c(5.384, 52.898), class = c("XY", "POINT",
"sfg")), structure(c(5.52, 52.458), class = c("XY", "POINT",
"sfg")), structure(c(5.752, 53.224), class = c("XY", "POINT",
"sfg")), structure(c(5.888, 52.703), class = c("XY", "POINT",
"sfg")), structure(c(5.873, 52.056), class = c("XY", "POINT",
"sfg")), structure(c(6.2, 53.413), class = c("XY", "POINT",
"sfg")), structure(c(6.259, 52.435), class = c("XY", "POINT",
"sfg")), structure(c(6.574, 52.75), class = c("XY", "POINT",
"sfg")), structure(c(6.585, 53.125), class = c("XY", "POINT",
"sfg")), structure(c(6.657, 52.069), class = c("XY", "POINT",
"sfg")), structure(c(6.399, 53.575), class = c("XY", "POINT",
"sfg")), structure(c(7.15, 53.196), class = c("XY", "POINT",
"sfg")), structure(c(6.891, 52.274), class = c("XY", "POINT",
"sfg")), structure(c(3.379, 51.381), class = c("XY", "POINT",
"sfg")), structure(c(3.596, 51.442), class = c("XY", "POINT",
"sfg")), structure(c(3.672, 51.379), class = c("XY", "POINT",
"sfg")), structure(c(3.622, 51.768), class = c("XY", "POINT",
"sfg")), structure(c(3.242, 51.505), class = c("XY", "POINT",
"sfg")), structure(c(3.998, 51.447), class = c("XY", "POINT",
"sfg")), structure(c(3.694, 51.657), class = c("XY", "POINT",
"sfg")), structure(c(3.861, 51.226), class = c("XY", "POINT",
"sfg")), structure(c(3.884, 51.527), class = c("XY", "POINT",
"sfg")), structure(c(4.006, 51.596), class = c("XY", "POINT",
"sfg")), structure(c(4.122, 51.992), class = c("XY", "POINT",
"sfg")), structure(c(4.193, 51.48), class = c("XY", "POINT",
"sfg")), structure(c(4.342, 51.449), class = c("XY", "POINT",
"sfg")), structure(c(4.313, 51.893), class = c("XY", "POINT",
"sfg")), structure(c(4.447, 51.962), class = c("XY", "POINT",
"sfg")), structure(c(4.926, 51.97), class = c("XY", "POINT",
"sfg")), structure(c(4.936, 51.566), class = c("XY", "POINT",
"sfg")), structure(c(5.146, 51.859), class = c("XY", "POINT",
"sfg")), structure(c(5.377, 51.451), class = c("XY", "POINT",
"sfg")), structure(c(5.707, 51.659), class = c("XY", "POINT",
"sfg")), structure(c(5.763, 51.198), class = c("XY", "POINT",
"sfg")), structure(c(5.762, 50.906), class = c("XY", "POINT",
"sfg")), structure(c(6.197, 51.498), class = c("XY", "POINT",
"sfg"))), class = c("sfc_POINT", "sfc"), precision = 0, bbox = structure(c(xmin = 3.242,
ymin = 50.906, xmax = 7.15, ymax = 53.575), class = "bbox"), crs = structure(list(
epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), n_empty = 0L)), row.names = c(NA,
-49L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(STN = NA_integer_), class = "factor", .Label = c("constant",
"aggregate", "identity")))
目前为止的代码
library( sf )
#Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
library( mapview )
#create boundary box
box <- st_bbox( stations ) %>% st_as_sfc()
#craete voronoi polygons
v <- st_voronoi(st_union( stations ) )
>Warning message:
>In st_voronoi.sfc(st_union(stations)) :
>st_voronoi does not correctly triangulate longitude/latitude data
我不知道这个警告对我的问题有没有重要意义。
到目前为止似乎还不错。。绘图如预期
mapview( list( v, box ) )
结果:
其中较小的矩形是边界框
box
当我想“剪辑”voronoi多边形时,问题就开始了
v
到边界盒
箱
...
st_crop( v, box )
> although coordinates are longitude/latitude, st_intersection assumes that they are planar
>Error in CPL_geos_op2(op, st_geometry(x), st_geometry(y)) :
>Evaluation error: TopologyException: no outgoing dirEdge found at 3.242 51.367318548387097.
问题
有人知道怎么把多边形夹进去吗
v
箱
?
sf
包,但当然所有的帮助是感激的!
