[GIS] Adding a Property (attribute) to a geometry in Shapely/Fiona

clipfionapython-2.7shapely

Once you have line and polygon geometries, intersection using Shapely is very straight forward. What I did after that is to integrate all the looped 'intersected' geometries into a MultiLineString.

Now my query is, while mapping the MultiLineString using Fiona, the output is a single geometry, with a single attribute.
If I could add individual attributes to each intersection geometry result, my MultiLineString would consist of multiple geometries.

How do I add a property or attribute to a LineString once the geometry is ready???

This problem is with reference to previous issues.

Step One – Vertical Lines in a Polygon Shapefile

Step Two – Clipping Line shapefiles to within the Extent of Polygon

The reason why I am doing this is because once I have the geometries of all the clipped lines, I need to extract the Node coordinates of each of them.

Best Answer

Shapely geometries are python classes, so you can simply set a new property in the object.

For my use-cases this has worked, though it may be worth checking for unintended side-effects depending on your scenario.

from shapely import geometry

p = geometry.Point(0, 0)
p.uid = 'unique point'
print(p.uid)  # prints: 'unique point'

with Fiona

import fiona

shape = fiona.open("my_shapefile.shp")
print shape.schema
{'geometry': 'LineString', 'properties': OrderedDict([(u'FID', 'float:11')])}
#first feature of the shapefile
first = shape.next()
print first # (GeoJSON format)
{'geometry': {'type': 'LineString', 'coordinates': [(0.0, 0.0), (25.0, 10.0), (50.0, 50.0)]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'FID', 0.0)])}

with PyShp

import shapefile

shape = shapefile.Reader("my_shapefile.shp")
#first feature of the shapefile
feature = shape.shapeRecords()[0]
first = feature.shape.__geo_interface__  
print first # (GeoJSON format)
{'type': 'LineString', 'coordinates': ((0.0, 0.0), (25.0, 10.0), (50.0, 50.0))}

with ogr:

from osgeo import ogr

file = ogr.Open("my_shapefile.shp")
shape = file.GetLayer(0)
#first feature of the shapefile
feature = shape.GetFeature(0)
first = feature.ExportToJson()
print first # (GeoJSON format)
{"geometry": {"type": "LineString", "coordinates": [[0.0, 0.0], [25.0, 10.0], [50.0, 50.0]]}, "type": "Feature", "properties": {"FID": 0.0}, "id": 0}

conversion to Shapely geometry (with the shape function)

from shapely.geometry import shape

shp_geom = shape(first['geometry']) # or shp_geom = shape(first) with PyShp
print shp_geom
LINESTRING (0 0, 25 10, 50 50)
print type(shp_geom)
<class 'shapely.geometry.linestring.LineString'>

computations
save the resulting shapefile

[GIS] Multi-part Geometrical Intersection error

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

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

Best Answer

Related Solutions

[GIS] Reading shapefile in Python

with Fiona

with PyShp

with ogr:

[GIS] Multi-part Geometrical Intersection error

New

Related Question