I have two tables of roads - the same roads, from different sources, and am trying to compare how similar their positions are. I'm wanting to buffer each road in one of the feature classes and calculate the % overlap between it and the corresponding road from the other feature class. Could obviously do this manually but since there are hundreds of roads would rather use Python. The thing is, I'm new to it. Any help on how to automate using the intersect tool between objects from different tables would be very welcome. The objects will have an attribute in common.

  • 3
    Very closely related: gis.stackexchange.com/questions/49923/…. Perhaps this is the same question as yours, Sally? – whuber Feb 16 '13 at 13:40
  • 1
    What GIS software are you using? – RyanDalton Feb 16 '13 at 18:26
  • 1
    What version of your GIS software are you using? – PolyGeo Feb 17 '13 at 4:42
  • I'm using ArcGIS. Thanks Seba, I think there's a danger of calculating intersections between the wrong roads though, if one of them has been mapped poorly. That's why I'd like to specify which objects in the feature class I'd like to perform the intersection on, can't figure out how to do this though. – user15261 Feb 17 '13 at 16:07

I am assuming you are using ArcGIS 10.1 with this answer. I also assume that there is a 1:1 relationship between the two feature classes, that they are both in a geodatabase, that they have the same spatial reference, and that they are already sorted in the same order. I assume that if the buffers of two pairs of features do not intersect at all that you do not care about that information being output. If any of this is not true then some adjustments will need to be made.

This script will create a pair of buffered polygon feature classes as well as an intersection feature class and will attribute it with the percent overlap between pairs of features as well as FIDs for joining back to the buffered polygon feature classes (or your original lines if they share the same FIDs). The pairing of features is only made by simultaneously iterating over both of the buffered polygon feature classes (using itertools.izip) and assuming that they are pre-sorted and 1:1. I do have an assertion that will result in an error message if the values of the common key field you mentioned are not equal.

The if __name__ == "__main__": section is where you configure the script using your own data. This is a common way to test functions in a Python module.

import arcpy, os, itertools

def computeIntersectOverlapPairs(feature_class_1, key_field_1, feature_class_2, key_field_2, output_feature_class, buffer_distance):
    desc = arcpy.Describe(feature_class_1)
    buffer1 = arcpy.Buffer_analysis(feature_class_1, os.path.join(os.path.dirname(output_feature_class), "{0}_Buffer".format(os.path.basename(feature_class_1))), buffer_distance)
    buffer2 = arcpy.Buffer_analysis(feature_class_2, os.path.join(os.path.dirname(output_feature_class), "{0}_Buffer".format(os.path.basename(feature_class_2))), buffer_distance)
    createOutputFC(output_feature_class, desc.SpatialReference)
    with arcpy.da.SearchCursor(buffer1, ["OID@", "SHAPE@", key_field_1]) as rows1:
        with arcpy.da.SearchCursor(buffer2, ["OID@", "SHAPE@", key_field_2]) as rows2:
            with arcpy.da.InsertCursor(output_feature_class, ["SHAPE@", "FID_1", "FID_2", "OVERLAP"]) as outrows:
                for row1, row2 in itertools.izip(rows1, rows2):
                    assert row1[2] == row2[2], "Key values must be equal! (FID1: {0}, FID2 {1}, Key1: {2}, Key2: {3})".format(row1[0], row2[0], row1[2], row2[2])
                    shape1 = row1[1]
                    shape2 = row2[1]
                    if not shape1.disjoint(shape2):
                        intersect_shape = shape1.intersect(shape2, 4)
                        percent_overlap = (intersect_shape.area / (shape1.area + shape2.area)) * 100
                        outrows.insertRow((intersect_shape, row1[0], row2[0], percent_overlap))

def createOutputFC(path, spatial_reference=None):
    arcpy.CreateFeatureclass_management(os.path.dirname(path), os.path.basename(path), "POLYGON", spatial_reference=spatial_reference)
    arcpy.AddField_management(path, "FID_1", "LONG")
    arcpy.AddField_management(path, "FID_2", "LONG")
    arcpy.AddField_management(path, "OVERLAP", "DOUBLE")

if __name__ == '__main__':
    arcpy.env.overwriteOutput = True
    inputfc1 = r"C:\GISData\test8.gdb\atlantic_hurricanes_2000_lines"
    inputfc2 = r"C:\GISData\test8.gdb\atlantic_hurricanes_2000_lines_shifted"
    outputfc = r"C:\GISData\test8.gdb\atlantic_hurricanes_2000_lines_overlap"
    field1 = "EVENTID"
    field2 = "EVENTID"
    bufferdist = 0.1
    computeIntersectOverlapPairs(inputfc1, field1, inputfc2, field2, outputfc, bufferdist)

Visual analysis

Input Lines: Input Lines

Buffered Lines: Buffered Lines

Overlapping Area (~12% in this case): Overlapping Area


You could:

  1. calculate buffers for the two shapes (choose "dissolve field" as NONE, so every feature/road gets its own buffer...
  2. intersect the two buffer-shapes
  3. add field in resulting shapes attribute table
  4. calculate the area of this shape
  5. use the calculated area to create an overlap-percentage ranking (the smaller the area, the higher the overlap)

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.