Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

Thanks guys,

enter image description here

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
# georgec@atgis.com.au; info@atgis.com.au
# 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"
share|improve this question
    
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. –  user3461 Aug 15 '11 at 13:29
    
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. –  GeorgeC Aug 15 '11 at 23:08
    
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. –  user3461 Aug 16 '11 at 17:18
    
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. –  GeorgeC Aug 17 '11 at 5:55

1 Answer 1

up vote 1 down vote accepted

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/

share|improve this answer
    
Thanks Brian...I tried it but it still doesn't work. Pls see updated post above. –  GeorgeC Aug 15 '11 at 23:01
    
got it to work. Please see the completed script added at the end of the original question. –  GeorgeC Aug 16 '11 at 5:37

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.