topic Re: Stand alone script fails using spatial analyst tools in Python Questions

Stand alone script fails using spatial analyst tools

ElizabethLewis — Mon, 05 Mar 2012 08:18:38 GMT

Hi,

I'm trying to create a stand alone script to automate watershed analysis. So far it works up until it gets to the spatial analyst
tools but then fails. Can anyone show me where I have gone wrong please? Here is my code:

import arcpy
import os
import time
import sys
from arcpy import env
from arcpy.sa import *

env.workspace = "H:/PhD/PythonLessons/Scripts/testIt"
catchmentName = "CatchmentName"
catchmentRaster = catchmentName + ".tif"
catchmentRasterProj = catchmentName + "proj.tif"
catchment100 = catchmentName + "100.tif"
catchment100Fill = catchmentName + "100Fill.tif"
catchment100Flow = catchmentName + "100Flow.tif"
catchment100Acc = catchmentName + "100Acc.tif"
catchment100PP = catchmentName + "100PP.tif"
catchment100WS = catchmentName + "100WS.tif"
folder = "H:/PhD/PythonLessons/Scripts/testIt"
folderAdd = "H:/PhD/PythonLessons/Scripts/testIt/"
outlet = "outlet.shp"
projection = "C:/Program Files (x86)/ArcGIS/Desktop10.0/Coordinate Systems/Projected Coordinate Systems/National Grids/Europe/British National Grid.prj"
# identify guage
#create pourpoint shapefile
# identify catchment outline
# select appropriate testIt

# have a folder with all of the little ascii testIt for a catchment. This folder will be called *catchment name*
print "checking extensions"

try:
 arcpy.CheckExtension("Spatial") == "Available"
 arcpy.CheckOutExtension("Spatial")
 print "done"
except:
 print "piddle"
 

print "converting the ascii files to raster files"

asciiList = os.listdir(folder)    

for x in range(len(asciiList)-1, -1, -1):
 if not asciiList.endswith(".asc"):
   asciiList.pop(x)


for file in asciiList:  
 print folderAdd + file
 print folderAdd + file[:-4]
 
 try:
  arcpy.ASCIIToRaster_conversion(folderAdd+ file, folderAdd+ file[:-4] + ".tif", "INTEGER")
 except:
  print "Processing " + file + "FAILED"
  
 time.sleep(3)
 
tifList = os.listdir(folder)

print tifList

for x in range(len(tifList) -1, -1, -1):
    if not tifList.endswith(".tif"):
        tifList.pop(x)

print tifList

tifString = ""
for tif in tifList:
 tifString += tif + ";"
 
tifString = tifString[:-1]

print tifString

print "mosaic to new raster"

try:
 arcpy.MosaicToNewRaster_management(tifString, folderAdd, catchmentRaster, "#", "16_BIT_SIGNED", "#", "1" , "MEAN")
except:
 print "Error with mosaic to new raster"
 
time.sleep(3)


print "Giving it the right coordinates"
try:
 arcpy.DefineProjection_management(folderAdd + catchmentRaster, projection)
 arcpy.ProjectRaster_management(folderAdd + catchmentRaster, folderAdd + catchmentRasterProj, projection, "BILINEAR", "#", "#","#",projection)
except: 
 print "Error giving it the right coordinates"

time.sleep(3)
 
print "Changing resolution"

try:
 outAggreg = Aggregate(catchmentRasterProj, 10, "MEDIAN", "TRUNCATE", "DATA")
 outAggreg.save(folderAdd + catchment100)
except:
 print "Error changing resolution"
 
time.sleep(3)
 
print "Filling sinks"

try:
 outFill = Fill(catchment100)
 outFill.save(folderAdd + catchment100Fill)
except:
 print "error filling sinks"
 
print "Calculating fow direction"

time.sleep(3)

try:
 outFlowDirection = FlowDirection(catchment100Fill, "NORMAL")
 outFlowDirection.save(folderAdd + catchment100Flow)
except:
 print "Error calculating flow direction"
 
 
print "Calculating flow accumulation"

time.sleep(3)

try:
 outFlowAccumulation = FlowAccumulation(catchment100Flow)
 outFlowAccumulation.save(folderAdd + catchment100Acc)
except:
 print "Error calculating flow accumulation"

print "Creating pour point snap"

time.sleep(3)

try:
 outSnapPour = SnapPourPoint(outlet, catchment100Acc, 5, "FID")
 outSnapPour.save(folderAdd + catchment100PP)
