[GIS] Comparing extents of rasters using ArcPy

arcgis-10.0arcgis-desktoparcpyerror-000824spatial-analyst

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\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"

Best Answer

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/

Related Solutions

ArcPy – Inserting DATAFRAME_ELEMENT Value into Dynamically Generated Shapefile

I'm a bit lost here, so my answer may be off base. When you say you want the ExtentID, that appears to be the name of the Dataframe you're exporting the extent of. If this is the case (which is what your code is doing to my best guess), you don't need the [x] index when printing it. :

 for el in arcpy.mapping.ListLayoutElements(mxd2, "DATAFRAME_ELEMENT"):
     ExtentIDStore = str(el.name)
     #print x, y, z <--- You don't need to print x y and z. They are unrelated at this point.
     print ExtentIDStore #used to see which item to add to the shape file

     print ExtentIDStore[x] # <- This will just print an indexed character in the string

By the time you reach the loop above, your code is only concerned with the Dataframe Element in question, so there is no need to do any indexing for storage. ExtendIDStore should be what you add to the shapefile row.

Your error which results at this line:

row.strField = str(x)

Is probably because you are not referencing a field that exists. What comes after the row reference is the name of the Field. In this case, it should be:

row.ExtentID = ExtentIDStore

File name you can get via the filePath property of the MapDocument class and the os.path.basename method, i.e.:

filepath = mxd2.filePath # Gets full path of MXD used in your dataframe loop
filename = os.path.basename(filepath) # Extracts the file name and extension
row.FileRef = filename

If you don't want the extension, use filename.strip(".mxd") to remove it.

Scale would be:

scale = el.scale # Scale of the dataframe referenced in the loop
row.Scale = scale

And that should fill in the rows you're looking for.

For the new problem, I think you are misusing the insert cursor. The second argument should be a spatial reference, which cannot be gleaned via an MXD. I would recommend creating the cursor with only one argument (the output file or dfStore). You also forgot the last (and most important) line:

rows.insertRow(row)

Without that your data never gets put into the table. Put that down at the end of the row loop so that the data is actually inserted after all the variables are set.

For more information about accessing Dataframe, layer, and map properties (scale, name, current position), check out the mapping module. The scale property of the Dataframe class can be seen on the classes page of the same name.

Python GDAL – Getting Extent of Georeferenced Rasters and Output to Polygon Shapefile

You almost had it right with your first attempt. What you are doing is trying to call the extent property with a filename as its parameter, when you need to construct a Raster object with that parameter.

In practical terms, that means:

extent1 = arcgisscripting.Raster('stg1_05.jpg').extent

Although it is usually better practice to break it down into two steps:

raster1 = arcgisscripting.Raster('stg1_05.jpg')
extent1 = raster1.extent

Note the lack of parentheses on extent, this is because it is a property rather than a method.

Best Answer

Related Solutions

ArcPy – Inserting DATAFRAME_ELEMENT Value into Dynamically Generated Shapefile

Python GDAL – Getting Extent of Georeferenced Rasters and Output to Polygon Shapefile

Related Question