Latest Contributions by ChrisSnyder

Off topic, but... This is true at 9.3 but not at 10.0 ESRI never said all their SA tools support native read/wriite, and of course never indicated specifically which tools still rely on native grid. In my experience, many of the v10 hydrologic tools (fill, flowdir, etc.) are still much faster in GRID format. Unless someday ESRI says - "now in ArcGIS v23.1 all the SA tools are native FGDB read/write and are now all faster when used in FGDB format than GRID format" I will keep using good ole' GRID format as my default raster format of choice. ... View more

Here a sniipit of some Python code that builds "accumulated" elevation, slope, curvature, etc. grids - these grids are used as part of a DEM-based stream perdiction model. Note: flowDirLclGrd = a flow direction grid derived from a filled DEM. flowAcc1Grd = root + "\\" + areaFolderName + "\\" + origOidFolderName + "\\flowacc1" gp.FlowAccumulation_sa(flowDirLclGrd, flowAcc1Grd, "", "FLOAT"); showGpMessage() #Process: Builds a flow accumulation grid - accumulated acreage flowAcc2Grd = root + "\\" + areaFolderName + "\\" + origOidFolderName + "\\flowacc2" somaExp = "FlowAccumulation(" + flowDirLclGrd + ") * " + str(int(gp.cellsize) * int(gp.cellsize)) + " / 43560" gp.SingleOutputMapAlgebra_sa(somaExp, flowAcc2Grd); showGpMessage() #Process: Builds a flow accumulation grid - of estimated CFS using only watershed area and precip pcpGrd = root + "\\" + areaFolderName + "\\" + origOidFolderName + "\\pcp" gp.Resample_management(precipGrd, pcpGrd, gp.cellsize, "BILINEAR"); showGpMessage() flowAcc3Grd = root + "\\" + areaFolderName + "\\" + origOidFolderName + "\\flowacc3" somaExp = "flowaccumulation(" + flowDirLclGrd + ", " + pcpGrd + ", 'FLOAT') * " + str(int(gp.cellsize) * int(gp.cellsize)) + ") / 12 ) / 31556926" #31556926 = seconds/year gp.SingleOutputMapAlgebra_sa(somaExp, flowAcc3Grd); showGpMessage() #Process: Reclassifies any 0s flowAcc1Grd to 1 (which makes it so that we can use division (aka can't divide a grid by another grid if the denominator has a 0 in it) flowAcc5Grd = root + "\\" + areaFolderName + "\\" + origOidFolderName + "\\flowacc5" somaExp = "con(" + flowAcc1Grd + " == 0, 1, " + flowAcc1Grd + ")" gp.SingleOutputMapAlgebra_sa(somaExp, flowAcc5Grd); showGpMessage() #Process: Builds an accumulated elevation grid (uses the fil2LclGrd) accElvGrd = root + "\\" + areaFolderName + "\\" + origOidFolderName + "\\acc_elev" somaExp = "flowaccumulation(" + flowDirLclGrd + ", " + fil2LclGrd + ", 'FLOAT') / " + flowAcc5Grd gp.SingleOutputMapAlgebra_sa(somaExp, accElvGrd); showGpMessage()   #Process: Builds an accumulated slope grid accSlpGrd = root + "\\" + areaFolderName + "\\" + origOidFolderName + "\\acc_slp" somaExp = "flowaccumulation(" + flowDirLclGrd + ", " + slpGrd + ", 'FLOAT') / " + flowAcc5Grd gp.SingleOutputMapAlgebra_sa(somaExp, accSlpGrd); showGpMessage() #Process: Builds an accumulated curvature grid - fancy stuff is to deal with negative values accCurGrd = root + "\\" + areaFolderName + "\\" + origOidFolderName + "\\acc_cur" somaExp = "(flowaccumulation(" + flowDirLclGrd + ", max(" + curGrd + ", 0), 'FLOAT') - flowaccumulation(" + flowDirLclGrd + ", max(-1 * " + curGrd + ", 0), 'FLOAT')) / " + flowAcc5Grd gp.SingleOutputMapAlgebra_sa(somaExp, accCurGrd); showGpMessage() #Process: Builds an accumulated profile grid accProGrd = root + "\\" + areaFolderName + "\\" + origOidFolderName + "\\acc_pro" somaExp = "(flowaccumulation(" + flowDirLclGrd + ", max(" + proGrd + ", 0), 'FLOAT') - flowaccumulation(" + flowDirLclGrd + ", max(-1 * " + proGrd + ", 0), 'FLOAT')) / " + flowAcc5Grd gp.SingleOutputMapAlgebra_sa(somaExp, accProGrd); showGpMessage() #Process: Builds an accumulated plan grid accPlnGrd = root + "\\" + areaFolderName + "\\" + origOidFolderName + "\\acc_pln" somaExp = "(flowaccumulation(" + flowDirLclGrd + ", max(" + plnGrd + ", 0), 'FLOAT') - flowaccumulation(" + flowDirLclGrd + ", max(-1 * " + plnGrd + ", 0), 'FLOAT')) / " + flowAcc5Grd  gp.SingleOutputMapAlgebra_sa(somaExp, accPlnGrd); showGpMessage() ... View more

