Python: How to List Polygon Intersections in QGIS

Header

You know what I love about programming? Once you know a little bit, you can Frankenstein your way out of a lot of problems.

Recently, I wanted to have a way to check two shapefiles in QGIS to see if their features overlap. I know what you’re thinking—there’s an intersect tool for that. But it just creates a new layer of the intersected areas. I wanted a list of which specific features were intersecting which other specific features in two separate layers.

Google pointed me at this thread on StackExchange which has the bones for setting up loop that’ll go through each feature in two separate layers.

Step 1: Set Up the Preliminary Variables

The print layers lines are optional to determine which layers are being run through the loop. Take it from someone who spent an hour trying to figure out why it wasn’t working, those lines can be real damn helpful if you’ve accidentally selected the wrong layers.

Also, from what I’ve been able to gather, the layer number corresponds with the layer’s position in the stack. Start at the bottom with O and work your way up. Deselected layers won’t register.

Now for the loop innards. In this case, we just want QGIS to grab the names of the polygons and add them to their respective lists if they intersect. This thread explains how to access the attributes where the names of your shapes live.

Step 2: Put the Attribute Values into the Lists

Once that’s done, we can go back to Automate the Boring Stuff With Python by Al Sweigart to write it into a .txt file. Remember that if you want to print a variable that is a string to the Python console, you need to place it after a % within quotation marks. But if you’re writing it into a txt file, you can just plop the variable in and concatenate any other string values with single quotes.

Step 3

Problem solved. Stitched together from three scripts online.

If all goes according to plan, the end result will look something like this:

The Result

The Code:

mapcanvas = iface.mapCanvas()  #Put the map in a variable called mapcanvas
layers = mapcanvas.layers()  #Put the layers in a variable called layers
pselections=[]  #Declare a Polygon 1 List
pmselections=[]  #Declare a Polygon 2 through 5 list

for w in layers[0].getFeatures():  #For each feature in the first layer. That will be called W
for s in layers[1].getFeatures(): #Go through each layer in the second layer. That will be called S
if s.geometry().intersects(w.geometry()): #If the geometry of S and W intersect
pmselections.append(s.attributes()[0])  #Get the Name in Attribute Column 0
pselections.append(w.attributes()[0])  #Put those into the Lists
break #only one or less intersection are possible

testfile = open('C:\\Users\\Desktop\\IntersectTest.txt', 'w')
i=0
for x in pselections:
print "%s" %x
print " intersects %s" %pmselections[i]
testfile.write(x + ' intersects ' + pmselections[i] + '\n')
i+=1
testfile.close()

print "Done"

I’m still attempting to wrangle the data into a csv file, but that’s a different week’s problem.

Advertisements
Tagged , , , , , , , ,

3 thoughts on “Python: How to List Polygon Intersections in QGIS

  1. ARJ says:

    What is Qgis?

    Like

  2. ARJ says:

    ha sounds interesting

    Liked by 1 person

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: