[GIS] Error with Python Script using Shapely

fionapythonshapely

My professor gave us a basic Python script that demonstrated how Shapely works with a shapefile, however it's not working for me. Here is the code:

from shapely.geometry import mapping, shape
from fiona import collection

with collection("hawaii_vector_data.shp","r") as input:
  for polygon in input:
    c = polygon['properties']['class']
    if c == 0 or c == None:
      continue
    s = shape(polygon['geometry'])
    print s

and when I run it I get this Traceback error:

Traceback (most recent call last):
  File "test.py", line 6, in <module>
    c = polygon['properties']['class']
KeyError: 'class'

Since I am new to using Shapely I have no idea what the error might be. Any advice would be great.

Best Answer

Fiona uses dictionaries to read shapefiles

1) the keys of the dictionary = schema of the shapefile

shapes = fiona.open("polygons.shp")
# schema of the shapefile -> the keys of the dictionary
print shapes.schema
{'geometry': 'Polygon', 'properties': OrderedDict([(u'id', 'int:10')])}
# read first feature of the shapefile
first = shapes.next()
print first
{'geometry': {'type': 'Polygon', 'coordinates': [[(151116.87238259654, 135890.8706318218), (153492.19971554304, 134793.3055883224), (153934.50204650551, 133892.31935858406), (152623.97662143156, 131811.86024627919), (150903.91200102202, 130894.49244872745), (149347.66305874675, 132991.33312884573), (149151.08424498566, 134383.76639298678), (151116.87238259654, 135890.8706318218)]]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', 2)])}

2) the geometry key

print first['geometry'] # GeoJSON format
{'type': 'Polygon', 'coordinates': [[(151116.87238259654, 135890.8706318218), (153492.19971554304, 134793.3055883224), (153934.50204650551, 133892.31935858406), (152623.97662143156, 131811.86024627919), (150903.91200102202, 130894.49244872745), (149347.66305874675, 132991.33312884573), (149151.08424498566, 134383.76639298678), (151116.87238259654, 135890.8706318218)]]}
# the coordinates of the geometry
print ['geometry']["coordinates"]
[[(151116.87238259654, 135890.8706318218), (153492.19971554304, 134793.3055883224), (153934.50204650551, 133892.31935858406), (152623.97662143156, 131811.86024627919), (150903.91200102202, 130894.49244872745), (149347.66305874675, 132991.33312884573), (149151.08424498566, 134383.76639298678), (151116.87238259654, 135890.8706318218)]]

If you want to transform the geometry to a Shapely geometry

poly = shape(first['geometry']
print poly
POLYGON ((151116.8723825965 135890.8706318218, 153492.199715543 134793.3055883224, 153934.5020465055 133892.3193585841, 152623.9766214316 131811.8602462792, 150903.912001022 130894.4924487275, 149347.6630587468 132991.3331288457, 149151.0842449857 134383.7663929868, 151116.8723825965 135890.8706318218))

3) the attributes -> property key

print first['property']
OrderedDict([(u'id', 2)])
# value of attribute id
print fist['properties']['id']
2

4) you problem

with fiona.open("polygons.shp") as input:
    for polygon in input:
        c = polygon['properties']['id']
        s = shape(polygon['geometry'])
        print c, s
....

And you can see that the error in your script comes from the fact that there is no "class" key in the properties of your shapefile.

If you want to select some features:

with fiona.open("polygons.shp") as input:
    for polygon in input:
        if polygon['properties']['id'] > 0:
             s = shape(polygon['geometry'])
             print c, s

Related Solutions

Shapely Python – Looping Through a Number of Polygons

Some solutions according to the position of an element in a list:

polA = Polygon([(0,0), (3,0), (3,3), (0,3)])
polB = Polygon([(2,-1), (5,-1), (5,2), (2,2)])
polC = Polygon([(5,2), (8,2), (8,5), (5,5)]) 
collection = [polA, polB, polC]

Iterating by index:

for i in range(len(collection)-1):
    print collection[i], collection[i+1], collection[i].touches(collection[i+1])
POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) False
POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) True

Iterating with itertools (standard module), all combinations without repeated elements

import itertools
for pol in  itertools.combinations(collection, 2):
    print pol[0],pol[1], pol[0].touches(pol[1])
 POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) False
 POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) False
 POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) True

or with the suggestion of Mike T

for pol1, pol2 in  itertools.combinations(collection, 2):
     print pol1, pol2, pol1.touches(pol2)
 POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) False
 POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) False
 POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) True

Compare with

for pol in  itertools.permutations(collection, 2):
    print pol[0],pol[1], pol[0].touches(pol[1]) 
POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) False
POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) False
POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) False
POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) True
POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) False
POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) True

Look also at Python - Previous and next values inside a loop

from itertools import tee, islice, chain, izip
def previous_and_next(some_iterable):
    prevs, items, nexts = tee(some_iterable, 3)
    prevs = chain([None], prevs)
    nexts = chain(islice(nexts, 1, None), [None])
    return izip(prevs, items, nexts)
 for previous, item, nxt in previous_and_next(collection):
      print "Item is now", item, "next is", nxt, "previous is", previous

 Item is now POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) next is POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) previous is None
 Item is now POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) next is POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) previous is POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0))
 Item is now POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) next is None previous is POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1))

Python Shapely Fiona – Dissolving Polygons Based on Attributes with Python

The question is about Fiona and Shapely and the other answer using GeoPandas requires to also know Pandas. Moreover GeoPandas uses Fiona to read/write shapefiles.

I do not question here the utility of GeoPandas, but you can do it directly with Fiona using the standard module itertools, specially with the command groupby ("In a nutshell, groupby takes an iterator and breaks it up into sub-iterators based on changes in the "key" of the main iterator. This is of course done without reading the entire source iterator into memory", itertools.groupby).

Original Shapefile coloured by the STATEFP field

enter image description here

from shapely.geometry import shape, mapping
from shapely.ops import unary_union
import fiona
import itertools
with fiona.open('cb_2013_us_county_20m.shp') as input:
    # preserve the schema of the original shapefile, including the crs
    meta = input.meta
    with fiona.open('dissolve.shp', 'w', **meta) as output:
        # groupby clusters consecutive elements of an iterable which have the same key so you must first sort the features by the 'STATEFP' field
        e = sorted(input, key=lambda k: k['properties']['STATEFP'])
        # group by the 'STATEFP' field 
        for key, group in itertools.groupby(e, key=lambda x:x['properties']['STATEFP']):
            properties, geom = zip(*[(feature['properties'],shape(feature['geometry'])) for feature in group])
            # write the feature, computing the unary_union of the elements in the group with the properties of the first element in the group
            output.write({'geometry': mapping(unary_union(geom)), 'properties': properties[0]})

Result

enter image description here

Best Answer

Related Solutions

Shapely Python – Looping Through a Number of Polygons

Python Shapely Fiona – Dissolving Polygons Based on Attributes with Python

Related Question