2

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
13

Suppose we have two polygons (green and blue):

enter image description here

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)

enter image description here

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

enter image description here

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:

enter image description here

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:

enter image description here

  • 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

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.