nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

loop down through a directory tree.

PaulHuffman

I want to loop through a folder and clip all feature classes by a AOI poly and write the results to an output folder. How to I also have the loop search all subfolders for vectors and raster feature classes as well?

Trouver plus de publications taguées avec

Accepted answers

All comments

PaulHuffman

Fixed it. This is what satisfied the colleague, and was put into production. Lines 30-34 is how I check to see if the current worksheet is a Grid or a workspace containing grids. Had to go through this extra BS because I couldn't use da.walk.

#testclipgridstifsshapes.py
#test trying to walk a directory tree, find grids, and clip to a area of Interest (AOI)
# and make it run in 10.0 so without da.walk.
# Paul Huffman, Yakama Fisheries, 7/21/2014

# Import arcpy module
import arcpy
import os

# Set Geoprocessing environments
#arcpy.env.scratchWorkspace = "c:\\avdata\\PythonTest\\ClipLidar\\Scratch"
arcpy.env.overwriteOutput = True
clipshape = "c:\\avdata\\PythonTest\\ClipLidar\\AOIPoly_SP.shp"
InputFolder = "c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Input"

#find input folder bottom folder
InputString = InputFolder.split("\\")[-1]

#Find AOI extent
AOI = arcpy.mapping.Layer(clipshape)
AOIextent=AOI.getExtent()

#Intialize empty lists
rasterList=[]

#Find the input  Grid feature classes
for root, dirs, files in os.walk(InputFolder):
    for targetdir in dirs:
        targetworkspace = root + "\\" + targetdir
        #if the targetworkspace can be initialized as a Raster object, do nothing else
        try: 
            arcpy.Raster(targetworkspace)
         #if the targetworkspace fails when an attempt is made to initialize it as a raster object, search for rasters within it
        except:
            arcpy.env.workspace = targetworkspace
            rasterInDir = arcpy.ListRasters("*", "GRID")
            if len(rasterInDir) > 0:
                for raster in rasterInDir:
                    rasterList.append(targetworkspace + "\\" + raster)
            
##for raster in rasterList:
##        print raster

#start comparing extents and clipping
for name in rasterList:
    FileObj = arcpy.mapping.Layer(name)
    FileExtent = FileObj.getExtent()
    #print "Processing: " + name
   #if image overlap or touches AOI polygon,  then clip image, write to output
    Extent_Disjoint = str(FileExtent.disjoint(AOIextent))
    if Extent_Disjoint == 'False':
        print "Clipping: " + name
        foldername = os.path.dirname(name)
        OutputPath = foldername.replace(InputString,"Output")
        if not os.path.exists(OutputPath):
            os.makedirs(OutputPath)
        if len(os.path.basename(name)) > 10:    
            Shortname = (os.path.basename(name))[0:11]
            print "Grid name shortened to " + Shortname
            Clipname = os.path.splitext(Shortname)[0] + "cp"
        else:    
            Clipname = os.path.splitext(name)[0] + "cp"
        Clipbasename = os.path.basename(Clipname)
        print "Writing " + (os.path.join(OutputPath,Clipbasename))
        arcpy.Clip_management(name,"#",(os.path.join(OutputPath,Clipbasename)),clipshape,"#","ClippingGeometry")

