1

I have a Shapefile containing near 7000 pour points across the globe (GCS, WGS84). I'm using the following model to select the first feature point of the Shapefile, copy the selected feature in a new Shapefile, and then make its watershed. I need to do this because using the Watershed tool directly doesn't take account of the overlapping watersheds. Every thing is explained in my previous question.

enter image description here

The model works well, but it's painfully slow. The first two tools (Iterate and Copy) take less then a second to process, but the last tool (Watershed) takes up to 30 minutes for each watershed. In this question, I was advised to use Arcpy instead of Model Builder, but I've never used Arcpy or Python before.

So here's my question : Is there a site with templates for Arcpy where I can just fill in the blanks?

I'm looking at online tutorials, but since I'm in a hurry, having a pre-coded script would speed thinks up. All I need is a watershed script with a looping of the pour point used.

  • Please open a new question rather than editing the question with a follow-up question. – Aaron Jun 16 '15 at 1:47
2

From ModelBuilder:

Model > Export > To Python Script...

This will create a python script that does exactly what you specified in the model. You can then look at it and modify it as desired. I'm not sure you will get any speedup from this, since using python calls all the same tools as ModelBuilder, so the watershed tool will still take the same amount of time to execute.

  • 6
    This will fall down on models like this which use iterators. Export to Python Script doesn't do a 1:1 functional match of the model to a Python script. It's broken and as far as I know Esri doesn't plan on fixing it. You're safer copying individual tools as Python snippets from the results window and piecing together a script that way. – Jason Scheirer Jun 15 '15 at 17:12
  • You're right, I hadn't thought about that. The iterator will have to be replaced with a simple python loop. This is what is causing the Model Functions error you are receiving. You'll have to remove the line that say: arcpy.ImportToolbox("Model Functions"), as well as the line with IterateFeatureSelection. The remaining code (other than imports/license checkouts) should be enclosed in a loop which iterates over the features. – BruceDoh Jun 15 '15 at 20:22
1

If your model does what you need, Modelbuilder already has a built-in function that will sort of do what you're asking for. In the Modelbuilder window, go to Model > Export > To Python Script. That will convert your model steps to the equivalent Python commands in a script file. This process is automated by ArcGIS, so it doesn't necessarily export the most pretty code possible, and sometimes the code doesn't work without a little tweaking, but in your situation I'd try that first and run a few tests of the exported Python script to make sure it produces the same output you expect from your model and go from there.

1

These options may help or may not depending on your hardware configuration and license/tool availability.

(1) You could look into enabling 64-bit background processing to help things along. If these tools are able to utilize it, it could potentially shave some time off for you.

(2) You might be able to (depending on the availability of the tools and your organizations licensing) look into ArcGIS Pro since it is a native 64-bit software. I haven't played with this much since my organization doesn't have the appropriate license to let me, but if the tools are available for you, it may be a better option than standard ArcMap.

(3) If you export your model into a python script, you can try to run it outside of ESRI software in another python IDE such as PyScripter. I have found that a lot of times, scripts run faster outside of ESRI software. They go through and check your licensing for you so you need not worry about that. Just import arcpy and you should be good.

  • I'm trying option 3 (PyScripter) and I got the following error : "IOError: The toolbox file Model Functions was not found.". How can I fix this? Thank you! – Nahas Jun 15 '15 at 18:48
  • This might be a problem with the export to python script part. You should be able to call your toolbox functions via arcpy although the Model Functions one might be referenced differently. Have you scanned through your code to find anything that looks like it is referencing that? – Branco Jun 15 '15 at 19:20
1

The trick here to work with smaller extracts from massive flow direction raster. The script below does exactly this, growing processing extent from initial buffer (radius) around point equal to (1000*cell size). When no growth in catchment size found it proceeds to next pour point, saving output in wshed_0000.shp, wshed_0001.shp etc in outFolder folder :)

Script looks stupid because I've tried to use multiprocessor and failed. Couldn't force arcpy to recognise availability of spatial analyst extension in multiprocessor pool. This is why add this script to your toolbox with no parameters.

On my rather good machine it took 9 minutes to process 3 points (each took 5 iterations). I suggest try first with 3 points only and change nPoints to 7000, if you are happy with the performance. Also it is worth including extreme size catchments in this original subset to see if 1000 cells is a good 1st approximation, change it if it result in to many iterations for supposedly smallest catchment.

Good question! Good luck!

# Description:
# muliple watersheds
# ---------------------------------------------------------------------------
# Import arcpy module
import arcpy
from arcpy import env
from arcpy.sa import Watershed

# things to change
pourPoints=r'd:\scratch\pour_points.shp'
fDir=r'd:\scratch\fdir'
outFolder=r'd:\aerials\images'
radius=1000
nPoints=3
nProcessors=4

env.workspace = outFolder
env.pyramid='NONE'
env.rasterStatistics = 'NONE'
fDir=arcpy.Raster(fDir)
env.snapRaster = fDir
cSize=fDir.meanCellHeight
env.cellSize = cSize
ext=arcpy.Describe(fDir).extent
extract,wshed,rdlayer='xtrct','wshed','RL'

def function(inputs):
    fid,pnt=inputs
    suffix=str(fid).zfill(4)
    quer='"FID" = '+str(fid);
    thePoint='layer'+suffix
    arcpy.MakeFeatureLayer_management(pourPoints,thePoint,quer)
    outFile='wshed_'+suffix+'.shp'
    i,nTotal=radius,1
    while True:
        bFirst=cSize*i
        arcpy.AddMessage('Trying buffer %i' %bFirst)
        buf=pnt.buffer(bFirst)
        buf=buf.clip(ext)
        anExtent=buf.extent
        env.extent=anExtent
        envelope='%f %f %f %f' %(anExtent.XMin, anExtent.YMin, anExtent.XMax, anExtent.YMax,)
        arcpy.Clip_management(fDir,envelope,extract,buf)
        oneGrid=Watershed (extract, thePoint, "FID")
        oneGrid.save(wshed)
        lmax=arcpy.Raster(wshed)
        arcpy.MakeRasterLayer_management(lmax, rdlayer)
        tbl=arcpy.da.TableToNumPyArray(rdlayer,("value","count"))
        fLine=tbl[0]
        newTotal=fLine[1]
        if newTotal==nTotal: break
        nTotal=newTotal
        i=2*i
    arcpy.RasterToPolygon_conversion(lmax, outFile, "NO_SIMPLIFY", "Value")
    for elm in [extract,wshed]:
        arcpy.Delete_management(elm)
    arcpy.AddMessage ('Processed %i' %fid)
    return outFile
if __name__ == "__main__":
    g = arcpy.Geometry()
    geometryList=arcpy.CopyFeatures_management(pourPoints,g)
    for i in range(0,nPoints,nProcessors):
        bList,m = [],0
        for j in range(i,min(i+4,nPoints)):
            m+=1; bList.append([j,geometryList[j]])
        for elm in bList:
            out=function(elm)
  • Thanks! But, I got the following errors in the Output of PyScript (same in ArcGIS) : Traceback (most recent call last): line 21, in <module> env.snapRaster = fDir line 529, in set_ self[env] = val line 581, in __ setitem __ ret_ = setattr(self._gp, item, value) RuntimeError: Object : Error in accessing environment <snapRaster> – Nahas Jun 16 '15 at 15:42
  • Something wrong with your raster. Find it in catalog, recalculate extent and try again. From arcgis – FelixIP Jun 16 '15 at 19:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.