Как найти путь вдоль линии, пока она сначала не пересекает многоугольник?

Question

У меня есть путь в виде упорядоченного ряда точек, который пересекает многоугольник. Я хотел бы извлечь часть линии перед первым пересечением многоугольника.

Я попытался разделить путь, рассчитав разницу с многоугольником, но линия также разбивается на свои самопересечения (см. пример). Мне нужно извлечь полную часть пути от начала до тех пор, пока он сначала не пересечет многоугольник (синий квадрат в примере).

# A wonky line that intersects itself   
l = sf::st_linestring(cbind(cos((0:100) * pi / 50), sin((0:100) * pi / 15 )))
# A polygon that intersects the line
p = sf::st_polygon(list(cbind(c(-.3, -.3, -.2, -.2, -.3), c(-.3, -.2, -.2, -.3, -.3))))

# Visualisation of the problem
plot(l)
plot(p, add = TRUE, col = "blue")

# Taking the first fragment of the difference does not work as the path intersects and divides itself
d = sf::st_difference(l, p)
plot(sf::st_linestring(d[[1]]), add = TRUE, col = "red")

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

Answer 1

Самопересекающаяся строка строки не является «простой» (см. st_is_simple(l)), даже если она «действительна» (см. st_is_valid(l)), и sf здесь предоставляет вам простую функцию. Что нам нужно сделать, так это работать с простыми LINESTRING и перестроить объект более длинной линии ...

Сначала отразим волнистую линию с многоугольником - d - это многоканальная строка:

> d = st_difference(l,p)

Если мы преобразуем его в sfc и затем приведем его к LINESTRING, мы получим детали в отдельной геометрии LINESTRING:

> dl = st_cast(st_sfc(d),"LINESTRING")

, и их будет 8:

> length(dl)
[1] 8

Итак, начиная с первого, сколько нам нужно следовать, пока мы не достигнем многоугольника? Давайте получим пересечения сегментов с многоугольником - это разреженный список - элементы имеют длину 0, если пересечения нет, больше 0, если есть пересечение:

> sint = st_intersects(dl,p)

Первый элемент этого списка , который не равен нулю - первая строка строки, которая попадает в многоугольник:

> min(which(lengths(sint)>0))
[1] 3

Это означает, что наша линия до многоугольника - это первые три строки строки:

> dl[1:3]
Geometry set for 3 features 
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: -0.220294 ymin: -0.9945219 xmax: 1 ymax: 0.9945219
epsg (SRID):    NA
proj4string:    NA
LINESTRING (1 0, 0.9980267 0.2079117, 0.9921147...
LINESTRING (0.9510565 0.8660254, 0.9297765 0.95...
LINESTRING (0.309017 -0.8660254, 0.309017 -0.86...

Вы можете объединить их в одну особенность и сюжет:

> plot(l)
> plot(p,add=TRUE,col="red")
> topoly = st_union(dl[1:3])
> plot(topoly, add=TRUE, lwd=3)