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.
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.
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.
Problem solved. Stitched together from three scripts online.
If all goes according to plan, the end result will look something like this:
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.
What is Qgis?
LikeLike
It’s a free open-source mapping software! Kind of like Photoshop + Google Earth, but you can play with it in Python.
LikeLiked by 1 person
ha sounds interesting
LikeLiked by 1 person
Sounds complicated but interesting. It intrigued me to read more 🙂
LikeLiked by 1 person
Awww, thanks Bernadette! I’m hoping to write more on QGIS in the coming months.
LikeLike