1

I need a python script that seeks to cycle through a number of viewpoints, deriving an AGL raster in each case (using arcpy.Viewshed_3d) and progressively summing the AGL rasters into a single cumulative raster. I have written a Python script to do this but am hitting a problem with the chunk of code that iteratively combines the individual AGL rasters into a single output.

The script takes as input a vector viewpoint layer and a DEM and it seems to work insofar as it loops as it should and generates a final output. However comparison with a manual 'testbed' dataset generated for a small sample of points reveals that the resultant values do not match. Trials with single points worked fine (i.e. for a single point the resultant AGL layer generated by the script was identical to one manually generated) but in the case of two points the output seemed to be the sum of the 2nd AGL raster twice (rather than AGL1 + AGL2). I created a debug version of the code using elif for a small number of points (and saving the temporary raster each time with a unique filename) and it worked fine - the AGL layers generated at each iteration matched the manual versions as did the final AGL output. However, this is not really practical for 100K+ viewpoints.

The faulty chunk of Code is given below:

# set up cursor to step through viewpoint file one at a time
searcher = arcpy.da.SearchCursor(inputVector, ["OID@"])

try:
        for row in searcher:
            # create layer
            arcpy.MakeFeatureLayer_management(inputVector, "viewLayer")
            # set query
            query = ' "OBJECTID" = ' + str(row[0])
            #select point
            arcpy.SelectLayerByAttribute_management("viewLayer", "NEW_SELECTION", query)
            # set viewshed parameters
            inRaster = inputDem
            outViewshed = "outView"
            outagl = "hiding"
            # calculate viewshed
            arcpy.Viewshed_3d(inRaster, "viewLayer", outViewshed, 1, "CURVED_EARTH", 0.13, "outAgl")
            # convert nodata values to zero in agl output to allow union of inputs when adding rasters
            hide1 = Raster("outAgl")
            hidecon = Con(IsNull(hide1), 0, hide1)
            arcpy.AddMessage(query) #debug line

            #create/update summed agl raster
            if row[0] == 1:
                arcpy.AddMessage("first time through loop") #debug line
                tempraster = hidecon #create the raster that will be used to sum 

            else:
                arcpy.AddMessage("looping") # debug line
                tempraster = tempraster + hidecon #add latest AGL to the summed raster

        # save the final output summed raster
        tempraster.save(outputAgl)

        #clean up any non-temporary working rasters
        arcpy.Delete_management("outAgl")
        arcpy.Delete_management("outView")

        #clean up cursor
        del searcher    

Apologies - I am a very amateur python coder.

