1

I'm new to ArcPy and to Python in general.

I have a huge raster dataset of all over the World (~1000 tiles .tif) and i'd like to automatically copy the tiles which are intersecting with 24 polygons (one shapefile containing polygons for 24 cities in Africa) in order to keep only the tiles that partially or completely contains my cities. I try using the script below. But when I do, he gives me all the tiles that are in between those cities (so pretty much all tiles in Africa). Probably because he uses the total extent of my shapefile, and not the polygons that are in that shapefile.

Does anyone knows how I could correct this so I get as result only the tiles that are intersecting with the 24 polygons inside that shapefile?

import os, sys, arcpy

AOI     = 
"C:/Users/jbihin/Desktop/testGufPython/To_Jeremie/Cities24_50km.shp"  # the 
area of interest to find the rasters for
rFolder = "D:/GUF/GUF04" # the folder containing the rasters to search from
oFolder = "C:/Users/jbihin/Desktop/testGufPython/To_Jeremie/gufSelect" # the         
folder to copy to

desc = arcpy.Describe(AOI) # get the extent of the AOI
sExt = desc.extent

arcpy.env.workspace = rFolder

for ThisRas in arcpy.ListRasters():
    rDesc = arcpy.Describe(ThisRas)
    rExt  = rDesc.extent
# check if this extent is related spatially
# by using not disjoint
    if sExt.disjoint(rExt):
        arcpy.AddMessage("Raster %s is outside" % (ThisRas))
else:
    arcpy.AddMessage("Raster %s overlaps" % (ThisRas))
    outFile = os.path.join(oFolder,ThisRas)
    arcpy.Copy_management(os.path.join(rFolder,ThisRas),outFile)
1

The issue comes from the fact that the extent of the AOI shapefile is the coordinate bounds of the ENTIRE area covered, as you mentioned, and is not restricted to the shapes of the individual features.

You'll likely want to add another loop for getting the AOI extents for each feature in your AOI shapefile, then using that list to compare within the FOR loop of your raster tiles.

Using your existing code (and assuming no issues with the data), that could look like this:

import os, sys, arcpy

AOI     = 
"C:/Users/jbihin/Desktop/testGufPython/To_Jeremie/Cities24_50km.shp"  # the 
area of interest to find the rasters for
rFolder = "D:/GUF/GUF04" # the folder containing the rasters to search from
oFolder = "C:/Users/jbihin/Desktop/testGufPython/To_Jeremie/gufSelect" # the         
folder to copy to

#Get the extent of each polygon feature individually from AOI, in a list
listExt = []
with arcpy.da.SearchCursor(AOI, 'SHAPE@') as cursor:
    for row in cursor:
        listExt.append(row[0].extent)

arcpy.env.workspace = rFolder

for ThisRas in arcpy.ListRasters():
    rDesc = arcpy.Describe(ThisRas)
    rExt  = rDesc.extent
    for sExt in listExt:
    # check if this extent is related spatially
    # by using not disjoint
        if sExt.disjoint(rExt):
            arcpy.AddMessage("Raster %s is outside" % (ThisRas))
        else:
            arcpy.AddMessage("Raster %s overlaps" % (ThisRas))
            outFile = os.path.join(oFolder,ThisRas)
            arcpy.Copy_management(os.path.join(rFolder,ThisRas),outFile)

Or, without Data Analysis tools, like this:

import os, sys, arcpy

AOI     = 
"C:/Users/jbihin/Desktop/testGufPython/To_Jeremie/Cities24_50km.shp"  # the 
area of interest to find the rasters for
rFolder = "D:/GUF/GUF04" # the folder containing the rasters to search from
oFolder = "C:/Users/jbihin/Desktop/testGufPython/To_Jeremie/gufSelect" # the         
folder to copy to

#Get the extent of each polygon feature individually from AOI, in a list
listExt = []
cursor = arcpy.SearchCursor(AOI)
for row in cursor:
    ext_row = row.getValue('Shape').extent  #NOTE:Geometry field name must match shapefile table   
    listExt.append(ext_row)
del cursor
arcpy.env.workspace = rFolder

for ThisRas in arcpy.ListRasters():
    rDesc = arcpy.Describe(ThisRas)
    rExt  = rDesc.extent
    for sExt in listExt:
    # check if this extent is related spatially
    # by using not disjoint
        if sExt.disjoint(rExt):
            arcpy.AddMessage("Raster %s is outside" % (ThisRas))
        else:
            arcpy.AddMessage("Raster %s overlaps" % (ThisRas))
            outFile = os.path.join(oFolder,ThisRas)
            arcpy.Copy_management(os.path.join(rFolder,ThisRas),outFile)
2
  • Thank you very much! I think this should solve the problem. Sadly, I'm still on arcgis 10.0, so the tool arcpy.da.SearchCursor is not working on my computer. Do you have any idea how to adapt the code with a "arcpy.SearchCursor" format? I've tried some different codes but I'm not comfortable yet with this topic and it does not work. – JBihin Mar 8 '18 at 9:45
  • See the updated solution, this should work. In the 'cursor' call, update the name "Shape" to match whatever the geometry field in your shapefile is (you can go into the layer properties to quickly see this). – AlecZ Mar 8 '18 at 17:55

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.