Problem with If , For statements

[KONPETROV]

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)

[DanPatterson_Retired]

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?

[KONPETROV]

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

[BillChappell]

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

[KONPETROV]

I tried it Mr. Chappell but I only get back two viewsheds and then I have again the same error

[BillChappell]

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

[KONPETROV]

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.

[KONPETROV · Accepted solution]

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.

[BillChappell]

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"