#Find the input tifs and shapefiles
for root, dirs, files in os.walk(InputFolder):
    for name in files:
        if os.path.splitext(name)[1] == ".tif":
            FileObj = arcpy.mapping.Layer(os.path.join(root, name))
            FileExtent = FileObj.getExtent()
            #print "Processing: " + name
           #if image overlap or touches AOI polygon,  then clip image, write to output
            Extent_Disjoint = str(FileExtent.disjoint(AOIextent))
            if Extent_Disjoint == 'False':
                print "Clipping: " + name
                OutputPath = root.replace(InputString,"Output")
                if not os.path.exists(OutputPath):
                    os.makedirs(OutputPath)
                Clipname = os.path.splitext(name)[0] + "cp.tif"
                print "Writing " + (os.path.join(OutputPath,Clipname))
                arcpy.Clip_management((os.path.join(root, name)),"#",(os.path.join(OutputPath,Clipname)),clipshape,"#","ClippingGeometry")
        #for shapefiles, use a different clip, write to output. 
        if os.path.splitext(name)[1] == ".shp":
            print(os.path.join(root, name))
            FileObj = arcpy.mapping.Layer(os.path.join(root, name))
            FileExtent = FileObj.getExtent()
            Extent_Disjoint = str(FileExtent.disjoint(AOIextent))
            if Extent_Disjoint == 'False':
                print "Clipping: " + name
                OutputPath = root.replace(InputString,"Output")
                if not os.path.exists(OutputPath):
                    os.makedirs(OutputPath)
                Clipname = os.path.splitext(name)[0] + "cp"
                print "Writing " + (os.path.join(OutputPath,Clipname))
                arcpy.Clip_analysis((os.path.join(root, name)),clipshape,(os.path.join(OutputPath,Clipname)))

PaulHuffman

There is still something wrong with the way arcpy.ListRasters works to detect grids. Maybe not wrong, but unexpected. Line 12's output shows that ListRasters still finds my grids in their parent workspace and also the subdirectory for the grid that contains the info external files.

The code segment:

#Find the input  Grid feature classes
for root, dirs, files in os.walk(InputFolder):
    for targetdir in dirs:
        targetworkspace = root + "\\" + targetdir
        arcpy.env.workspace = targetworkspace
        rasterInDir = arcpy.ListRasters("*", "GRID")
        if len(rasterInDir) > 0:
            for raster in rasterInDir:
                rasterList.append(targetworkspace + "\\" + raster)
            
for raster in rasterList:
        print raster

The output of the print statement, in part:

c:\avdata\PythonTest\ClipLidar\TestData\Input\BareEarth\hillsha_tane1
..
..
c:\avdata\PythonTest\ClipLidar\TestData\Input\BareEarth\hillsha_tane1\hillsha_tane1

PaulHuffman

This is what I have now. Works with tifs and shapefiles, but I need to figure out how to detect and process grids as well. The script looks for files down through the tree below "c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Input" then writes the data into a copy of the input tree structure below "c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Output" To migrate this script to another datas et, the end user has to edit in the new clipshape AOI, the new input path, then also change the name of the top of the input tree "Input" in line 48. This is too tricky. I'd rather have one section at the top where all the editable parameters are found. Maybe I can generalize this by somehow getting the name of that bottom rung folder out of the InputFolder string and stick it in a variable.

Lines 27 - 39 are redundantly repeated again redundantly at line 43. Bad style. These lines could define a function, except for the differences in lines 37 and 39.

#testclip.py
#test trying to walk a directory tree, find rasters and shapefiles, and clip to a area of Interest (AOI)
# and make it run in 10.0 so without da.walk.
# Paul Huffman, Yakama Fisheries, 7/7/2014

# Import arcpy module
import arcpy
import os
#arcpy.CheckOutExtension("spatial")

# Set Geoprocessing environments
arcpy.env.scratchWorkspace = "c:\\avdata\\PythonTest\\ClipLidar\\Scratch"
arcpy.env.workspace = "c:\\avdata\\PythonTest\\ClipLidar\TestData"
arcpy.env.overwriteOutput = True
clipshape = "c:\\avdata\\PythonTest\\ClipLidar\\AOIPoly_SP.shp"
InputFolder = "c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Input"
#OutputFolder = "c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Output"

#Find AOI extent
AOI = arcpy.mapping.Layer(clipshape)
AOIextent=AOI.getExtent()

