Running Viewshed on one feature at time in shapefile?

Question

I have over 100 features in a shapefile and I want to create individual viewsheds for each and then run processes on each of those resulting viewsheds in ArcGIS using Python.

Is it possible to write a script using a cursor to run through each "row" and use that individual observer point to run the viewshed?

Then update the cursor to the next row and rerun it until all features have been exhausted?

Or would I have to first Split the shapefile into over 100 files and run each individually?

Also, is it possible to designate the filenames for the resulting viewsheds (or shapefiles) using a specific field?

Answer 1

I've never used the viewshed tool so I can't speak to the specifics of using that tool, however in regards to your suggestion in this question I will add an answer here to build on KHibma's solution.

The first thing we need to address is the format of your SQL query: "[OID] = count". As it stands, this query will fail since "count", as represented here, is a string and strings need to be surrounded by single quotes: "[OID] = 'count'". Now, the second reason this query will fail is that OID is an integer field and as such there will be no OID value equal to "count". This is where string formatting comes into play: '"OID" = {}'.format(i) or '"OID" = {0}'.format(i) for those using Python 2.6.x (Use quotation marks to enclose the fieldname when using a shapefile. See ArcGIS help.).

To avoid certain pitfalls of making queries on the OBJECTID field, we can nest the loop into a cursor to make sure that we are only querying OID values that exist:

arcpy.MakeFeatureLayer_management ("C:/data/pts.shp", "pts")

with arcpy.da.SearchCursor('pts',['OID']) as cursor:
    for row in cursor:
        oid = row.OID # or whatever your OBJECTID field is called
        arcpy.SelectLayerByAttribute_management ("pts", "NEW_SELECTION", '"OID" = {}'.format(oid))
        outViewshed = Viewshed("elevation","pts",2,"CURVED_EARTH",0.15)
        outViewshed.save("C:/sapyexamples/output/outvwshd"+oid)

Answer 2

OK, borrowing from the answer from @Jason, I was able to tweak it a bit to finally come up with a script that works great. I also switched it to use FID instead of OID. Either way, it will work.

arcpy.MakeFeatureLayer_management (inPoints, "pts")

with arcpy.da.SearchCursor('pts',['FID']) as cursor:
    for row in cursor:
        fid = str(row[0]) ### Must be string
        print fid
        arcpy.SelectLayerByAttribute_management ("pts", "NEW_SELECTION", '"FID" = {}'.format(fid))
        outViewshed = Viewshed(inDEM,"pts",1,"CURVED_EARTH",0.15)
        outViewshed.save("C:/Python/output/viewsheds/"+fid)

Perhaps this will help someone in the future. Thanks everyone!

Answer 3

You could do a cursor, but its more straight forward to use the select by attributes tool. Simply put it in a loop ... Here's some untested code (you'll probably have to tinker with the expression in the select: (note I just add the count to the output name, if you want to get the specific value from a field then you'll have to open a cursor)

arcpy.MakeFeatureLayer_management ("C:/data/pts.shp", "pts")
count = int(arcpy.GetCount_management("pts").getOutput(0)) 
i=0
while i < count:
   arcpy.SelectLayerByAttribute_management ("pts", "NEW_SELECTION", "[OID] = i")
   outViewshed = Viewshed("elevation","pts",2,"CURVED_EARTH",0.15)
   outViewshed.save("C:/sapyexamples/output/outvwshd"+count)
   i+=1

Or you could do it all in a model. The "advanced" model at the bottom of this topic explains how to do it. Just cut out the tabulate area part, and do just the viewsheds with the iterator. Using the model only tool "get field value" you can use inline variables to update the final output name. This would actually be my preferred way to tackle this problem (model over script)