17

Background

From a known Point, I require to establish the nearest surrounding "visible perimeter" against a table of MultiLineStrings, as shown on the diagram.

I've searched this site with a number of terms (e.g. minimum edge, minimum perimeter, nearest neighbour, clip, containing polygon, visibility, snap, cut nodes, ray-trace, flood fill, inner boundary, routing, concave hull) but cannot find any previous question which seems to match this scenario.

Diagram

  • The green circle is the known Point.
  • The black lines are the known MultiLineStrings.
  • The grey lines are an indication of a radial sweep from the known Point.
  • The red points are the nearest intersection of the radial sweep and the MultiLineStrings.

enter image description here

Parameters

  • The Point will never intersect the MultiLineStrings.
  • The Point will always be nominally centred within the MultiLineStrings.
  • The MultiLineStrings will never fully enclose the Point, therefore the perimeter will be a MultiLineString.
  • There will be a table containing approximately 1,000 MultiLineStrings (normally containing a single line of about 100 points).

Methodology Considered

  • Undertake a radial sweep by constructing a series of lines from the known Point (at, say, 1 degree increments).
  • Establish the nearest intersection point of each radial sweep line with the MultiLineStrings.
  • When one of the radial sweep lines doesn't intersect with any of the MultiLineStrings, this would indicate a gap in the perimeter which would be accommodated in the perimeter MultiLineString construction.

Summary

Whilst this technique will find the nearest intersections, it won't necessarily find all the closest perimeter node points, dependent on the resolution of the radial sweep. Can anyone recommend an alternative method to establish all the perimeter points or supplement the radial sweep technique with some form of buffering, sectoring or offsetting?

Software

My preference is to use SpatiaLite and/or Shapely for the solution but would welcome any suggestions which could be implemented using open source software.

Edit: Working Solution (based on answer by @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}})
  • Normally a radial sweep has no meaningful "resolution": it scans from one "event" to the next in order, where events consist of the original nodes of the polylines and of their mutual intersections (which can be found dynamically while sweeping around the original nodes). Its output will be perfectly accurate. – whuber Dec 30 '13 at 17:46
17

I have reproduced your example with shapefiles.

You can use Shapely and Fiona to solve your problem.

1) Your problem (with a shapely Point):

enter image description here

2) starting with an arbitrary line (with an adequate length):

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

enter image description here

3) using shapely.affinity.rotate to create the radii (rotating the line from the point, look also the Mike Toews's answer at Python, shapely library: is it possible to do an affine operation on shape polygon?):

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

enter image description here

4) now, using shapely:cascaded_union (or shapely:unary_union) to get a MultiLineString:

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

5) the same with the original lines (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) the intersection between the two multigeometries is computed and the result is saved to a 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}})

Result:

enter image description here

7) but, problem, if you use a a longer radius, the result is different:

enter image description here

8) And if you want to get your result, you need to select only the point with the shortest distance from the original point on a radius:

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

Final result:

enter image description here

I hope that is what you want.

  • 1
    Excellent answer - particularly informative with respect to Fiona usage for input/output via shapefiles. I've appended some code to my question which uses your answer and modified it to reduce the amount of intersection calculations required - thank you. – Rusty Magoo Dec 30 '13 at 21:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.