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.
ERROR MESSAGE when ExtractByMask is used.
Spatial Analyst license is AVAILABLE P:\2011\Job_154_PythonScript_for_AOI\Working\Orthophotomosaic\310000_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"