#Find the input feature classes
for root, dirs, files in os.walk(InputFolder):
    for name in files:
        if os.path.splitext(name)[1] == ".tif":
            FileObj = arcpy.mapping.Layer(os.path.join(root, name))
            FileExtent = FileObj.getExtent()
            print "Processing: " + name
           #if image overlap or touches AOI polygon,  then clip image, write to output
            Extent_Disjoint = str(FileExtent.disjoint(AOIextent))
            if Extent_Disjoint == 'False':
                print "Clipping: " + name
                OutputPath = root.replace("Input","Output")
                if not os.path.exists(OutputPath):
                    os.makedirs(OutputPath)
                Clipname = os.path.splitext(name)[0] + "cp.tif"
                print "Writing " + (os.path.join(OutputPath,Clipname))
                arcpy.Clip_management((os.path.join(root, name)),"#",(os.path.join(OutputPath,Clipname)),clipshape,"#","ClippingGeometry")
        #for shapefiles, use a different clip, write to output. 
        if os.path.splitext(name)[1] == ".shp":
            print(os.path.join(root, name))
            FileObj = arcpy.mapping.Layer(os.path.join(root, name))
            FileExtent = FileObj.getExtent()
            Extent_Disjoint = str(FileExtent.disjoint(AOIextent))
            if Extent_Disjoint == 'False':
                print "Clipping: " + name
                OutputPath = root.replace("Input","Output")
                if not os.path.exists(OutputPath):
                    os.makedirs(OutputPath)
                Clipname = os.path.splitext(name)[0] + "cp"
                print "Writing " + (os.path.join(OutputPath,Clipname))
                arcpy.Clip_analysis((os.path.join(root, name)),clipshape,(os.path.join(OutputPath,Clipname)))

PaulHuffman

I got this to work for at least tifs and shapefiles. I'm not sure why I had to use mapping.layer, line 20 and 29, before I could get an Extent of the data. Seems like the overlap tests, lines 32 - 47, could be down more efficiently. Is there just one test that would find if there is no touching or overlap, then skip that file. Where's the docs on that .overlaps stuff? Can't find it today.

Currently all my output lands in one folder. My next step needs to be writing the output to separate folders in the output directory, making the output tree structure look like the input tree. I'm guessing that will involve manipulating the strings of the full paths of the inputs.

#testclip.py
#test trying to walk a directory tree, find rasters and shapefiles, and clip to a area of Interest (AOI)
# and make it run in 10.0 so without da.walk.
# Paul Huffman, Yakama Fisheries, 7/7/2014

# Import arcpy module
import arcpy
import os
#arcpy.CheckOutExtension("spatial")

# Set Geoprocessing environments
arcpy.env.scratchWorkspace = "c:\\avdata\\PythonTest\\ClipLidar\\Scratch"
arcpy.env.workspace = "c:\\avdata\\PythonTest\\ClipLidar\TestData"
arcpy.env.overwriteOutput = True
clipshape = "c:\\avdata\\PythonTest\\ClipLidar\\AOIPoly_SP.shp"
InputFolder = "c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Input"
OutputFolder = "c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Output"

#Find AOI extent
AOI = arcpy.mapping.Layer(clipshape)
AOIextent=AOI.getExtent()

#Find the input feature classes
for root, dirs, files in os.walk(InputFolder):
    for name in files:
        Subset = 0
        if os.path.splitext(name)[1] == ".tif":
            #print os.path.join(root, name), '\n'
            FileObj = arcpy.mapping.Layer(os.path.join(root, name))
            FileExtent = FileObj.getExtent()
            print "Processing: " + name
            Extent_Overlap = str(FileExtent.overlaps(AOIextent))
            Extent_Within = str(FileExtent.within(AOIextent))
            Extent_Touches = str(FileExtent.touches(AOIextent))
            Extent_Crosses = str(FileExtent.crosses(AOIextent))
            if Extent_Overlap == 'True':
                Subset = 1
                print "Extent_Overlaps"
            elif Extent_Within == 'True':
                Subset = 1
                print "Extent_Within"
            elif Extent_Touches == 'True':
                Subset = 1
                print "Extent_Touches"
            elif Extent_Crosses == 'True':
                Subset = 1
                print "Extent_Crosses"
            if Subset == 1:
                arcpy.Clip_management((os.path.join(root, name)),"#",(os.path.join(OutputFolder,name)),clipshape,"#","ClippingGeometry")
        if os.path.splitext(name)[1] == ".shp":
            print(os.path.join(root, name))
            arcpy.Clip_analysis((os.path.join(root, name)),clipshape,(os.path.join(OutputFolder,name)))