except:
 print "Error snapping pour point"
 
time.sleep(3)

print "Creating watershed"

try:
 outWatershed = Watershed(catchment100Flow, catchment100PP)
 outWatershed.save(folderAdd + catchment100WS)
except:
 print "Error creating watershed"
 
time.sleep(3)
 
print "Creating files for SHETRAN prepare"

try:
 arcpy.RasterToASCII_conversion(catchment100, folderAdd + catchmentName + "100DEM.txt")
 arcpy.RasterToASCII_conversion(catchment100WS, folderAdd + catchmentName + "100Mask.txt")
except:
 print "error creating files for SHETRAN prepare"
 
time.sleep(3)
 
print "Done!"

raw_input("close to exit")

Thanks very much

Re: Stand alone script fails using spatial analyst tools

DanPatterson_Retired — Mon, 05 Mar 2012 10:35:48 GMT

have you tried appending _sa to all of the spatial analyst tools as you have done with the _management etc tools?

Re: Stand alone script fails using spatial analyst tools

FabianBlau — Mon, 05 Mar 2012 10:51:38 GMT

Please post the error message.
Appending _sa should not help with your way to import the functions, but try it anyway. Try also to call arcpy.Fill_sa(...) etc.. I think this would be also better readable.

Re: Stand alone script fails using spatial analyst tools

ElizabethLewis — Mon, 05 Mar 2012 12:28:35 GMT

Hi,

I've tried arcpy.Aggregate_sa(...) but it still doesn't work. I don't get an error message apart from my own:

"Giving it the right coordinates"
"Changing resolution"
"Error changing resolution"
"error filling sinks"
"Calculating fow direction"
"Error calculating flow direction"
"Calculating flow accumulation"
"Error calculating flow accumulation"
"Creating pour point snap"
"Error snapping pour point"
"Creating watershed"
"Error creating watershed"
"Creating files for SHETRAN prepare"
"error creating files for SHETRAN prepare"
"Done!"

Any other ideas?

Thanks

Re: Stand alone script fails using spatial analyst tools

FabianBlau — Mon, 05 Mar 2012 13:05:36 GMT

It would be better to print the arcgis-errormessages. Dont use try/except in this way. When any error occurs your script dont stop. The following code will crash also (in your case).
Do you have an arcinfo-license? aggregate needs it.

Re: Stand alone script fails using spatial analyst tools

ElizabethLewis — Sat, 11 Dec 2021 19:23:12 GMT

I've changed the script so that it now looks like this:

try:
 outAggreg = arcpy.Aggregate_sa(catchmentRasterProj, 10, "MEDIAN", "TRUNCATE", "DATA")
 arcpy.outAggreg.save_sa(folderAdd + catchment100)
except Exception as e:
    print e.message

The error message is: 'module' object has no attribute 'Aggregate_sa@

Thanks

Re: Stand alone script fails using spatial analyst tools

FabianBlau — Tue, 06 Mar 2012 14:08:22 GMT

I tried aggregate in my python window.

This didnt work (also the AttributeError :mad:):
arcpy.Aggregate_sa('ws', 'Aggrega_ws2')

This was ok:
arcpy.sa.Aggregate("ws",5,"MAXIMUM","EXPAND","DATA")

Re: Stand alone script fails using spatial analyst tools

ElizabethLewis — Tue, 06 Mar 2012 14:17:52 GMT

arcpy.sa.Aggregate(...) gave this error message:

ERROR 999999: Error executing function.
failed to open raster dataset
Failed to execute <aggregate>.

Also, I'm trying to do this as a stand alone script and not in the python window.

Thanks!

Note: Tool_sa syntax invalid at 10.x

curtvprice — Sat, 11 Dec 2021 19:23:15 GMT

have you tried appending _sa to all of the spatial analyst tools as you have done with the _management etc tools?

Dan, this changed at 10.x. The Spatial Analyst licensed functions need to be access the new python-esque way:

arcpy.sa.Con(...)
from arcpy.sa import *
out = Con(..)

If you have a 3D license, you can access the raster tools included in 3D using the _3d suffix.

Re: Stand alone script fails using spatial analyst tools

curtvprice — Sat, 11 Dec 2021 19:23:18 GMT

I'm trying to do this as a stand alone script and not in the python window.

Should not matter.

My guess is your raster is getting saved somewhere (scratch workspace) besides where you think it is. Your script seems to imply that your processing results are expected in "H:/PhD/PythonLessons/Scripts/testIt" but you haven't set the scratch workspace to that location.

