ошибка при обрезке sf voronoi полигонов к граничной рамке

Question

Я пытаюсь «обрезать» вороной-полигоны (созданные с помощью 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

Проблемы начинаются, когда я хочу «обрезать» вороной-полигоны v в граничную рамку 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?

Я предпочитаю оставаться в пакетах sf, но, конечно, вся помощь приветствуется!

Answer 1

Возвращение от st_voronoi - это GEOMETRY COLLECTION:

> v
Geometry set for 1 feature 
geometry type:  GEOMETRYCOLLECTION

, и операции с ними могут быть немного странными.Разделение их на деколлированные части обычно помогает:

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