PaulHuffman

I was trying to get this script done yesterday, but puzzled why half of it doesn't work. Because ESRI is doing site maintenance, I couldn'tget back in to ask this question. Client still waiting.

I was trying to get a script to walk down a tree of LIDAR data, find any rasters and feature classes, and clip them to an area of interest. My test data tree just has some tifs in one folder, some shapefiles in another to make things simpler. And I need to make this run in ArcGIS 10.0, so da.walk and da.Cursor are not available.

When my script hit the first tif file, the AOI polygon didn't overlap it at all and I got an error message that it didn't overlap. (Failed to create raster dataset. The reason could be: The clip feature is outside the raster extent.

Failed to execute (Clip).) I was going to put off solving the problem of testing for overlap until I could get the file finding and clipping figured out. So, following an example at http://gis.stackexchange.com/questions/13571/arcpy-compare-extents-of-rasters, I added the commented out section, lines 26 to 49, to check for overlap of extent for the tifs. However, the script hits line 48, and says name 'Subset' is not defined. Also none of the print statements in the section are executed, so it looks like none of the extent tests worked. What went wrong? If I just run the script down to the first print statement, all the tifs are found and printed.

In the next section, an if statement looks for just shapefiles, and successfully clips them and writes the output to a new output tree. No check for overlap is done for the shapefiles but I probably should to eliminate empty output shapefile and to save execution time.

Hey, how do you insert code in this new forum?

#testclip.py

#test trying to walk a directory tree, find rasters and shapefiles, and clip to a area of Interest (AOI)

# and make it run in 10.0 so without da.walk.

# Paul Huffman, Yakama Fisheries, 7/7/2014

# Import arcpy module

import arcpy

import os

arcpy.CheckOutExtension("spatial")

# Set Geoprocessing environments

arcpy.env.scratchWorkspace = "c:\\avdata\\PythonTest\\ClipLidar\\Scratch"

arcpy.env.workspace = "c:\\avdata\\PythonTest\\ClipLidar\TestData"

arcpy.env.overwiteOutput = True

clipshape = "c:\\avdata\\PythonTest\\ClipLidar\\AOIPoly.shp"

InputFolder = "c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Input"

OutputFolder = "c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Output"

#Find AOI extent

AOI = arcpy.mapping.Layer(clipshape)

AOIextent=AOI.getExtent()

#Find the input feature classes

for root, dirs, files in os.walk(InputFolder):

for name in files:

if os.path.splitext(name)[1] == ".tif":

print os.path.join(root, name), '\n'

FileObj = arcpy.mapping.Layer(os.path.join(root, name))

FileExtent = FileObj.getExtent()

print "Processing: " + name

print FileExtent

## Extent_Overlap = str(FileExtent.overlaps(AOIextent))

## Extent_Within = str(FileExtent.within(AOIextent))

## Extent_Touches = str(FileExtent.touches(AOIextent))

## Extent_Crosses = str(FileExtent.crosses(AOIextent))

## if Extent_Overlap == 'True':

## Subset = 1

## print "Extent_Overlaps"

## elif Extent_Within == 'True':

## Subset = 1

## print "Extent_Within"

## elif Extent_Touches == 'True':

## Subset = 1

## print "Extent_Touches"

## elif Extent_Crosses == 'True':

## Subset = 1

## print "Extent_Crosses"

## #print(os.path.join(root, name))

## #arcpy.Clip_management(name,clipshape,(os.path.join(OutputFolder,name)))

## if Subset == 1:

