Trouver les segments de ligne les plus proches à pointer à l'aide de galbe?

Contexte

À partir d'un point connu, j'ai besoin d'établir le "périmètre visible" environnant le plus proche par rapport à une table de chaînes MultiLineStrings, comme indiqué sur le diagramme.

J'ai cherché sur ce site avec un certain nombre de termes (par exemple, bord minimum, périmètre minimum, voisin le plus proche, clip, contenant un polygone, visibilité, accrochage, nœuds coupés, trace de rayon, remplissage par inondation, limite intérieure, routage, coque concave) mais ne trouve aucune question précédente qui semble correspondre à ce scénario.

Diagramme

• Le cercle vert est le point connu.
• Les lignes noires sont les MultiLineStrings connues.
• Les lignes grises indiquent un balayage radial depuis le point connu.
• Les points rouges sont l'intersection la plus proche du balayage radial et des MultiLineStrings.

Paramètres

• Le Point ne coupera jamais les MultiLineStrings.
• Le Point sera toujours nominalement centré dans les MultiLineStrings.
• Les MultiLineStrings n'engloberont jamais complètement le Point, par conséquent le périmètre sera un MultiLineString.
• Il y aura un tableau contenant environ 1 000 MultiLineStrings (contenant normalement une seule ligne d'environ 100 points).

Méthodologie envisagée

• Effectuez un balayage radial en construisant une série de lignes à partir du point connu (par exemple, par incréments de 1 degré).
• Déterminez le point d'intersection le plus proche de chaque ligne de balayage radial avec les chaînes MultiLineStrings.
• Lorsqu'une des lignes de balayage radiales ne coupe aucune des chaînes MultiLineStrings, cela indiquerait un espace dans le périmètre qui serait logé dans la construction MultiLineString du périmètre.

Sommaire

Bien que cette technique trouve les intersections les plus proches, elle ne trouvera pas nécessairement tous les points de nœud de périmètre les plus proches, en fonction de la résolution du balayage radial. Quelqu'un peut-il recommander une méthode alternative pour établir tous les points du périmètre ou compléter la technique de balayage radial avec une forme de tampon, de sectorisation ou de compensation?

Logiciel

Ma préférence est d'utiliser SpatiaLite et / ou Shapely pour la solution, mais je serais heureux de recevoir toute suggestion qui pourrait être mise en œuvre à l'aide d'un logiciel open source.

Edit: Working Solution (basé sur la réponse 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}})

J'ai reproduit votre exemple avec des fichiers de formes.

Vous pouvez utiliser Shapely et Fiona pour résoudre votre problème.

1) Votre problème (avec un galbe Point):

2) en commençant par une ligne arbitraire (avec une longueur adéquate):

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

3) en utilisant shapely.affinity.rotate pour créer les rayons (en faisant tourner la ligne à partir du point, regardez également la réponse de Mike Toews sur Python, bibliothèque bien faite: est-il possible de faire une opération affine sur le polygone de forme? ):

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) maintenant, en utilisant bien faite: cascaded_union (ou bien faite: unary_union ) pour obtenir un MultiLineString:

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

5) la même chose avec les lignes d'origine (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) l'intersection entre les deux multigeométries est calculée et le résultat est enregistré dans un fichier de formes:

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

Résultat:

7) mais, problème, si vous utilisez un rayon plus long, le résultat est différent:

8) Et si vous voulez obtenir votre résultat, vous devez sélectionner uniquement le point avec la distance la plus courte du point d'origine sur un rayon:

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

Résultat final:

J'espère que c'est ce que tu veux.