Znalezienie najbliższych segmentów linii do punktu za pomocą zgrabnego?

tło

Ze znanego punktu wymagam ustalenia najbliższego otaczającego „widocznego obwodu” względem tabeli MultiLineStrings, jak pokazano na schemacie.

Przeszukałem tę stronę z wieloma terminami (np. Minimalna krawędź, minimalny obwód, najbliższy sąsiad, klip, zawierający wielokąt, widoczność, przyciąganie, wycinanie węzłów, ray-trace, wypełnienie zalewowe, wewnętrzna granica, routing, wklęsły kadłub), ale nie mogę znaleźć żadnego poprzedniego pytania, które wydaje się pasować do tego scenariusza.

Diagram

• Zielony okrąg to znany Punkt.
• Czarne linie to znane linie MultiLineStrings.
• Szare linie oznaczają promieniowe przemiatanie od znanego punktu.
• Czerwone punkty to najbliższe przecięcie promieniowego przeciągnięcia i MultiLineStrings.

Parametry

• Punkt nigdy nie będzie przecinał MultiLineStrings.
• Punkt zawsze będzie nominalnie wyśrodkowany w MultiLineStrings.
• MultiLineString nigdy nie zamknie całkowicie Punktu, dlatego obwód będzie MultiLineString.
• Pojawi się tabela zawierająca około 1000 MultiLineStrings (zwykle zawierająca jedną linię około 100 punktów).

Rozważona metodologia

• Wykonaj przemiatanie promieniowe, budując serię linii ze znanego punktu (w, powiedzmy, co 1 stopień).
• Ustal najbliższy punkt przecięcia każdej promieniowej linii przeciągnięcia za pomocą MultiLineStrings.
• Gdy jedna z promieniowych linii przeciągnięcia nie przecina się z żadną z linii MultiLineString, oznaczałoby to przerwę na obwodzie, która byłaby umieszczona w obwodowej konstrukcji MultiLineString.

streszczenie

Chociaż ta technika znajdzie najbliższe przecięcia, niekoniecznie znajdzie wszystkie najbliższe punkty węzłów obwodowych, w zależności od rozdzielczości przesunięcia promieniowego. Czy ktoś może polecić alternatywną metodę ustalenia wszystkich punktów na obwodzie lub uzupełnić technikę przemiatania promieniowego jakąś formą buforowania, podziału na sektory lub kompensacji?

Oprogramowanie

Wolę używać SpatiaLite i / lub Shapely do rozwiązania, ale chętnie przyjmę wszelkie sugestie, które mogłyby zostać zaimplementowane przy użyciu oprogramowania typu open source.

Edycja: rozwiązanie robocze (na podstawie odpowiedzi @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}})

Odtworzyłem twój przykład z plikami kształtu.

Możesz użyć Shapely i Fiona, aby rozwiązać swój problem.

1) Twój problem (zgrabny Point):

2) zaczynając od dowolnej linii (o odpowiedniej długości):

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

3) używając shapely.affinity.rotate do utworzenia promieni (obracając linię od punktu, spójrz również odpowiedź Mike'a Toewsa na Python, foremna biblioteka: czy można wykonać operację afiniczną na wielokącie kształtu? ):

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

4) teraz, używając shapely: cascaded_union (lub shapely: unary_union ), aby uzyskać MultiLineString:

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

5) to samo z oryginalnymi liniami (plik kształtu)

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) obliczane jest przecięcie dwóch multigeometrii, a wynik zapisywany jest do pliku kształtu:

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

Wynik:

7) ale problem, jeśli użyjesz dłuższego promienia, wynik jest inny:

8) A jeśli chcesz uzyskać wynik, musisz wybrać tylko punkt o najmniejszej odległości od pierwotnego punktu na promieniu:

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

Ostateczny wynik:

Mam nadzieję, że tego właśnie chcesz.