Is there a site with templates for Arcpy where I can just fill in the blanks?

Question

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.

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.

Answer 1

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.

Answer 2

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.

Answer 3

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.

Answer 4

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)