Adding total visible pixels for each observer location from created viewsheds

AnneHarrison

Hi,

I have used iteration in ModelBuilder to create multiple viewsheds for each observer location in a point feature class. However, I now need to transfer the value for visible pixels back to the observer feature class, so that I can then compare the total visible areas of different observer points. I can't join the tables as there is no common field in the raster and points files.

Does anyone know a way of automating this without manually having to insert the values into a new field in the points attributes table?

Thanks!

XanderBakker

Hi Anne,

I think that you can use python for this. You would probably have to convert the model to python and clean it up a bit.

I suppose you loop through the observer points and for each point you are calculating a viewshed raster. This raster contains 1 for visible and 0 for not visible. Using numpy it is easy to get the sum of the raster. In case you have NoData cells, the numpy array needs to be correct for this. The code is very simple:

import numpy, arcpy
from numpy import *
myArray = arcpy.RasterToNumPyArray("Viewshed1")
# avoid erroneous sum due to NoData cells
myArray2 = where(myArray==1,ones(myArray.shape),zeros(myArray.shape))
visibleCells = myArray2.sum()

In the output field of the point featureclass (observers), you can use a simple updatecursor (and a whereclause) to set the calculated value.

The script could look something like this:

import arcpy, os, numpy
from numpy import *

# settings of your data
fcpoints = r'C:\Project\_Forums\Viewshed\test.gdb\sources'
fldVisible = 'VisiblePix' # will contain the number of visible pixels
dem = r'C:\Project\_Forums\cost\fgdb\test.gdb\elevutm16n'

# ouput workspace and name for viewshed
outputWS = r'C:\Project\_Forums\Viewshed\test.gdb'
prefixVS = 'Viewshed'
arcpy.env.workspace = outputWS
arcpy.env.overwriteOutput = True

# grab a SA license
if arcpy.CheckExtension("Spatial") == "Available":
    arcpy.CheckOutExtension("Spatial")
else:
    exit("no SA license available")

# loop through points
fields = ("OID@",fldVisible)
with arcpy.da.UpdateCursor(fcpoints, fields) as cursor:
    for row in cursor:
        oid = row[0]

        # define output name (of oid = 1, raster will be called Viewshed1
        outVSname = outputWS + os.sep + prefixVS + str(oid)

        # create featurelayer for current observer point
        whereclause = "OBJECTID = {0}".format(oid)
        FLsources = "sources_Layer"
        arcpy.MakeFeatureLayer_management(fcpoints,FLsources,whereclause,outputWS)

        # calculate viewshed
        arcpy.gp.Viewshed_sa(dem,FLsources,outVSname,"1","FLAT_EARTH","0,13")

        # determine number of visible cells
        myArray = arcpy.RasterToNumPyArray(outVSname)
        myArray2 = where(myArray==1,ones(myArray.shape),zeros(myArray.shape))
        visPixels = myArray2.sum()
        print "There are {0} visible pixels for point with OID={1}".format(visPixels,oid)

        # update value in point featureclass
        row[1] = visPixels
        cursor.updateRow(row)
        del myArray, myArray2

del row

# release license
arcpy.CheckInExtension("Spatial")

If you have any questions, please let me know.

Kind regards,

Xander