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.
UPDATED/WORKING CODE
# Subset Generator# Author: George Corea, Atherton Tablelands GIS# georgec@atgis.com.au; info@atgis.com.au# Licence:#AddMessage, AddWarning, AddErrorimport arcinfoimport arcpy, glob, os, sys, stringfrom arcpy import envfrom arcpy.sa import *path = os.getcwd()os.chdir(path)class LicenseError(Exception): passtry: 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 AOIextentSR = 0x=0Subset=0FileArea=0f = 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"
أكثر...
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).
# Subset Generator# Author: George Corea, Atherton Tablelands GIS# georgec@atgis.com.au; info@atgis.com.au# Licence:#AddMessage, AddWarning, AddErrorimport arcinfoimport arcpy, glob, os, sys, stringfrom arcpy import envfrom arcpy.sa import *path = os.getcwd()os.chdir(path)class LicenseError(Exception): passtry: 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 AOIextentSR = 0x=0Subset=0FileArea=0f = 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"
أكثر...