So here’s a quick thing to note before we get started. Python will throw a shit-fit if you do not have everything properly formatted. Indentation matters. I can’t figure out how to make WordPress format my tabs, so if you do copy and paste the below code, indent the loops to avoid errors.
The Boring Details*:
The intersect analysis takes two input features and generates an output with only the overlap between the two. For example, I have a set of points in random locations around California. That’s input #1. I want to know which county they are located within. The shapefile with all the counties in California—that’s input #2. The output will be a shapefile with the points from input #1 including the fields in the attribute table from the overlapping county in input #2. It defaults to a point type because that’s the lowest dimension geometry of our inputs.
*It is not absolutely necessary to understand all that before you get started—I didn’t really get how it worked until after I started running the script on different shapefiles. Just know that the intersect analysis makes an attribute table that combines the fields of both input features. Like so:
The Script
The script requires two parameters. The first is going to be the name of your points, the second will be the polygons (counties). So when we write the Python, we can assume two string inputs from the user.
point = arcpy.GetParameterAsText(0) polygon = arcpy.GetParameterAsText(1)
Then we need to set up the intersect analysis:
arcpy.Intersect_analysis([point, polygon],r"c:\\Users\\Desktop\\Test.shp")
The two inputs need to be enclosed in brackets. The filepath parameter at the end is a temporary location for ArcGIS to store the results of the analysis. I just have it mapped to my desktop because I’m lazy, but it could go anywhere.
Onward. This is how you get Python to create the text document. I also have that document mapped to my desktop. For the purposes of this script, I’ve assigned that to variable f. PythonForBeginners has a really good introduction to file handling if you want to read more about how that works.
import os f = open(r"c:\\Users\\Desktop\\TestFile.txt","w")
Now, you’re probably wondering how to access the data in the attribute table. For that, you’ll use a cursor. You can loop through all the rows in your newly created intersect shapefile and have ArcPy spit out specific fields.
with arcpy.da.SearchCursor("c:\\Users\\Desktop\\Test.shp",['WAYPT', 'COUNTY']) for row in cursor: f.write(str(row) + '\n')
SearchCursor needs the names of the fields you wish to keep. In the attribute table that I have, I want the fields labelled ‘WAYPT’ and ‘COUNTY’.
Within the for loop, we’ll then have it write each row into the text file within f, and we’ll use ‘\n’ to get a new line after every loop.
To finish off the script—I want to remove the temporary shapefile we created for the intersect results. No reason to keep crowding my disk space. Since I know what’s that’s called, I can just use the os.remove to get rid of them. Keep in mind, shapefiles are actually a cluster of files. You have to delete all of them.
f.close os.remove(r"c\\Users\\Desktop\\Test.cpg") os.remove(r"c\\Users\\Desktop\\Test.dbf") os.remove(r"c\\Users\\Desktop\\Test.prj") os.remove(r"c\\Users\\Desktop\\Test.sbn") os.remove(r"c\\Users\\Desktop\\Test.shx") os.remove(r"c\\Users\\Desktop\\Test.shp") os.remove(r"c\\Users\\Desktop\\Test.shp.xml")
Once you have to script in a .py file (the copy/paste version is below), you’ll need to add it to a toolbox. Here’s a quick tutorial if you’ve never done that before. Just remember to set your two parameters to strings before you close out.
Before you run the script, you’ll need to have both the point and polygon layers loaded. Then when you run it, just plug in the names of the layers to run the analysis.
Hit Ok, and you’ll have a brand new text document telling you which points are within which counties. Like this:
Eventually I’m hoping to update this so that it exports specifically into a CSV format. A little prettier to deal with, but this at least gives you a snapshot of your intersect data. More ArcPy soon.
The Copy/Paste Version
import arcpy point = arcpy.GetParameterAsText(0) polygon = arcpy.GetParameterAsText(1) arcpy.Intersect_analysis([point, polygon],r"c:\\Users\\Desktop\\Test.shp") import os f = open(r"c:\\Users\\Desktop\\TestFile.txt","w") with arcpy.da.SearchCursor("c:\\Users\\Desktop\\Test.shp",['Field1', 'Field2']) for row in cursor: f.write(str(row) + '\n') f.close os.remove(r"c\\Users\\Desktop\\Test.cpg") os.remove(r"c\\Users\\Desktop\\Test.dbf") os.remove(r"c\\Users\\Desktop\\Test.prj") os.remove(r"c\\Users\\Desktop\\Test.sbn") os.remove(r"c\\Users\\Desktop\\Test.shx") os.remove(r"c\\Users\\Desktop\\Test.shp") os.remove(r"c\\Users\\Desktop\\Test.shp.xml")
I’m in a little dilemma, which I hope you can resolve.
Is it possible to create a Python script that sequentially processes all files in the folder it exists, using argparse?
I’m trying to implement Canny Edge Detection, with the biggest hurdle being processing input JPEG files.
LikeLike
AH, that sounds wildly advanced for me–I haven’t used argparse yet, but I’ll definitely look into it!
LikeLike
I’ve had the problem of the indentation before as well. It looks like you went with the pre tag, which is what I was going to suggest. I have noticed that there are time that WP is just being a b*tch and won’t work with you. Often, you just need to restart your browser and go back to editing your post.
LikeLike
Yeah, it’s a monumental pain in the ass reformatting every time I copy and paste from WordPress into a .PY file. So many errors. Glad to know others are frustrated by it too.
LikeLiked by 1 person
Just a suggestion, but maybe the approach needs to be to code and then put it into the tutorial. That’s what I do for mine. At least when the code is small. Otherwise I’ll just link a large code block to GitHub.
LikeLiked by 1 person