Encontrando segmentos de linha mais próximos para apontar usando o shapely?

fundo

De um ponto conhecido, eu preciso estabelecer o "perímetro visível" ao redor mais próximo de uma tabela de MultiLineStrings, conforme mostrado no diagrama.

Pesquisei este site com vários termos (por exemplo, borda mínima, perímetro mínimo, vizinho mais próximo, clipe, contendo polígono, visibilidade, snap, nós de corte, rastreamento de raios, preenchimento de inundação, limite interno, roteamento, casco côncavo), mas Não é possível encontrar nenhuma pergunta anterior que pareça corresponder a esse cenário.

Diagrama

• O círculo verde é o ponto conhecido.
• As linhas pretas são as MultiLineStrings conhecidas.
• As linhas cinzas são uma indicação de uma varredura radial do ponto conhecido.
• Os pontos vermelhos são a interseção mais próxima da varredura radial e das MultiLineStrings.

Parâmetros

• O ponto nunca cruzará as MultiLineStrings.
• O ponto sempre será nominalmente centralizado nas MultiLineStrings.
• As MultiLineStrings nunca incluirão totalmente o Point; portanto, o perímetro será uma MultiLineString.
• Haverá uma tabela contendo aproximadamente 1.000 MultiLineStrings (normalmente contendo uma única linha de cerca de 100 pontos).

Metodologia considerada

• Realize uma varredura radial construindo uma série de linhas a partir do ponto conhecido (em, digamos, incrementos de 1 grau).
• Estabeleça o ponto de interseção mais próximo de cada linha de varredura radial com as MultiLineStrings.
• Quando uma das linhas de varredura radial não se cruza com nenhuma das MultiLineStrings, isso indica uma lacuna no perímetro que seria acomodada na construção do MultiLineString de perímetro.

Sumário

Embora essa técnica encontre as interseções mais próximas, não encontrará necessariamente todos os pontos de nó de perímetro mais próximos, dependendo da resolução da varredura radial. Alguém pode recomendar um método alternativo para estabelecer todos os pontos de perímetro ou complementar a técnica de varredura radial com alguma forma de amortecimento, setorização ou compensação?

Programas

Minha preferência é usar SpatiaLite e / ou Shapely para a solução, mas gostaria de receber sugestões que possam ser implementadas usando software de código aberto.

Editar: Solução de Trabalho (com base na resposta 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}})

Eu reproduzi seu exemplo com shapefiles.

Você pode usar Shapely e Fiona para resolver seu problema.

1) O seu problema (com um formato Point):

2) começando com uma linha arbitrária (com um comprimento adequado):

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

3) usando shapely.affinity.rotate para criar os raios (girando a linha a partir do ponto, veja também a resposta de Mike Toews na biblioteca bem torneada do Python: é possível fazer uma operação afim no polígono da 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)]

4) agora, usando shapely: cascaded_union (ou shapely: unary_union ) para obter um MultiLineString:

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

5) o mesmo com as linhas originais (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) a interseção entre as duas multigeometrias é calculada e o resultado é salvo em um arquivo de forma:

 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:

7) mas, problema, se você usar um raio mais longo, o resultado será diferente:

8) E se você deseja obter o resultado, precisa selecionar apenas o ponto com a menor distância do ponto original em um raio:

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:

Espero que seja isso que você quer.