¿Encontrar segmentos de línea más cercanos para apuntar usando bien proporcionado?

17

Antecedentes

Desde un punto conocido, requiero establecer el "perímetro visible" circundante más cercano contra una tabla de MultiLineStrings, como se muestra en el diagrama.

He buscado en este sitio con una serie de términos (por ejemplo, borde mínimo, perímetro mínimo, vecino más cercano, clip, que contiene polígono, visibilidad, ajuste, nodos de corte, trazado de rayos, relleno de inundación, límite interno, enrutamiento, casco cóncavo) pero No puedo encontrar ninguna pregunta previa que parezca coincidir con este escenario.

Diagrama

  • El círculo verde es el punto conocido.
  • Las líneas negras son las conocidas MultiLineStrings.
  • Las líneas grises son una indicación de un barrido radial desde el punto conocido.
  • Los puntos rojos son la intersección más cercana del barrido radial y las MultiLineStrings.

ingrese la descripción de la imagen aquí

Parámetros

  • El punto nunca se cruzará con MultiLineStrings.
  • El punto siempre estará centrado nominalmente dentro de MultiLineStrings.
  • MultiLineStrings nunca encerrará completamente el Punto, por lo tanto, el perímetro será MultiLineString.
  • Habrá una tabla que contiene aproximadamente 1,000 MultiLineStrings (que normalmente contiene una sola línea de aproximadamente 100 puntos).

Metodología considerada

  • Realice un barrido radial construyendo una serie de líneas desde el Punto conocido (en, digamos, incrementos de 1 grado).
  • Establezca el punto de intersección más cercano de cada línea de barrido radial con MultiLineStrings.
  • Cuando una de las líneas de barrido radial no se cruza con ninguna de las MultiLineString, esto indicaría un espacio en el perímetro que se acomodaría en la construcción del perímetro MultiLineString.

Resumen

Si bien esta técnica encontrará las intersecciones más cercanas, no necesariamente encontrará todos los puntos del nodo perimetral más cercanos, dependiendo de la resolución del barrido radial. ¿Alguien puede recomendar un método alternativo para establecer todos los puntos del perímetro o complementar la técnica de barrido radial con alguna forma de amortiguación, sectorización o compensación?

Software

Prefiero usar SpatiaLite y / o Shapely para la solución, pero agradecería cualquier sugerencia que pueda implementarse utilizando un software de código abierto.

Editar: Solución de trabajo (basada en la respuesta de @gene)

from shapely.geometry import Point, LineString, mapping, shape
from shapely.ops import cascaded_union
from shapely import affinity
import fiona

sweep_res = 10  # sweep resolution (degrees)
focal_pt = Point(0, 0)  # radial sweep centre point
sweep_radius = 100.0  # sweep radius

# create the radial sweep lines
line = LineString([(focal_pt.x,focal_pt.y), \
                   (focal_pt.x, focal_pt.y + sweep_radius)])

sweep_lines = [affinity.rotate(line, i, (focal_pt.x, focal_pt.y)) \
               for i in range(0, 360, sweep_res)]

radial_sweep = cascaded_union(sweep_lines)

# load the input lines and combine them into one geometry
input_lines = fiona.open("input_lines.shp")
input_shapes = [shape(f['geometry']) for f in input_lines]
all_input_lines = cascaded_union(input_shapes)

perimeter = []
# traverse each radial sweep line and check for intersection with input lines
for radial_line in radial_sweep:
    inter = radial_line.intersection(all_input_lines)

    if inter.type == "MultiPoint":
       # radial line intersects at multiple points
       inter_dict = {}
       for inter_pt in inter:
           inter_dict[focal_pt.distance(inter_pt)] = inter_pt
       # save the nearest intersected point to the sweep centre point
       perimeter.append(inter_dict[min(inter_dict.keys())])

    if inter.type == "Point":
       # radial line intersects at one point only
       perimeter.append(inter)

    if inter.type == "GeometryCollection":
       # radial line doesn't intersect, so skip
       pass