## arcpy.Clip_management((os.path.join(root, name)),"#",(os.path.join(OutputFolder,name)),clipshape,"#","ClippingGeometry")

## if os.path.splitext(name)[1] == ".shp":

## print(os.path.join(root, name))

## #arcpy.Clip_management((os.path.join(root, name)),"#",(os.path.join(OutputFolder,name)),clipshape,"#","ClippingGeometry")

## arcpy.Clip_analysis((os.path.join(root, name)),clipshape,(os.path.join(OutputFolder,name)))

End of code.

I also found the python glob method. It also can find all the tifs by searching through the tree like:

#Find the input feature classes

print (glob.glob("c:\\avdata\\PythonTest\\ClipLidar\\TestData\\Input\\*\\*.tif"),'\n')

##for root, dirs, files in os.walk(InputFolder):

## for name in files:

## if os.path.splitext(name)[1] == ".tif":

## print os.path.join(root, name), '\n'

where the glob statement returns the same list as the next four commented out statements. However, I eventually want to parse out the folders in which the files were found so I can reconstruct the input tree structure in the output tree. Don't know I can do that with the glob returns. Actually I don't know how to to that yet with the walk returns either.

ChrisSnyder

Some comments:

1. os.walk does get a you a handle on the directory name that is being parsed, so you could look for if dirName[-4:] == ".gdb" and if found, set the workspace to it, and issue arcpy.listfeatureclasses(). However, if I were you I wouldn't waste time developing code in anything < than v10.1 (since there is so many cool new toys in v10.1 SP1 like .da cursors, and .da.walk).

2. How about just intersecting the extent rectangle of your clip layer vs. the extent rectangles of the FCs that are returned by *.walk? If there is no intersection, don't clip. This extent comparison is easier/faster to do in v10.1 BTW.

PaulHuffman

It looks like I need to do it in 10.0 unless I travel there with a > 10.0 laptop installation. And I might have to do that because some of that data is in geodatabases, so os.walk will not find them.

Something else has come up. Many of the input feature classes are way outside the AOI clip. Too much time is being taken up by handling these out of range FCs and the output files are just empty. Is there a way of avoiding the out of range data in the clipping loop? It looks like I could first test for intersection with the clip polygon with Select layer by Location but I'm having a hard time telling if that tool is available at 10.0. I don't have 10.0 installed anywhere at my current location, and I was unable to find 10.0 help online.

ChrisBrannin__GISP

I want to loop through a folder and clip all feature classes by a AOI poly and write the results to an output folder. How to I also have the loop search all subfolders for vectors and raster feature classes as well?

This will walk all directories you want, including GDB's. You can also add the if to remove directories you're not interested in crawling.

for dirpath, dirnames, filenames in arcpy.da.Walk(in_workspace, datatype="FeatureClass",type="All"):
    if "Remove_Folder" in dirnames:
        dirnames.remove('Remove_Folder')
    for filename in filenames:
        if "{Field Name}" in [f.name for f in arcpy.ListFields(os.path.join(dirpath, filename))]:
            arcpy.Clip_analysis(os.path.join(dirpath, filename), clip_features, out_feature_class, xy_tolerance)

markdenil

Or, if you are using an older (<= 10.1) Arc,
you can use the python os.walk

for root, dirs, files in os.walk(basePathVariable):
        print "%s, " % (root)
        print str(sum(getsize(join(root, name)) for name in files))
        print " bytes, %s files\n" % (len(files))

prints the number and size of non-directory files in each directory and subdirectory.

Not exactly what you want to do, but you can get where you want to go from there...

ChrisSnyder

If you have v10.1+, you can use arcpy.da.walk:

http://resources.arcgis.com/en/help/main/10.1/index.html#//018w00000023000000

Rejoindre la conversation

Les membres connectés peuvent publier, suivre les mises à jour, et plus encore. Nouveau ici ? Inscrivez-vous gratuitement.

Liens rapides

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.