A feature layer (what you are calling a virtual layer I think) is nothing more than a pointer to a feature class on disk. As a best practice, always add the field(s) to the feature class 1st, then make a feature layer from the feature class. The eventual users of the script/model will not have write access to the original features, so the tool fails since it can not write to the original table. What you probably want to do in this case is copy the original feature classes to a new workspace where the users will have write access. A very good option (if the datasets aren't too big) is the in_memory workspace: http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//002w0000005s000000.htm ... View more

The MakeQueryTable tool can do this: http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00170000006r000000.htm The syntax is a little squirly. Be sure to read the documentation. Also, be sure to include the Shape field if you want a spatial layer returned, otherwise it'll just return a table view. This might help too: http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Examples_of_queries_with_the_Make_Query_Table_tool/0017000000v9000000/ ... View more

Wow - Very nice! I'd have to take a while to digest all that... I agree, there is a significant amount of overhead splitting things up, managing the separate processes (multiprocess seems to make it easier), and then squishing it all back together. Not for the faint hearted... If you are interested, here's my stone age attempt at parallel processing using os.spawnv (this code is almost 5 years old now!). I have a large process that I run every 3 months or so for my bosses, just a big overlay of our land base (riparian areas, forest inventory, harvest areas, habitat concerns, etc.) that we use for reporting, forest modeling, etc. Originally, due to the size of all the layers involved, the union wouldn't run at all, so I create a "tile" feature class (composed of ~ 40 little rectangular polygons, and then would run the overlay for each tile - one at a time. Then I figured, hey, I could run many tiles at once (but just had to write some code to regulate the timing of it all). Instatnt paralell process! Anyway, it's pretty stone age (uses .txt files to comunicate!), but it works. The "run_union_slave" script does the actual geoprocessing, this "master" script below just manages the overall parallel process. # Description # ----------- # This master script regulates the execution of the script run_union_slave.py (for each tile) # Author: Chris Snyder, WA Department of Natural Resources, chris.snyder(at)wadnr.gov import sys, string, os, shutil, time, traceback, glob try:     #Defines some functions     def showPyMessage():         try:             print >> open(logFile, 'a'), str(time.ctime()) + " - " + str(message)             print str(time.ctime()) + " - " + str(message)         except:             pass     def launchProcess(tileNumber):         global message         global numberOfProcessors         message = "Processing tile #" + str(tileNumber); showPyMessage()         #Added this to address bug with indexTiles.shp used in the overlay tools in v9.2+ (can't run more than 1 overlay process at once)         newTempDir = r"C:\temp\gptmpenvr_" + time.strftime('%Y%m%d%H%M%S')         os.mkdir(newTempDir)         os.environ["TEMP"] = newTempDir         os.environ["TMP"] = newTempDir         message = "Set TEMP and TMP variables to " + newTempDir; showPyMessage()         #End bug fix...         parameterList = [pythonExePath, slaveScript, root, str(tileNumber), yyyymmdd]         if tileNumber <= indxPolyFeatureCount:             if tileNumber == indxPolyFeatureCount or numberOfProcessors == 1:                 message = "Waiting for tile #" + str(indxPolyFeatureCount) + " to finish..."; showPyMessage()                 os.spawnv(os.P_WAIT, pythonExePath, parameterList)             else:                 os.spawnv(os.P_NOWAIT, pythonExePath, parameterList)         time.sleep(1); message = "Waiting a few seconds"; showPyMessage()         #Specifies the root variable, makes the logFile variable, and does some error checking...     dateTimeStamp = time.strftime('%Y%m%d%H%M%S')     root = sys.argv[1] #r"E:\1junk\overlay_20060531"     yyyymmdd = sys.argv[2] #this is passed from the master script since the tiles may be processed over a 2+ day period     if os.path.exists(root)== False:         print "Specified root directory: " + root + " does not exist... Bailing out!"         sys.exit()     scriptName = sys.argv[0].split("\\")[len(sys.argv[0].split("\\")) - 1][0:-3] #Gets the name of the script without the .py extension      logFile = root + "\\log_files\\" + scriptName + "_" + dateTimeStamp[:8] + ".log" #Creates the logFile variable     if os.path.exists(root + "\\log_files") == False: #Makes sure log_files exists         os.mkdir(root + "\\log_files")     if os.path.exists(logFile)== True:         os.remove(logFile)         message = "Deleting log file with the same name and datestamp... Recreating " + logFile; showPyMessage()     workspaceDir = root + "\\index_tiles"     if os.path.exists(workspaceDir)== True:         message = "Overlay directory: " + workspaceDir + " already exist... Deleting and recreating " + workspaceDir; showPyMessage()         shutil.rmtree(workspaceDir)     else:         message = "Creating index tiles directory: " + workspaceDir; showPyMessage()     os.mkdir(workspaceDir)             #Process: Finds python.exe     pythonExePath = ""     for path in sys.path:         if os.path.exists(os.path.join(path, "python.exe")) == True:             pythonExePath = os.path.join(path, "python.exe")     if pythonExePath != "":         message = "Python.exe file located at " + pythonExePath; showPyMessage()     else:         message = "ERROR: Python executable not found! Exiting script...."; showPyError(); sys.exit()     #Determines the number of processors on the machine...     numberOfProcessors = int(os.environ.get("NUMBER_OF_PROCESSORS"))     maxNumberOfProcessors = 3 #The maximum number of processors you want to use     if numberOfProcessors > maxNumberOfProcessors:         numberOfProcessors = maxNumberOfProcessors     #Determines the number of tiles in the index layer (tile numbers are assumed to be sequential)     #Note: I didn't use gp.getcount so that this script wouldn't use the gp at all (and eat up more memory)     indxPolyCountFile = glob.glob(root + "\\log_files\\therearethismanytiles_*.txt")     if len(indxPolyCountFile) == 0:         message = "Can't find index polygon count .txt file in " + root + "\\log_files" + "! Exiting script..."; showPyMessage()         sys.exit()     indxPolyFeatureCount = int(indxPolyCountFile[0].split("_")[-1][0:-4])     if indxPolyFeatureCount == 0:         message = "Index polygon count is 0! Exiting script..."; showPyMessage()         sys.exit()         tileNumber = 1     numberOfProcesses = 0     slaveScript = r"\\Snarf\am\div_lm\ds\gis\ldo\current_scripts\run_union_slave_v93.py"     tilesThatAreProcessingList = []     tilesThatFinishedList = []     tilesThatBombedList = []         #Process: This while loop initialy launches <numberOfProcessors> instances of the slave script      while numberOfProcesses < numberOfProcessors:         numberOfProcesses = numberOfProcesses + 1         tilesThatAreProcessingList.append(tileNumber)         launchProcess(tileNumber)         tileNumber = tileNumber + 1     #Process: This while loop checks for tiles that finished or bombed, if there are, new slave scripts are launched     while len(tilesThatBombedList) + len(tilesThatFinishedList) < indxPolyFeatureCount:         isDoneList = glob.glob(root + "\\log_files\\" + slaveScript.split("\\")[len(slaveScript.split("\\")) - 1][0:-3] + "_isalldone_*") #makes a list of the .txt file created by the slave script         doneBombedList = glob.glob(root + "\\log_files\\" + slaveScript.split("\\")[len(slaveScript.split("\\")) - 1][0:-3] + "_bombed_*")         if len(isDoneList) == 0 and len(doneBombedList) == 0: #if there are no .txt files, wait for x number of seconds             time.sleep(5)         else: #else if there are .txt files...             if len(isDoneList) > 0: #if there are tiles that are "_isalldone_"                 for isDoneItem in isDoneList: #for each .txt file indicating completion...                     tileThatIsDone = int(isDoneItem[isDoneItem.rfind("_") + 1:isDoneItem.rfind(".")])                     os.remove(root + "\\log_files\\" + slaveScript.split("\\")[len(slaveScript.split("\\")) - 1][0:-3] + "_isalldone_" + str(tileThatIsDone) + ".txt")                     message = "Tile #" + str(tileThatIsDone) + " is done!!!"; showPyMessage()                     tilesThatFinishedList.append(tileThatIsDone)                     tilesThatAreProcessingList.remove(tileThatIsDone)                     tilesThatAreProcessingList.append(tileNumber)                     launchProcess(tileNumber)                     tileNumber = tileNumber + 1                     time.sleep(5)             if len(doneBombedList) > 0: #if there are tiles that "_bombed_"                   for doneBombedItem in doneBombedList: #for each .txt file indicating failure...                     tileThatBombed = int(doneBombedItem[doneBombedItem.rfind("_") + 1:doneBombedItem.rfind(".")])                     os.remove(root + "\\log_files\\" + slaveScript.split("\\")[len(slaveScript.split("\\")) - 1][0:-3] + "_bombed_" + str(tileThatBombed) + ".txt")                     message = "Tile #" + str(tileThatBombed) + " bombed!!!"; showPyMessage()                     tilesThatBombedList.append(tileThatBombed)                     tilesThatAreProcessingList.remove(tileThatBombed)                     tilesThatAreProcessingList.append(tileNumber)                     launchProcess(tileNumber)                     tileNumber = tileNumber + 1                     time.sleep(5)             message = "Tiles that are currently processing: " + str(tilesThatAreProcessingList); showPyMessage()             message = "Tiles that are done: " + str(tilesThatFinishedList); showPyMessage()             message = "Tiles that bombed: " + str(tilesThatBombedList); showPyMessage()             time.sleep(5)     if len(tilesThatBombedList) > 0:         message = "ERROR - these tiles failed to process: " + str(tilesThatBombedList); showPyMessage()         sys.exit()             message = "ALL DONE!"; showPyMessage()     print >> open(root + "\\log_files\\" + scriptName + "_isalldone.txt", 'a'), scriptName + "_isalldone.txt" except:     message = "\n*** LAST GEOPROCESSOR MESSAGE (may not be source of the error)***"; showPyMessage()     message = "\n*** PYTHON ERRORS *** "; showPyMessage()     message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()     message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage() ... View more

Read this: http://forums.arcgis.com/threads/17956-Using-centroid-geometry-to-get-values-of-intersecting-feature?p=58625&viewfull=1#post58625 and the posts below it. Basically, we came to the same conclusion: You can't directly store/use the geometry object in a Python dictionary. ... View more

If you create a simple GRID you will not have any kind of mistake. I think team04 meant that by using good ole' grid format as an input/output format (not a FGDB raster), the problem was avoided. Just an FYI too, but generally, reading writting grid format in SA tools is way faster than FGDB raster format. If you use FGDB format stuff, behind the scenes it is formats to grid, runs the SA tool algorithm, and then, writes the grid output back to FGDB raster fromat. There is a lot of time that is potentially lost in all that reformating... ... View more

Just select your code text and click the little "#" button (3rd to the left of the YouTube button). ... View more

In this case I am able to split the data beforehand, then define a python module that does the processing on the section of data that is passed to it; it's these things that run in parallel. Let me know if you want more info - I can prepare a sample code for you - there's not much around on the internet for this kind of thing! Stacy, I'd love to see some samples. Ocassionally, I have the need to write large number-crunching stuff in Python (traversing large arrays or dictionary-type structures), and I'd love to explore ways to speed it up the processing. Youre' right, it's pretty hard to get practical examples of using some of this stuff! ... View more

2) You can pass a list objects (maybe tuples?) directly into the v93 gp tools Looks like this statement (at least in v9.3.1) only applies to geometry objects - not on disk FCs, so I wasn't completely wrong... For example: inputList = [r"D:\csny490\test.gdb\townships", r"D:\csny490\test.gdb\regions"] outputFC =  r"D:\csny490\test.gdb\output" gp.intersect_analysis(inputList, outputFC) errors out, but inputList = [r"D:\csny490\test.gdb\townships", r"D:\csny490\test.gdb\regions"] outputFC =  r"D:\csny490\test.gdb\output" gp.intersect_analysis(";".join(inputList), outputFC) works fine. ... View more

