Encuentre los segmentos de línea más cercanos al punto

15

Background

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

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

Diagram

  • El círculo verde es el Punto conocido.
  • Las líneas negras son las MultiLineStrings conocidas.
  • 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 MultiLineStrings.

Parámetros

  • ElpuntonuncaentrecruzarálosMultiLineStrings.
  • ElPuntosiempreestarácentradonominalmentedentrodelasMultiLineStrings.
  • LasMultiLineStringsnuncaencierrancompletamenteelPunto,porlotanto,elperímetroseráunaMultiLineString.
  • Habráunatablaquecontieneaproximadamente1,000MultiLineStrings(quenormalmentecontieneunasolalíneadeaproximadamente100puntos).

Metodologíaconsiderada

  • EmprendaunbarridoradialconstruyendounaseriedelíneasdesdeelPuntoconocido(en,porejemplo,enincrementosde1grado).
  • EstablezcaelpuntodeintersecciónmáscercanodecadalíneadebarridoradialconMultiLineStrings.
  • CuandounadelaslíneasdebarridoradialnoseintersectaconningunadelasMultiLineStrings,estoindicaríaunabrechaenelperímetroqueseacomodaríaenlaconstruccióndelperímetroMultiLineString.

Summary

Sibienestatécnicaencontrarálasinterseccionesmáscercanas,nonecesariamenteencontrarátodoslospuntosdelnododelperímetromáscercanos,dependiendodelaresolucióndelbarridoradial.¿Alguienpuederecomendarunmétodoalternativoparaestablecertodoslospuntosdelperímetroocomplementarlatécnicadebarridoradialconalgúntipodeamortiguamiento,sectorizaciónocompensación?

Software

MipreferenciaesusarSpatiaLitey/oShapelyparalasolución,peroagradeceríacualquiersugerenciaquepuedaimplementarseutilizandosoftwaredecódigoabierto.

Editar:Solucióndetrabajo([email protected])

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}})     
pregunta Rusty Magoo 30.12.2013 - 13:44

1 respuesta

15

He reproducido tu ejemplo con shapefiles.

Puede usar Shapely y Fiona para resolver tu problema.

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

2)comenzandoconunalíneaarbitraria(conunalongitudadecuada):

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

3)utilizando shapely.affinity.rotate para crear el radio (girando el línea desde el punto, busque también la respuesta de Mike Toews en Python, biblioteca bien formada: Es posible realizar 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)]

4)ahora,usando shapely: cascaded_union (o shapely:unary_union ) para obtener una 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:

7)pero,problema,siusaunradiomáslargo,elresultadoesdiferente:

8) Y si desea obtener su resultado, debe seleccionar solo el 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:

Espero que sea lo que quieres.

    
respondido por el gene 30.12.2013 - 17:59

Lea otras preguntas en las etiquetas