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()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")desc = arcpy.Describe(fcpoints)
if desc.hasOID:
fldOID = desc.OIDFieldName
else:
exit("featureclass does not have an Object ID field") fld = arcpy.AddFieldDelimiters(fcpoints,fldOID)
whereclause = "{1} = {0}".format(oid,fld)import arcpy, os, numpy
from numpy import *
# settings of your data
fcpoints = r'G:\xxx\xxxx\xxxx\xxxx\xxx\xxxx\xxx\clip16points.shp'
fldVisible = 'VisiblePix' # will contain the number of visible pixels
dem = r'G:\xxx\xxxx\xxxx\xxxx\xxx\xxxx\xxx\combraster2'
# ouput workspace and name for viewshed
outputWS = r'C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.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")
# get the oid field
desc = arcpy.Describe(fcpoints)
if desc.hasOID:
fldOID = desc.OIDFieldName
else:
exit("featureclass does not have an Object ID field")
# 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 a whereclause
fld = arcpy.AddFieldDelimiters(fcpoints,fldOID)
whereclause = "{1} = {0}".format(oid,fld)
# create featurelayer for current observer point
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))
VisiblePix = myArray2.sum()
print "There are {0} visible pixels for point with oid={1}".format(VisiblePix,oid)
# update value in point featureclass
row[1] = VisiblePix
cursor.updateRow(row)
del myArray, myArray2
del row
# release license
arcpy.CheckInExtension("Spatial")>>> import arcpy, os, numpy
... from numpy import *
...
... # settings of your data
... fcpoints = r'C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb\clip16points'
... fldVisible = 'VisiblePix' # will contain the number of visible pixels
... dem = r'C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb\clipdem30mWGS1'
...
... # ouput workspace and name for viewshed
... outputWS = r'C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.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")
... desc = arcpy.Describe(fcpoints)
... if desc.hasOID:
... fldOID = desc.OIDFieldName
... else:
... exit("featureclass does not have an Object ID field")
...
... # 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 = "fldOID = {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")