Hi Anne,I notice a few more things. I wrote the code based on rasters and points FC in a file geodatabase. You are using grids and a point shapefile. One of the errors is "No valid Observation points present.". This is likely the cause of the wrong whereclause that is being build. The whereclause is using the field OBJECTID (standard in fgdb fc's) and the field is not delimited. To make this a little more gentle the code should be changed a bit:1) use the describe methode to obtain the oid field (if it has one, if not exit the procedure):desc = arcpy.Describe(fcpoints)
if desc.hasOID:
fldOID = desc.OIDFieldName
else:
exit("featureclass does not have an Object ID field")
2) the whereclause should be build differently: fld = arcpy.AddFieldDelimiters(fcpoints,fldOID)
whereclause = "{1} = {0}".format(oid,fld)
This will result probably in a whereclause like: "FID" = 0So the final code in your case would be something like: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")
Before you run it please check you input dem and environment settings as explained in my previous post.It this still doesn't work, try to run the code with a local copy of your dem. The 'G' drive most likely maps to a network folder and I have seen issues with running scripts with grids on network locations that failed, when running it with a local copy actually worked.Kind regards,Xander