Tengo una ruta como una serie ordenada de puntos que se cruza con un polígono. Me gustaría extraer la parte de la línea antes de la primera intersección del polígono.

Intenté dividir la ruta calculando la diferencia con el polígono, pero la línea también se divide en sus intersecciones (ver ejemplo). Necesito extraer la parte completa de la ruta desde el principio hasta que se cruce por primera vez con el polígono (cuadrado azul en el ejemplo).

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

En el ejemplo, la ruta está segmentada por todas las intersecciones (incluso en sí misma), por lo que la primera parte de la ruta no se extiende hasta el polígono. Sospecho que hay una función en el paquete sf específicamente para mi propósito, pero aún no la he encontrado.

r sf
1
randr 23 oct. 2019 a las 18:33

1 respuesta

La mejor respuesta

Una cadena de líneas auto-intersectadas no es "simple" (vea st_is_simple(l)) aunque sea "válida" (vea st_is_valid(l)), y sf aquí le brinda una característica simple. Lo que tenemos que hacer es trabajar con LINESTRINGs simples y reconstruir la característica de línea más larga ...

Primero diferencie la línea ondulada con el polígono - d es entonces una multilínea:

> d = st_difference(l,p)

Si lo convertimos a sfc y luego lo convertimos en LINESTRING obtenemos las partes en geometrías de LINESTRING separadas:

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

Y hay 8 de ellos:

> length(dl)
[1] 8

Entonces, comenzando con el primero, ¿cuántos tenemos que seguir hasta que lleguemos al polígono? Obtengamos las intersecciones de los segmentos con el polígono; esta es una lista dispersa; los elementos tienen una longitud de 0 si no hay intersección, mayor que 0 si hay una intersección:

> sint = st_intersects(dl,p)

El primer elemento de esa lista que no es cero es la primera cadena de línea que llega al polígono:

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

Lo que significa que nuestra línea hasta el polígono son las primeras tres cadenas de líneas:

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

Puede unirlos en una sola característica y trazar:

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

enter image description here

0
Spacedman 26 oct. 2019 a las 17:29