For best results when using Spatial Analyst tools its best when they are sent to the same folder:

from arcpy.env import ENV
ENV.workspace = r"E:\work"
ENV.scratchWorkspace = ENV.workspace

Setting the current and scratch to the same folder allows temporary grids to be renamed instead of copied - this can save a lot of time with big rasters.

Re: Stand alone script fails using spatial analyst tools

ElizabethLewis — Sat, 11 Dec 2021 19:23:20 GMT

Ok, I set

env.scratchWorkspace = env.workspace

and tried this set up

try:
        outAggreg = Aggregate(catchmentRasterProj, 10, "MEDIAN", "TRUNCATE", "DATA")
 outAggreg.save(catchment100)
except Exception as e:
    print e.message

for my spatial analyst tools.

It now gives the error 'failed to open raster dataset'. I'm not sure why it won't open it. I checked that it would open in ArcMap 10 and it does that just fine. I've also checked that the file name is right. I checked the license and that I'd imported the right modules.

Also, when I run the, as well as an error it creates a new folder within the workspace called 'aggrega_catc1' but I thought that setting the scratch workspace to the env workspace would stop this?

Thanks again

Re: Stand alone script fails using spatial analyst tools

curtvprice — Sat, 11 Dec 2021 19:23:23 GMT

It now gives the error 'failed to open raster dataset'. I'm not sure why it won't open it.

To see if it's a path issue, make sure it exists. (If you don't specify a full path, it will look in env.workspace for it.)

if not arcpy.exists(catchmentRasterProj): 
    raise Exception, "%s does not exist" % catchmentRasterProj

In the same manner, the output is also written to the current workspace, so instead of

outFill.save(folderAdd + catchment100Fill)

you can just specify

outFill.save(catchment100Fill)

as well as an error it creates a new folder within the workspace called 'aggrega_catc1' but I thought that setting the scratch workspace to the env workspace would stop this?

This is a temporary raster; all tools write a temporary raster to the scratch workspace when raster tools run. As I said the point of setting the workspaces the same is so when the tool successfully runs this temp raster is renamed (instead of copied) to your output raster when you execute outraster.save.

Re: Stand alone script fails using spatial analyst tools

ElizabethLewis — Thu, 08 Mar 2012 07:55:19 GMT

I tried arcpy.Exists and it can find the file so it's not my file path that's wrong. Aggregate still won't open the file raising the error 'failed to open raster dataset'.

Thanks for explaining the temporary files.

Re: Stand alone script fails using spatial analyst tools

DorothyMortenson — Thu, 15 Mar 2012 20:12:57 GMT

I had been trying to do something similar with the SnapPointPoint. I am converting from the gp to arcpy.
What I found was I was trying to create a grid and it didn't like it.
So I created a filegeodatabase for my results to go into. This seems to have solved my issue.
All of my original data is in GRID from arcgis 9.3 I plan to convert it, but I wanted to get the python scripts working first.

Re: Stand alone script fails using spatial analyst tools

curtvprice — Thu, 15 Mar 2012 20:46:59 GMT

I had been trying to do something similar with the SnapPointPoint. I am converting from the gp to arcpy.
What I found was I was trying to create a grid and it didn't like it.
So I created a filegeodatabase for my results to go into. This seems to have solved my issue.
All of my original data is in GRID from arcgis 9.3 I plan to convert it, but I wanted to get the python scripts working first.

The grid data format is very path sensitive - you need to be very careful to avoid long grid names or names that do not follow the very restrictive grid/coverage naming rules. On the plus side, the grid data format is usually faster to process (though this depends on the tool).

Re: Stand alone script fails using spatial analyst tools

ElizabethLewis — Sat, 11 Dec 2021 19:23:26 GMT

Hi,

I've finally resolved this issue. It turns out that the arc modules wouldn't run on the H drive that I was using so I saved everything to the C drive instead and it all ran fine.

Here's the code for future reference:

import arcpy
import os
import time
import sys
from arcpy import env
from arcpy.sa import *

