Find common borders of features and get their attributes using fiona and shapely

Question

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?

Answer 1

As Fiona works with dictionaries, it is easy.

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

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)])

Answer 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