将 sf voronoi 多边形裁剪到边界框时出错 [英] error when clipping sf voronoi polygons to boundary box
本文介绍了将 sf voronoi 多边形裁剪到边界框时出错的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我正在尝试将 voronoi-polygons(使用 sf-package 创建)剪辑"到边界框...但它向我抛出了一个我无法定义的错误.我对 R 的空间世界不是很有经验.
I am trying to 'clip' voronoi-polygons (created with the sf-package) to a boundary box... But it throws me an error that I cannot define. I'm not very experienced in the spatial-worldof 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 ) )
在创建 voronoi 多边形时，我收到警告
On creating the voronoi-polygons, I get a warning
>Warning message:
>In st_voronoi.sfc(st_union(stations)) :
>st_voronoi does not correctly triangulate longitude/latitude data
不确定此警告对我的问题是否重要.到目前为止看起来还不错.. 绘图符合预期
Not sure if this warning is of any importance to my problem. So far so good it seems.. plotting is as expected
mapview( list( v, box ) )
结果:
其中较小的矩形是边界框
where the smaller rectangle is the boundary box box
当我想将 voronoi 多边形 v '剪切'到边界框
Problems begin when I want to 'clip' the voronoi-polygons v to the boundary box box...
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 中的多边形裁剪到 box 中的多边形/框吗?
Anyone got an idea on how to clip the polygons in v to the polygon/box in box?
我更喜欢留在 sf 包中，但当然感谢所有帮助！
I prefer to stay within the sf packages, but of course all help is appreciated!
推荐答案
st_voronoi 返回的是一个 GEOMETRY COLLECTION:
> v
Geometry set for 1 feature
geometry type: GEOMETRYCOLLECTION
和这些操作可能有点奇怪.将它们拆分成未收集的部分通常会有所帮助:
and operations with those can be a bit weird. Splitting them into their decollected parts usually helps:
> vparts = st_collection_extract(v)
> vcrop = st_crop(vparts, box)
although coordinates are longitude/latitude, st_intersection assumes that they are planar
> plot(vcrop,col=NA)
这篇关于将 sf voronoi 多边形裁剪到边界框时出错的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
