[GIS] Multi-part Geometrical Intersection error

fionaintersectionpython-2.7shapely

I am working with intersection of polygon and line segments (continued from a previous question).

When a Line intersects with a Polygon, the result is a LineString;
When a LineRing intersects with a Line, the result is a MultiPoint.

While in most cases, the intersection is only in two points, the case I have encountered is multiple intersections (in figure).

Multiple intersections of a single Line Segment

The indicated yellow lines are the ones that intersect only at two points, and thus result in no errors. But when the program begins intersection of the LineString and LinearRing where there is multiple intersections, the following error shows up. 

NotImplementedError: Multi-part geometries do not themselves provide the array interface

Would there be any way to work around this? I am not at liberty to share the code, but any specific, additional information can be provided.

To provide some additional context, I am appending each intersection resultant MultiPoint and MultiLineString to a 'list' of previous results, respectively. So while the previous results were only a pair of XY coordinates, multiple intersections give multiple results. The line of code that results in the error is when I try to convert the 'list' of LineString (intersection) results to a MultiLineString using MultiLineString()

For instance, multi_line is the list of intersections resulting from line.intersection(polygon). This list needs to written to a MultiLineString for creating a new shapefile. 

MultipleLine = MultiLineString(multi_line)

Error

Best Answer

If you use Shapely, transform your polygon into a LinearRing or a LineString

from shapely.geometry import shape
import fiona
Multilines = MultiLineString([shape(line['geometry']) for line in fiona.open("lines.shp")])
Poly = shape(fiona.open("one_poly.shp").next()['geometry'])
Multilines.intersection(Poly)
<shapely.geometry.multilinestring.MultiLineString object at 0x1093285d0>

The intersection are polylines as you say

enter image description here

Now transform the polygon:

ring = LineString(list(poly.exterior.coords))
# or
ring = LinearRing(list(poly.exterior.coords))
Multilines.intersection(ring)
<shapely.geometry.multipoint.MultiPoint object at 0x109328990>

enter image description here

New

Without your complete script, I can help you but if I understand your question you want to convert a 'list' of intersections results to a MultiLineString or a MultiPoint, but you don't need this list.

If your Shapefile contains multiple lines, a solution is to open it with a for loop, adding the lines to a list and compute the intersections with another for loop which gives you your list.

But, knowing that the result of the intersection of a MultiLineString with a Polygon is a MultiLineString and the intersection of a MultiLineString with a LinearRing is a MultiPoint, it is easiest to open your shapefile with:

Directly a MultiLineString

multilines1 = MultiLineString([shape(line['geometry']) for line in fiona.open("lines.shp")])

A list of LineStrings

multilines2 = [shape(line['geometry']) for line in fiona.open("lines.shp")]

A list of the coordinates of the LineStrings

multilines3 = [line['geometry']['coordinates'] for line in fiona.open("lines.shp")]

And with:

ring =  LineString(list(Poly.exterior.coords))

You compute directly the intersections as MultilineStrings or MultiPoints with

multi_line = Multilines1.intersection(Poly)  # MultiLineString
multi_point = Multilines1.intersection(ring) # MultiPoint

or 

multi_line = MultiLineString(Multilines2).intersection(Poly) # MultiLineString
multi_point = MultiLineString(Multilines2).intersection(ring) # MultiPoint
multi_line = MultiLineString(Multilines3).intersection(Poly) # MultiLineString
multi_line = MultiLineString(Multilines3).intersection(ring) # MultiPoint

If you have multiple polygons, use a MultiPolygon

multipoly = MultiPolygon([shape(poly['geometry']) for poly in fiona.open("multipolygons.shp")])
multi_line = Multilines1.intersection(multipoly)  # MultiLineString
multiring = MultiLineString([list(shape(poly['geometry']).exterior.coords) for poly in fiona.open("multipolygons.shp")])
multi_point = Multilines1.intersection(multiring)  # MultiPoint
Related Question