2

What I have:

  • 1/3 Arc Second DEM for the extent of my AOI
  • Feature class containing ~16 observers with fields for 'observer_offset,' 'outer_radius,' 'horizontal start angle,' and 'horizontal end angle'

What I am trying to accomplish:

  • Create an individual single-band geodesic viewshed for each individual observer in the feature class in tif format.

What I am attempting to use:

  • Viewshed2() in ArcGIS Pro 3.0.3

Code I am attempting to write:

import arcpy
from arcpy import env
from arcpy.sa import *


# set local variables
inRaster = r"phase_2\phase_2_dem"
inObservers = r"C:\Users\Austin\OneDrive\Documents\SynthosGIS\Geolinks\geolinks_fullrun\geolinks_fullrun\geolinks_fullrun.gdb\test_test_test"
var = None
analysisType = "FREQUENCY"
verticalError = "0 meters"
outAnalysisRelationTable = None
refractCoeff = 0.13
surfaceOffset = "0 meters"
observerElevation = None
observerOffset = "observer_off"
innerRadius = None
innerIs3D = "GROUND"
outerRadius = "Tower_Range_"
outerIs3D = "GROUND"
horizStartAngle = "hor_start_angle"
horizEndAngle = "hor_end_angle"
vertUpperAngle = 90
vertLowerAngle = -90
analysisMethod = "ALL_SIGHTLINES"

with arcpy.da.SearchCursor(inObservers, ['Account_Name']) as vwshd_itterator:
    for row in vwshd_itterator:
        account = str(row[0])
        query = f"Account_Name = '{account}'"
        arcpy.SelectLayerByAttribute_management(inObservers, "NEW_SELECTION", query)
            
        out_raster = arcpy.sa.Viewshed2(inRaster, inObservers, var, analysisType,
                         verticalError, outAnalysisRelationTable, refractCoeff,
                         surfaceOffset, observerElevation, observerOffset,
                         innerRadius, innerIs3D, outerRadius, outerIs3D,
                         horizStartAngle, horizEndAngle, vertUpperAngle,
                         vertLowerAngle, analysisMethod)
        out_raster.save(r"C:\Users\Austin\OneDrive\Documents\SynthosGIS\Geolinks\geolinks_fullrun\geolinks_fullrun\vwshd_" + account + ".tif")

What I Need:

  • This is the first time that I have used both the Viewshed2 script and the DA SearchCursor function. I am not sure if I have integrated the two functions correctly, or what I will need to do make it run. As I type, I won't be able to try out any scripts for a couple days as I am away from my workstation. Any input to this script, so that I can hit the ground running when I get back?

Update 1: RESOLVED

  • I have updated the code to what is currently showing. I am currently running into an issue with the query for the selection in the for loop. Am I implementing the f-string incorrectly?

Update 2: UNRESOLVED

  • I have achieved a function that will run. However, the select by attribute function does not appear to be looping through features individually. Instead it selects all features on the layer at once and is creating a separately named file reprocessing the same viewshed over and over again.
16
  • 1
    You don't seem to be doing anything with input_feature_class inside the loop at all.
    – Vince
    Commented Mar 24, 2023 at 18:02
  • @Vince, I edited to use the appropriate variable. I had been working between a couple example scripts to learn the syntax in this and completely just missed that I was using the wrong variable. Commented Mar 24, 2023 at 18:23
  • 1
    A DA Cursor doesn't change the selection environment, and altering the selection environment in the middle of cursor use is fraught, at best. Instead, grab the OIDs as a cursor-generated list outside the loop, then use SelectByAttribute using OID inside the loop.
    – Vince
    Commented Mar 24, 2023 at 18:45
  • So, I would use the search cursor to to .append() the OID to a new empty list, then use something like this? ''' for s in new_list: arcpy.management.SelectLayerByAttribute('inObservers', 'NEW_SELECTION', where 'OID' == s Viewshed2() ''' Commented Mar 24, 2023 at 19:41
  • Well, something like, oidList = [r[0] for arcpy.da.SearchCursor(inObservers, ['OID@'])]... The OID column is different by source datatype, so it might be """"FID" = {:d}""".format(s)
    – Vince
    Commented Mar 24, 2023 at 21:02

1 Answer 1

1

After a great deal of back and forth with Vince on throwing this together, I am happy to say I have reached a functional final result. The below successfully translates my feature class into a functional feature layer, the iterates over that layer and stores an attribute field into a list to use as selection criteria to loop the viewshed2 tool through on each individual feature.

import arcpy
from arcpy import env
from arcpy.sa import *


# set local variables
inRaster = r"phase_2\phase_2_dem"
inObservers = r"C:\Users\Austin\OneDrive\Documents\SynthosGIS\Geolinks\geolinks_fullrun\geolinks_fullrun\geolinks_fullrun.gdb\fc_phase_2_towers"
var = None
analysisType = "FREQUENCY"
verticalError = "0 meters"
outAnalysisRelationTable = None
refractCoeff = 0.13
surfaceOffset = "0 meters"
observerElevation = None
observerOffset = "observer_off"
innerRadius = None
innerIs3D = "GROUND"
outerRadius = "Tower_Range_"
outerIs3D = "GROUND"
horizStartAngle = "hor_start_angle"
horizEndAngle = "hor_end_angle"
vertUpperAngle = 90
vertLowerAngle = -90
analysisMethod = "ALL_SIGHTLINES"
inLayer = 'pts'

arcpy.management.MakeFeatureLayer(inObservers, inLayer)

account = [str(row[0]) for row in arcpy.da.SearchCursor(inLayer, ['Account_Name'])]
print(account)

for a in account:
    query = f"Account_Name = '{a}'"
    arcpy.SelectLayerByAttribute_management(inLayer, "NEW_SELECTION", query)
    out_raster = Viewshed2(inRaster, inLayer, var, analysisType,
                         verticalError, outAnalysisRelationTable, refractCoeff,
                         surfaceOffset, observerElevation, observerOffset,
                         innerRadius, innerIs3D, outerRadius, outerIs3D,
                         horizStartAngle, horizEndAngle, vertUpperAngle,
                         vertLowerAngle, analysisMethod)

    out_raster.save(r"C:\Users\Austin\OneDrive\Documents\SynthosGIS\Geolinks\geolinks_fullrun\test_itterate\vwshd_" + str(a) + ".tif")
    print('completed ' + str(a))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.