Comparing extents of rasters using ArcPy?

Question

I am writing a script that takes an input of an Area of Interest (AOI) and then walks through the base directories to find any file that has an extent that intersects with the AOI, clip it and generate a text output file. I've got most of it working except comparing the raster extent to the AOI extent. It is not giving errors but just doesn't work...

Any advise?

I am also having issues with trying to clip the extent to a polygon, it clips but maintains the MAX/MIN and doesn't clip to the mask. I have tried three methods as in the comments in the script.

Extract by Mask says I don't have Spatial Analyst but it is on.

enter image description here

ERROR MESSAGE when ExtractByMask is used.

Spatial Analyst license is AVAILABLE P:2011円\Job_154_PythonScript_for_AOI\Working\Orthophotomosaic310000円_8105000.tif 310000 8104000 311000 8105000 NaN NaN NaN NaN outside extent Traceback (most recent call last): File "P:2011円\Job_154_PythonScript_for_AOI\Working\SubsetGenerator.py", line 80, in outExtractByMask = ExtractByMask(File, AOI) File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 6726, in ExtractByMask in_mask_data) File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Utils.py", line 47, in swapper result = wrapper(*args, **kwargs) File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 6722, in wrapper out_raster) File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing_base.py", line 474, in return lambda *args: val(*gp_fixargs(args)) ExecuteError: Failed to execute. Parameters are not valid. ERROR 000824: The tool is not licensed. Failed to execute (ExtractByMask).

UPDATED/WORKING CODE

# Subset Generator
# Author: George Corea, Atherton Tablelands GIS
# [email protected]; [email protected]
# Licence:
#AddMessage, AddWarning, AddError
import arcinfo
import arcpy, glob, os, sys, string
from arcpy import env
from arcpy.sa import *
path = os.getcwd()
os.chdir(path)
class LicenseError(Exception):
 pass
try:
 if arcpy.CheckExtension("Spatial") == "Available":
 arcpy.CheckOutExtension("Spatial")
 print "Spatial Analyst license is AVAILABLE"
 else:
 # raise a custom exception
 #
 raise LicenseError
 env.workspace = path+'\\junk'
except LicenseError:
 print "Spatial Analyst license is unavailable"
except:
 print arcpy.GetMessages(2)
RootDirectories = [r'P:2011円\Job_154_PythonScript_for_AOI\Working\Orthophotomosaic']
RasterTypes = ['tif','jpg', 'ecw'] #boolean, valuelist
#AOI = arcpy.mapping.Layer('AOI.shp')
#AOIextent=AOI.getExtent()
AOI = 'AOI.tif'
AOIobj=arcpy.Raster(AOI)
AOIextent=AOIobj.extent
#print AOIextent
SR = 0
x=0
Subset=0
FileArea=0
f = open(path+'\\AOI_Clip\\AOI_SubsetGenerator.txt', 'a')
f.write("UniqueID, baseName, extension, dataType, path, \
OrigName, NewName, size(uncompressed bytes), area(m2)"+"\n")
f.close()
for RootDirectory in RootDirectories:
 for root, dirs, files in os.walk(RootDirectory):
 for RasterType in RasterTypes:
 FileList = [os.path.join(root, fi) for fi in files \
 if fi.endswith(RasterType)]
 for File in FileList:
 FileDesc=arcpy.Describe(File)
 SR = FileDesc.spatialReference
 if SR.name != "Unknown":
 FileObj = arcpy.Raster(File)
 FileExtent = FileObj.extent
 print "Processing: " + File
 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"
 else:
 #print "OUTSIDE extent"+str(Subset)
 Subset =0
#Process files that are subset
 if Subset == 1:
 #print Subset
 x=x+1
 outRaster=path+'\\AOI_Clip\\'+str(FileDesc.baseName)+"^AOI_Clip"+str(x)+"."+str(FileDesc.extension)
 try:
 outExtractByMask = ExtractByMask(File, 'AOI.shp')
 outExtractByMask.save(outRaster)
 print "Created: " + str(outRaster)
 except:
 print "Raster is outside of extent"
 #FileArea=float(arcpy.GetRasterProperties_management(File, "ROWCOUNT")*arcpy.GetRasterProperties_management(File, "COLUMNCOUNT")*arcpy.GetRasterProperties_management(File, "CELLSIZEY")*2)
 FileArea=\
 (string.atol(str(arcpy.GetRasterProperties_management(outRaster, "ROWCOUNT")))\
 *string.atol(str(arcpy.GetRasterProperties_management(outRaster, "COLUMNCOUNT"))))\
 *(string.atol((str(arcpy.GetRasterProperties_management(outRaster, "CELLSIZEY")))\
 *string.atol((str(arcpy.GetRasterProperties_management(outRaster, "CELLSIZEX"))))))
 print FileArea
 f = open(path+'\\AOI_Clip\\AOI_SubsetGenerator.txt', 'a')
 f.write(str(x)+","+str(FileDesc.baseName)+","\
 +str(FileDesc.extension)+","\
 +str(FileDesc.dataType)+","\
 +str(FileDesc.path)+","\
 +str(SR.name)+","\
 +str(outRaster)+","\
 +str(FileObj.uncompressedSize)+","\
 +str(FileArea)
 +"\n")
 f.close()
 else:
 print "No Subsets"
 print "Changing Directory to: " + str(dirs)
arcpy.CheckInExtension("Spatial")
print "Checked Out extensions and Exiting"

• user3461: Are the rasters GeoTiffs, or are they images with world reference files (e.g. each .TIF has corresponding .TFW file)? If they use world files, then you MAY be able to get extents by making calculations from the info in the world file and the dimensions (width and height in pixels) in the image file.
• GeorgeC: Thanks. They have TFW's. Calculating the extent isn't an issue. I can create a polygon around each tif. The issue is "how do I compare the extents?". So if the extent of the TIF crosses the AOI, it's reference is added to the csv and it is cut to the AOI polygon. The whole purpose is to automatically select any images that are entirely within or touches/crosses the AOI so that we can supply these to clients with only the AOI as an input.
• user3461: Sorry, I wasn't clear that I meant was to make the calculation results a vector that is the footprint of the image. Once that is done, an intersection of the footprint and AOI vectors is easy.
• GeorgeC: Thanks...i didn't want to go the route of processing the vector frootprint as arcpy seemed to have an extent comparison method but the logic of it was different from what I thought. I had to first save the result as a true/false and then if true run the other commands. Best.

Answer 1

Try using the ExtractByMask tool again. Many times, even though ArcMap or ArcCatalog show that you have the Spatial Analyst extension checked out, the script will not recognize it. You can explicitly check out the Spatial Analyst extension at the beginning of your script and then it should work.

Check out this link to ArcGIS help:

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/CheckOutExtension/000v0000003q000000/

• GeorgeC: Thanks Brian...I tried it but it still doesn't work. Pls see updated post above.
• GeorgeC: got it to work. Please see the completed script added at the end of the original question.