Jörg, per your original post, how about: gp.Intersect_analysis([o1,o2],output) ... View more

Well, live and learn. I guess that's why I look at these forum posts... Contrary to my beliefs 5 minutes ago: 1) You can use geomety objects directly in v9.3 (thought that was new in v10) 2) You can pass a list objects (maybe tuples?) directly into the v93 gp tools, and apparently no longer have to use string parameters. #snipit of ESRI example code inGeom1 = gp.createobject("geometry", "polygon", ary) inGeom2 = gp.createobject("geometry", "polygon", ary2) inGeoms = [inGeom1, inGeom2] outGeom = gp.createobject("geometry") newOutGeom = gp.union_analysis(inGeoms, outGeom) Good to know... ... View more

Seeing that Arcpy isn't a native Python object, would there be any benifit (or would it even work) to send an arcpy.Whatever commands to Parallel Python, Pickle, etc.? The only way I've ever gotten a "parallel" process to actually work with gp/arcpy objects is using os.spawnv and/or subprocess. ... View more

You are attempting to pass an actual Python list into the Append tool. While this seems like it should work, it doesn't. You need to format the list as a big long string. For example: fcString = "" for fc in fcList:    fcString = fcString + fc + ";" arcpy.Append_management(fcString[:-1], outLocation + os.sep + emptyFC, "NO_TEST") Another way: arcpy.Append_management(str(fcList).replace(",",";")[1:-1], outLocation + os.sep + emptyFC, "NO_TEST") ... View more

The issue you are having is that you need to list EVERY field in the table and indicate "ORIGINAL_NAME NEW_NAME HIDDEN/VISIBLE" Maybe something like this: keepFieldList = ("FIELD1","FIELD2") fieldInfo = "" fieldList = gp.listfields(layerPath) for field in fieldList:     if field.name in keepFieldList:         fieldInfo = fieldInfo + field.Name + " " + field.name + " VISIBLE;"     else:         fieldInfo = fieldInfo + field.Name + " " + field.Name + " HIDDEN;" gp.MakeFeatureLayer_management(layerPath, "feature_layer", "", "", fieldInfo[:-1]); showGpMessage() ... View more