env.workspace = "C:/temp/testIt"
env.scratchWorkspace = env.workspace
catchmentName = "CatchmentName"
catchmentRaster = catchmentName + ".tif"
catchmentRasterProj = catchmentName + "proj.tif"
catchment100 = catchmentName + "100.tif"
catchment100Fill = catchmentName + "100Fill.tif"
catchment100Flow = catchmentName + "100Flow.tif"
catchment100Acc = catchmentName + "100Acc.tif"
catchment100PP = catchmentName + "100PP.tif"
catchment100WS = catchmentName + "100WS.tif"
folder = "C:/temp/testIt"
folderAdd = "C:/temp/testIt/"
outlet = "outlet.shp"
projection = "C:/Program Files (x86)/ArcGIS/Desktop10.0/Coordinate Systems/Projected Coordinate Systems/National Grids/Europe/British National Grid.prj"
# identify guage
#create pourpoint shapefile
# identify catchment outline
# select appropriate testIt

# have a folder with all of the little ascii testIt for a catchment. This folder will be called *catchment name*
print "checking extensions"

try:
 arcpy.CheckExtension("Spatial") == "Available"
 arcpy.CheckOutExtension("Spatial")
 print "done"
except:
 print "Didn't work!"
 

print "converting the ascii files to raster files"

asciiList = os.listdir(folder)    

for x in range(len(asciiList)-1, -1, -1):
 if not asciiList.endswith(".asc"):
   asciiList.pop(x)


for file in asciiList:  
 print folderAdd + file
 print folderAdd + file[:-4]
 
 try:
  arcpy.ASCIIToRaster_conversion(folderAdd+ file, folderAdd+ file[:-4] + ".tif", "INTEGER")
 except:
  print "Processing " + file + "FAILED"
  
 time.sleep(3)
 
tifList = os.listdir(folder)

print tifList

for x in range(len(tifList) -1, -1, -1):
    if not tifList.endswith(".tif"):
        tifList.pop(x)

print tifList

tifString = ""
for tif in tifList:
 tifString += tif + ";"
 
tifString = tifString[:-1]

print tifString

print "mosaic to new raster"

try:
 arcpy.MosaicToNewRaster_management(tifString, folderAdd, catchmentRaster, "#", "16_BIT_SIGNED", "#", "1" , "MEAN")
except:
 print "Error with mosaic to new raster"
 
time.sleep(3)


print "Giving it the right coordinates"
try:
 arcpy.DefineProjection_management(catchmentRaster, projection)
 arcpy.ProjectRaster_management(catchmentRaster, catchmentRasterProj, projection, "BILINEAR", "#", "#","#",projection)
except Exception as e:
    print e.message


time.sleep(3)
 
print "Changing resolution"

try:
 outAggreg = Aggregate(folderAdd + catchmentRasterProj, 10, "MEDIAN", "TRUNCATE", "DATA")
 outAggreg.save(catchment100)
except Exception as e:
    print e.message
 
time.sleep(3)
 
print "Filling sinks"

try:
 outFill = Fill(catchment100)
 outFill.save(catchment100Fill)
except Exception as e:
    print e.message
 
print "Calculating fow direction"

time.sleep(3)

try:
 outFlowDirection = FlowDirection(catchment100Fill, "NORMAL")
 outFlowDirection.save(catchment100Flow)
except Exception as e:
    print e.message
 
 
print "Calculating flow accumulation"

time.sleep(3)

try:
 outFlowAccumulation = FlowAccumulation(catchment100Flow)
 outFlowAccumulation.save(catchment100Acc)
except Exception as e:
    print e.message

print "Creating pour point snap"

time.sleep(3)

try:
 outSnapPour = SnapPourPoint(outlet, catchment100Acc, 5, "FID")
 outSnapPour.save(catchment100PP)
except Exception as e:
    print e.message
 
time.sleep(3)

print "Creating watershed"

try:
 outWatershed = Watershed(catchment100Flow,catchment100PP)
 outWatershed.save(catchment100WS)
except Exception as e:
    print e.message
 
time.sleep(3)
 
print "Creating files for SHETRAN prepare"

try:
 arcpy.RasterToASCII_conversion(catchment100, catchmentName + "100DEM.txt")
 arcpy.RasterToASCII_conversion(catchment100WS, catchmentName + "100Mask.txt")
except Exception as e:
    print e.message
 
time.sleep(3)
 
print "Done!"

raw_input("close to exit")

Thanks for the help!

Re: Stand alone script fails using spatial analyst tools

MichaelGleason — Fri, 21 Jun 2013 15:41:36 GMT

I'm having a very similar problem over a year later -- arcpy.sa modules fail with 999999 errors on the S: drive, but run without issue on C:. I'm experiencing this at ArcGIS 10.0 SP 3 and ArcGIS 10.1 SP 1.

This seems like it should really be a bug report.