5

I would like to find neighboured features in a shapefile using fiona and shapely. As shown in this answer https://gis.stackexchange.com/a/95595/62250 it's easy to get common borders between features.

Now that I have the common border of two features, how do I know which to which two features this border belongs?

How can I read out the attributes of those features?

5

As Fiona works with dictionaries, it is easy.

With the example of Finding the common borders between polygons in the same shapefile

enter image description here

import fiona
from shapely.geometry import shape, mapping
# load the shapefile as an iterable
polys = fiona.open("polygons.shp")

Now, with the itertools module:

import itertools
for pol1,pol2 in  itertools.combinations(polys, 2):
    if shape(pol1['geometry']).touches(shape(pol2['geometry'])):
        print  "Polygon ",pol1['properties']['id'],  " touches Polygon " , pol2['properties']['id']

Polygon 1 touches Polygon 4
Polygon 2 touches Polygon 4
Polygon 3 touches Polygon 4

And directly:

schema = {'geometry': 'LineString','properties':  OrderedDict([(u'poly1', 'int:10'), (u'poly2', 'int:10')])}
with fiona.open('intersection2.shp','w','ESRI Shapefile', schema) as output:
   with fiona.open("intersections.shp") as polys:
      for pol1,pol2 in  itertools.combinations(polys, 2):
          geompol1 = shape(pol1['geometry'])
          geompol2 = shape(pol2['geometry'])
          if geompol1.touches(geompol2):
             output.write({'geometry':mapping(geompol1.intersection(geompol2)),  'properties': {'poly1':pol1['properties']['id'],'poly2':pol2['properties']['id']}})

Control

 for line in fiona.open("intersection2.shp"):
     print line['properties']
 OrderedDict([(u'poly1', 1), (u'poly2', 4)])
 OrderedDict([(u'poly1', 2), (u'poly2', 4)])
 OrderedDict([(u'poly1', 3), (u'poly2', 4)])

enter image description here

  • Thanks! There's a typo in the last part of code, it should be print line['properties'] instead of print lines['properties']. – maz Jun 13 '16 at 10:58
2

Using the answer above (https://gis.stackexchange.com/a/197780/62250) I modified it to work with different crs and calculate the length of the intersection borders:

import fiona
from shapely.geometry import shape, mapping
import itertools
# load the shapefile as an iterable
polys = fiona.open("polygons.shp")

Writing the shapefile:

newschema = {'geometry': 'LineString','properties': {'id1':'str', 'id2':'str', 'length':'float'}}
with fiona.open('intersection.shp','w', crs=polys.crs, driver='ESRI Shapefile', schema=newschema) as e:
    for pol1,pol2 in  itertools.combinations(polys, 2):
        geompol1 = shape(pol1['geometry'])
        geompol2 = shape(pol2['geometry'])
        if geompol1.touches(geompol2):
            id1 = pol1['properties']['id']
            id2 = pol2['properties']['id']
            length = shape(mapping(geompol1.intersection(geompol2))).length
            e.write({'geometry':mapping(geompol1.intersection(geompol2)),  'properties': {'id1':id1, 'id2':id2, 'length':length}})
            print  "Polygon ", pol1['properties']['id'],  " touches Polygon ", pol2['properties']['id'], ", length of border: ", length

Testing the intersection shapefile:

for line in fiona.open("intersection.shp"):
    geom = shape(line['geometry'])
    print "Polygon ", line['properties']['id1'], " touches Polygon ", line['properties']['id2'], ", length of border: ", geom.length
  • geompol1.intersection(geompol2).convex_hull.length may work better if the intersection happens to be MultiPoint, and you are looking for a border (with at least two distinct points) – prusswan Jan 2 at 4:01

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.