UPDATE: Below is a debug version (once again - apologies for untidy code) that replaced the loop with else and elif and saved the intermediate AGL outputs. This worked insofar as the final AGL (outputAGL) was exactly what I was after and matched the test dataset (but it is not a practical solution given the number of AGLs I need to sum (1000+)

# import system., OS and arcpy modules
import sys, os, arcpy

# import geoprocessing environments
from arcpy import env

# import Spatial Analyst and Data Access modules
from arcpy.sa import *
from arcpy.da import *

#check out 3D analyst extension
arcpy.CheckOutExtension("3D")

# enable file overwrite & ensure output extent is max. of inputs
arcpy.env.overwriteOutput = True
arcpy.env.extent = "MAXOF"

# get input files and output filename
inputVector = arcpy.GetParameterAsText(0)
inputDem = arcpy.GetParameterAsText(1)
outputAgl = arcpy.GetParameterAsText(2)

# set up cursor to step through viewpoint file one at a time
searcher = arcpy.da.SearchCursor(inputVector, ["OID@"])

try:
    for row in searcher:
        # create layer
        arcpy.MakeFeatureLayer_management(inputVector, "viewLayer")
        # set query
        query = ' "OID" = ' + str(row[0])
        #select point
        arcpy.SelectLayerByAttribute_management("viewLayer", "NEW_SELECTION", query)
        # set viewshed parameters
        inRaster = inputDem
        outViewshed = "outView"
        outagl = "hiding"
        # calculate viewshed
        arcpy.Viewshed_3d(inRaster, "viewLayer", outViewshed, 1, "CURVED_EARTH", 0.13, "outAgl")
        # convert nodata values to zero in agl output to allow union of inputs when adding rasters
        hide1 = Raster("outAgl")
        hidecon = Con(IsNull(hide1), 0, hide1)
        # debug line
        arcpy.AddMessage(query)

        #create/update summed agl raster
        if row[0] == 1:
            # debug line
            arcpy.AddMessage("first time through loop")
            tempraster = hidecon
            tempraster.save("D:/python_2016/loop1")

        # debug chunk
        elif row[0] == 2:
            arcpy.AddMessage("second time..")
            # tempHelp = tempraster
            tempraster = tempraster + hidecon
            tempraster.save("D:/python_2016/loop2")

        elif row[0] == 3:
            arcpy.AddMessage("third time...")
            tempraster = tempraster + hidecon
            tempraster.save("D:/python_2016/loop3")

        elif row[0] == 4:
            arcpy.AddMessage("fourth time...")
            tempraster = tempraster + hidecon
            tempraster.save("D:\python_2016\loop4")
        # debug chunk end

        else:
            # debug line
            arcpy.AddMessage("looping")
            tempHelp = tempraster
            tempraster = tempHelp + hidecon
            ## saver = "delete_me"
            ## tempraster.save(saver)

    # save the final output summed raster
    tempraster.save(outputAgl)

    #clean up any non-temporary working rasters
    arcpy.Delete_management("outAgl")
    arcpy.Delete_management("outView")
    ## arcpy.Delete_management(saver)

    #clean up cursor
    del searcher

except:
#issue error message
arcpy.AddMessage("The Script crashed and burned sadly")
# report any errior message generated
arcpy.AddMessage(arcpy.GetMessages())
1

Finally resolved the problem this morning by adding tempraster.save() to persist the accumulating raster between loop iterations.

            #create/update summed agl raster
        if row[0] == 1:
            tempraster = hidecon #create the raster that will be used to sum
            tempraster.save()

        else:
            tempraster = tempraster + hidecon #add latest AGL to the summed raster
            tempraster.save()

    # save the final output summed raster
    tempraster.save(outputAgl)

The resulting code worked like a dream. Many thanks to Felix for helping me to focus in on exactly where the problem was - your help was very much appreciated, Mark

0

Hard to tell what is wrong, because a code is a bit untidy. So I modified it slightly:

try:  arcpy
except:  import arcpy
import traceback, os, sys
from arcpy import env
env.overwriteOutput = True
from arcpy.sa import *

inputVector=r'D:\Scratch\Walkover.gdb\rpoints'
inRaster = r'D:\Aerials\Backups\arc2mdem'
outAgl="outAgl"
outViewshed="outViewshed"
outagl = "hiding"

try:
  def showPyMessage():
    arcpy.AddMessage(str(time.ctime()) + " - " + message)
  searcher=[row[0] for row in arcpy.da.TableToNumPyArray(inputVector,"OID@")]
  for row in searcher:
    query = ' "OID" = ' + str(row)
    arcpy.MakeFeatureLayer_management(inputVector, "viewLayer",query)
    arcpy.Viewshed_3d(inRaster, "viewLayer", outViewshed, 1, "CURVED_EARTH", 0.13, outAgl)
    hidecon = Con(IsNull(outAgl), 0, outAgl)
    if row==1:
      tempraster = hidecon
    else:
      tempraster = tempraster + hidecon
  tempraster.save(r'D:\Aerials\nztm\toriver')

except:
  message = "\n*** PYTHON ERRORS *** "; showPyMessage()
  message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
  message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()            

Test results:

enter image description here

I struggle to understand why you are totalling AGLs

  • Many Thanks Felix - I am just back in to work so will give it a try. Mark – Mark Gillings Jun 24 '16 at 6:52
  • I've run a new series of tests against a test dataset (manually calculated set of AGLs for 4 widely spaced viewpoints that were then summed using map algebra). Felix's tidy code generates the same result as my original with neither matching the test dataset. Instead of generating an AGL at each iteration and progressively adding it to a 'sum' raster, both seem to be multiplying the final AGL generated (i.e. for 4 viewpoints AGL#4) by the number of times it loops (i.e. viewpoints). So in my test the code generates a final product that is 4 x AGL#4 (instead of AGL#1 + AGL#2 + AGL#3 + AGL#4). – Mark Gillings Jun 24 '16 at 11:46
  • I have updated the original question to include a debug version of the code I produced in an attempt to resolve the issue. This uses elif and saves the interim AGL for each point. In this case the final summed output matched the test dataset exactly (i.e. succeeded in summing AGL#1+AGL#2+AGL#3+AGL#4 into a single raster). – Mark Gillings Jun 24 '16 at 11:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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