# Finding differences between two shapefiles in Python

I'm trying to find out whether there were any border changes between two yearly shapefiles. I've heard a lot of praise for Python's open source GIS tools like Shapely/Fiona, and thought I would try those, but I'm at a loss for how to go about implementing this. Is Shapely poorly suited for this task, and, if so, what tool would be better? Any advice you have would be greatly appreciated. Thanks!

• What about using the object.equals method? What else have you tried? – Fezter Dec 10 '13 at 22:43
• Do you mean with "border changes" the extent of the shapefile? – ustroetz Dec 10 '13 at 23:09
• @ustroetz - sorry, I should have explained better. Each shapefile is a map of school districts, and I want to see whether any of those districts change their borders between the two shapefiles. So it is not the extent of the shapefile, but rather each individual district's borders. – maxliving Dec 11 '13 at 16:29
• @Fezter - My understanding is that the object.equals method compares two polygons, but what I want is to compare a whole set of polygons to another set (see the clarification of my question above). But it's quite possible I'm misunderstanding something here. – maxliving Dec 11 '13 at 16:31
• What you could do is get all the vertices of the two shapefiles and compare them. If they are not identical, there were changes. Would that work? – ustroetz Dec 11 '13 at 17:23

Suppose we have two polygons (green and blue):

They are not equal (as Fetzer suggest):

``````green.equals(blue)
False
``````

and

``````blue.equals(green)
False
``````

And we can can determine the difference (in red):

``````blue.difference(green)
``````

and

``````green.difference(blue)
``````

gives an empty geometry

Thus, you can use a supplementary condition:

``````if not poly1.difference(poly2).is_empty:
process
``````

And if you want to find the nodes that have been modified:

`````` S1 = set(list(blue.difference(green).exterior.coords)
S2 = set(list(blue.exterior.coords)
S3 = set(list(green.exterior.coords)
``````

S1 - S2 and S1 - S3 gives the points (two here blue and red):

and the distance:

``````point1.distance(point2)
``````

## New : compare multiple polygons:

Here is one solution:
For that, I use Fiona to open the polygon shapefiles and save a resulting shapefile with differences:

``````import fiona
from shapely.geometry import shape
green = fiona.open("poly1.shp")
blue = fiona.open("poly2.shp")
# test the function difference between green and blue shapefiles
[not shape(i['geometry']).difference(shape(j['geometry'])).is_empty for i,j in zip(list(green),list(blue))]
[False, False, False, True]
# control
for geom in [shape(i['geometry']).difference(shape(j['geometry'])) for i,j in zip(list(green),list(blue))]:
print geom
GEOMETRYCOLLECTION EMPTY
GEOMETRYCOLLECTION EMPTY
GEOMETRYCOLLECTION EMPTY
POLYGON ((-0.0806077747083538 0.6329375155131045, 0.0085568963219002 0.5081069760707490, -0.0816567708381215 0.6025166277498414, -0.1529885076623247 0.5437728444828506, -0.1292856235630944 0.6206937720158269, -0.0806077747083538 0.6329375155131045))
# test the function difference between blue and green shapefiles
[not shape(i['geometry']).difference(shape(j['geometry'])).is_empty for i,j in zip(list(blue),list(green))]
[True, False, False, False]
# control
for geom in [shape(i['geometry']).difference(shape(j['geometry'])) for i,j in zip(list(blue),list(green))]:
print geom
POLYGON ((0.2292711598746081 0.7386363636363635, 0.6691026645768023 0.6691026645768023, 0.2440830721003134 0.7205329153605015, 0.1074843260188087 0.3452978056426331, 0.2292711598746081 0.7386363636363635))
GEOMETRYCOLLECTION EMPTY
GEOMETRYCOLLECTION EMPTY
GEOMETRYCOLLECTION EMPTY
# thus you can write a resulting shapefile withe the differences
from shapely.geometry import mapping
schema = {'geometry': 'Polygon','properties': {'test': 'int'}}
with fiona.open('diff.shp','w','ESRI Shapefile', schema) as e:
for geom  in [shape(i['geometry']).difference(shape(j['geometry'])) for i,j in zip(list(green),list(blue)]:
if not geom.is_empty:
e.write({'geometry':mapping(geom), 'properties':{'test':1}})
for geom  in [shape(i['geometry']).difference(shape(j['geometry'])) for i,j in zip(list(blue),list(green))]:
if not geom.is_empty:
e.write({'geometry':mapping(geom), 'properties':{'test':2}})
``````

Result:

• This makes sense for a given polygon, but if I have two sets of polygons and want to see which ones are different, how do I do that? Unless I can get an identifier to match them across shapefiles (which may be a possibility, but I'd prefer to figure out a solution that doesn't require that). – maxliving Dec 15 '13 at 2:12
• see above for a solution – gene Dec 15 '13 at 12:43
• Thank you! This looks like it will do what I need. My only concern is how quickly/slowly it will run on my real data, but I'll test it out and see. – maxliving Dec 15 '13 at 14:24