Hello I have created this code to run a number of viewsheds for a group of points shp Pnt1, Pnt2... and a group of Dems, dem1, dem2 etc. For some reason I don't have any results but I don't know why. This is the code:
import arcpy, os from arcpy import env from arcpy.sa import * arcpy.env.overwriteOutput = True env.workspace = arcpy.GetParameterAsText(0) out = arcpy.GetParameterAsText(1) fc = arcpy.ListFeatureClasses("Pnt*", "Point") ras = arcpy.ListRasters("dem*", "GRID") point = "Pntclip_pol1" dem = "dempol1" i = 1 for shp in fc: for raster in ras: if (shp == point and raster == dem): inRaster = raster inObserverFeatures = shp outViewshed = Viewshed(inRaster, inObserverFeatures, "") outViewshed.save(out + "\\" + "view" + str(i)) i = int(i) + 1 point = "Pntclip_pol" + str(i) dem = "clippol" + str(i)
Solved! Go to Solution.
Your nested if statements seems to be the problem. In such cases, you would be advised to make a list containing just the pairs that you want to process. You formulate the pairs
[ [ A with rastA],
[ B with rastA],
[B with rastB]
etc ]
then run throught list of pairs. It is the difference between combinations and permutations? which do you want?
Ok I developed the code at that form but I have an error of ERROR 010004: All cells in Raster clippol1 have the NODATA value. Stop execution. ERROR 010067: Error in executing grid expression. Failed to execute (Viewshed).
This is the code where I tried to define extent also but I get back only one viewshed and not with the coding I want (line 11). And these errrors keep appear, I extracted these rasters from Modelbuilder I haven't set anything to NODATA so as to have problem
import arcpy, os from arcpy import env from arcpy.sa import * ink = "C:/DEFENCE_DATA/Results/ArcMap/Visibility/Calc_View" fc = arcpy.ListFeatureClasses("Pnt*", "Point") demi = "clippol1" i = 1 for seip in fc: arcpy.env.extent = ink + "/clippol" + str(i) outViewshed = Viewshed(demi, seip) outViewshed.save(ink + "/view" + str(i)) i = i + 1 dem = "clippol" + str(i)
I think that I am very close to finally make it run but I have to solve that with NODATA
You had DEMi and DEM below, as written it would always just use the first dem, your next dem is not used..
So your second viewshed may not have a point1 in it.
If your running this from the IDLE or pythonwin, I’d use print commands to see what your variables are, so if it errors off you have a better clue.
I wrote what I’d use in blue below. I also shortened your output string.
One big potential problem is if your typing in the workspace path, it’s so easy to confuse the forward and back slash.
Are the points and dems in the same projection?
Since I had some time, I rewrote the code, see below in green. It works, my e-mail is William.Chappell@its.ny.gov<mailto:William.Chappell@its.ny.gov> if you want the script I used, in case copy/paste messes with the indent formats.
Good Luck, Bill
I tried it Mr. Chappell but I only get back two viewsheds and then I have again the same error
Walk through the process manually, is do it with out the scripts, see if you get an error manually on the third one..
It may be the data or the name of the third one.
Bill Chappell
ITS Group - GIS Developer
800 North Pearl Street- Room 222
Menards, NY 12204
518-408-0185
William.Chappell@Its.ny.gov
Well manually it works perfectly but with the script I have no idea why not. I did a search and I saw that that error is usual and it doesn't have a logic in many circumstances. I have sent you 5 pairs of rasters and shp to try it if you have the time, because I don't think I am doing something wrong, this error is very strange and I believe propably has to do with the software.
Problem solved
import arcpy, os from arcpy import env from arcpy.sa import * arcpy.env.overwriteOutput = True arcpy.CheckOutExtension("Spatial") #arcpy.env.workspace = arcpy.GetParameterAsText(0) arcpy.env.workspace = "C:/Results/ArcMap/Visibility/Calc_View" ink = arcpy.env.workspace fc = arcpy.ListFeatureClasses("Pn*", "Point") ras = arcpy.ListRasters("cl*", "GRID") i = 0 for seip in fc: i= i + 1 demi = "clpol" + str(i) vect = "Pnclpol" + str(i) + ".shp" for rast in ras: if (seip == vect and rast == demi): outViewshed = Viewshed(rast, seip) outViewshed.save("C:/Results/ArcMap/Visibility/Calc_View/View" + str(i))
Error 010004 is very hard to get away and I confirm as in the post Raster NoDATA Value Error 010004 Driving me crazy! that renaming your files may be a solution to the problem also.
This ran for me using your data..
import arcpy, os
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
arcpy.env.workspace = "C:/test"
ink = "C:/test"
fcList = arcpy.ListFeatureClasses("Pnt*", "Point")
ras = "clippol"
i = 1
for fc in fcList:
ras = ras + str(i)
print "Point is: " + fc + " the DEM is: " + ras
outViewshed = Viewshed(ras, fc)
outViewshed.save("C:/test/out/outvwshd"+ str(i))
ras = "clippol"
i=i+1
print "Done"