# combine the nearest perimeter points into one geometry
solution = cascaded_union(perimeter)

# save the perimeter geometry
schema = {'geometry': 'MultiPoint', 'properties': {'test': 'int'}}
with fiona.open('perimeter.shp', 'w', 'ESRI Shapefile', schema) as e:
     e.write({'geometry':mapping(solution), 'properties':{'test':1}})
Rusty Magoo
fuente
Normalmente, un barrido radial no tiene una "resolución" significativa: escanea de un "evento" al siguiente en orden, donde los eventos consisten en los nodos originales de las polilíneas y sus intersecciones mutuas (que se pueden encontrar dinámicamente mientras se recorre el original) nodos) Su salida será perfectamente precisa.
whuber

Respuestas:

17

He reproducido tu ejemplo con shapefiles.

Puedes usar Shapely y Fiona para resolver tu problema.

1) Su problema (con un bien proporcionado Point):

ingrese la descripción de la imagen aquí

2) comenzando con una línea arbitraria (con una longitud adecuada):

from shapely.geometry import Point, LineString
line = LineString([(point.x,point.y),(final_pt.x,final_pt.y)])

ingrese la descripción de la imagen aquí

3) usando shapely.affinity.rotate para crear los radios (girando la línea desde el punto, mira también la respuesta de Mike Toews en Python, biblioteca bien formada: ¿es posible hacer una operación afín en el polígono de forma? ):

from shapely import affinity
# Rotate i degrees CCW from origin at point (step 10°)
radii= [affinity.rotate(line, i, (point.x,point.y)) for i in range(0,360,10)]

ingrese la descripción de la imagen aquí

4) ahora, usando shapely: cascaded_union (o shapely: unary_union ) para obtener un MultiLineString:

from shapely.ops import cascaded_union
mergedradii = cascaded_union(radii)
print mergedradii.type
MultiLineString

5) lo mismo con las líneas originales (shapefile)

import fiona
from shapely.geometry import shape
orlines = fiona.open("orlines.shp")
shapes = [shape(f['geometry']) for f in orlines]
mergedlines = cascaded_union(shapes)
print mergedlines.type
MultiLineString

6) se calcula la intersección entre las dos multigeometrías y el resultado se guarda en un shapefile:

 points =  mergedlines.intersection(mergedradii)
 print points.type
 MultiPoint
 from shapely.geometry import mapping
 schema = {'geometry': 'MultiPoint','properties': {'test': 'int'}}
 with fiona.open('intersect.shp','w','ESRI Shapefile', schema) as e:
      e.write({'geometry':mapping(points), 'properties':{'test':1}})

Resultado:

ingrese la descripción de la imagen aquí

7) pero, problema, si usa un radio más largo, el resultado es diferente:

ingrese la descripción de la imagen aquí

8) Y si desea obtener su resultado, debe seleccionar solo un punto con la distancia más corta desde el punto original en un radio:

points_ok = []
for line in mergeradii:
   if line.intersects(mergedlines):
       if line.intersection(mergedlines).type == "MultiPoint":
          # multiple points: select the point with the minimum distance
          a = {}
          for pt in line.intersection(merged):
              a[point.distance(pt)] = pt
          points_ok.append(a[min(a.keys())])
       if line.intersection(mergedlines).type == "Point":
          # ok, only one intersection
          points_ok.append(line.intersection(mergedlines))
solution = cascaded_union(points_ok)
schema = {'geometry': 'MultiPoint','properties': {'test': 'int'}}
with fiona.open('intersect3.shp','w','ESRI Shapefile', schema) as e:
     e.write({'geometry':mapping(solution), 'properties':{'test':1}})

Resultado final:

ingrese la descripción de la imagen aquí

Espero que sea lo que quieras.

gene
fuente
1
Excelente respuesta: particularmente informativa con respecto al uso de Fiona para entrada / salida a través de archivos de forma. He agregado un código a mi pregunta que usa su respuesta y lo modifiqué para reducir la cantidad de intersectioncálculos requeridos. Gracias.
Rusty Magoo