2

I would like to calculate the total length of a line network, similar to what is done in this QGIS tutorial.

Here is what I have:

import fiona
from shapely.geometry import LineString

length = 0.0
with fiona.open('/Users/felix/data/centerline.shp', 'r') as shp:
    print(shp.crs)
    for rec in shp:
        line = LineString(rec['geometry']['coordinates'])
        length += line.length

While this method works fine for some shapefiles, for others it produces this error:

{'datum': 'NAD27', 'proj': 'lcc', 'y_0': 0, 'lat_2': 35.46666666666667, 'lat_0': 33.5, 'no_defs': True, 'x_0': 609601.2192024384, 'lon_0': -118, 'lat_1': 34.03333333333333, 'units': 'us-ft'}
Traceback (most recent call last):
  File "/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/shapely/geometry/linestring.py", line 200, in geos_linestring_from_py
    array = ob.__array_interface__
AttributeError: 'list' object has no attribute '__array_interface__'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/felix/Documents/Scripts/python/county_workflow_automation/scratch.py", line 11, in <module>
    line = LineString(rec['geometry']['coordinates'])
  File "/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/shapely/geometry/linestring.py", line 49, in __init__
    self._set_coords(coordinates)
  File "/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/shapely/geometry/linestring.py", line 94, in _set_coords
    self._geom, self._ndim = geos_linestring_from_py(coordinates)
  File "/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/shapely/geometry/linestring.py", line 271, in geos_linestring_from_py
    assert n == 2 or n == 3
AssertionError

Can somebody help me interpret this error? Or, if possible, is there a more direct way to accomplish this task?

  • I think it is a projection related error..I just tested your code against my data, it worked well. – SIslam Nov 10 '16 at 6:25
2

What is your version of Shapely ?

Because the error comes from the lines 200-201 of the Shapely file linestring.py

    array = ob.__array_interface__
    assert len(array['shape']) == 2

With a valid LineString

print(ob)
LINESTRING (270818.4477341335 4458997.238866236, 270833.274660459 4458983.162670357, 270830.8347865067 4458975.280000665, 270822.3890689795 4458967.209648361, 270823.3274820381 4458959.702343893, 270822.7644342029 4458958.013200387)
array = ob.__array_interface__
print array
{'version': 3, 'data': <shapely.coords.c_double_Array_12 object at 0x103ba9158>, 'typestr': '<f8', 'shape': (6, 2)}

With no error

array = ob.array_interface
AttributeError: 'list' object has no attribute 'array_interface'

You can control if the LineStrings are valids

for line in fiona.open('lines_feet.shp'):
    if LineString(line['geometry']['coordinates']).__array_interface__:
        print("ok")
    ok
    ok
    ...

But the easiest way is to use the shape function of Shapely. The resulting unit is feet ('units': 'us-ft' in crs)

import fiona
length = 0.0
with fiona.open('lines_feet.shp', 'r') as shp:
print(shp.crs)
for line in shp:
   length += shape(line['geometry']).length
print("{:.2f} in {}".format(length,'feets'))
{'datum': 'NAD27', 'proj': 'lcc', 'y_0': 0, 'lat_2': 35.46666666666667, 'lat_0': 33.5, 'no_defs': True, 'x_0': 609601.2192024384, 'lon_0': -118, 'lat_1': 34.03333333333333, 'units': 'us-ft'}
264.14 in feets

or in one line

length= sum([shape(line['geometry']).length  for  line in fiona.open("lines_feet.shp")])
print("{:.2f} in {}".format(length,'feets'))
264.14 in feets
  • I am using shapely version 1.5.17. It appears as though I was attempting to pass MultiLineString geometries into the LineString constructor. Not sure how this happened, but your reference to shapely.geometry.shape() fixed my issue.Thank you very much for your response, gene. – Jess Nov 11 '16 at 1:47

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.