Comparing Raster Extents

JonathanMulder · ‎11-28-2013

Greetings,

I'm very new to Python having just taken ESRI's 3-day Python scripting course.

I'm trying to clip a large raster into smaller rasters that are the size of DOQQs. I found some code that compares raster extents, then if the extents intersect, it clips the raster.

My code has issues with this line: "If DOQQExtent.within(TargetRasterExtent.extent):" This is my complete code:

import arcpy

# Get directories of DOQQ images
in_raster = arcpy.GetParameterAsText(0)
in_directories = arcpy.GetParameterAsText(1)
out_raster = arcpy.GetParameterAsText(2)
TargetRasterExtent = arcpy.Describe(in_raster).extent
for CurrentDirectory in in_directories:
    DOQQRasters = arcpy.ListFiles("*.tif")
    for CurrentRaster in CurrentDirectory:
        DOQQExtent = arcpy.Describe(CurrentRaster).extent
        If DOQQExtent.within(TargetRasterExtent.extent):
#        if arcpy.Describe(CurrentRaster).extent.overlaps(TargetRasterExtent) == 'True':
            print "Extent Overlaps"

Any idea what's wrong with the line?

Thanks for any help and Happy Thanksgiving!

Jonathan Mulder
California Departement of Water Resources

JonathanMulder · ‎11-28-2013

Update.

I now realize I have to check out the Spatial Analyst before comparing extents. So I modified my code as follows:

import arcpy
arcpy.CheckOutExtension("spatial")

# Get directories of DOQQ images
in_raster = arcpy.GetParameterAsText(0)
in_directories = arcpy.GetParameterAsText(1)
out_raster = arcpy.GetParameterAsText(2)
TargetRasterExtent = arcpy.Describe(in_raster).extent
for CurrentDirectory in in_directories:
    DOQQRasters = arcpy.ListFiles("*.tif")
    for CurrentRaster in CurrentDirectory:
        DOQQExtent = arcpy.Describe(CurrentRaster).extent
        if DOQQExtent.within(TargetRasterExtent.extent):
#        if arcpy.Describe(CurrentRaster).extent.overlaps(TargetRasterExtent) == 'True':
            print "Extent Overlaps"

Unfortunately, it still doesn't work and returns the following error:

Executing: Script H:\Documents\GIS\BaseData\Raster\Imagery\DOQQs\39121\o39121b5nw_YubaCity.tif H:\Documents\GIS\BaseData\Raster\Imagery\DOQQs\39121 H:\Documents\GIS\BaseData\Raster\Imagery\DOQQs\39121
Start Time: Thu Nov 28 14:11:24 2013
Running script Script...

Traceback (most recent call last):
File "H:\Documents\GIS\PythonScripts\CutRasterByDOQQs.py", line 25, in <module>
    DOQQExtent = arcpy.Describe(CurrentRaster).extent
File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\__init__.py", line 1234, in Describe
    return gp.describe(value)
File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\geoprocessing\_base.py", line 374, in describe
    self._gp.Describe(*gp_fixargs(args, True)))
IOError: "H" does not exist

Failed to execute (Script).
Failed at Thu Nov 28 14:11:25 2013 (Elapsed Time: 0.84 seconds)

JonathanMulder · ‎11-28-2013

Update: I think I'm getting closer but still not getting the clips. I've modified my code as follows, but no clipped rasters are generated. I suspect it's in my line comparing extents and the Clip_management function doesn't fire off.

Any insight?

import arcpy
#from arcpy import env
#from arcpy.sa import*
arcpy.CheckOutExtension("spatial")

# Get name of raster mosaic
in_rastermosaic = arcpy.GetParameterAsText(0)
# Get directories of DOQQ images
in_directories = arcpy.GetParameterAsText(1)
# Get directories of outputted rasters
out_raster = arcpy.GetParameterAsText(2)

TargetRasterExtent = arcpy.Describe(in_rastermosaic).extent
for CurrentDirectory in in_directories:
    DOQQRasters = arcpy.ListFiles("*.tif")
    for CurrentRaster in DOQQRasters:
        DOQQExtent = arcpy.Describe(CurrentRaster).extent
        if DOQQExtent.within(TargetRasterExtent.extent):
            arcpy.Clip_management(in_rastermosaic,"#","1958_" + CurrentRaster,"#",0,"NONE")
    #        if arcpy.Describe(CurrentRaster).extent.overlaps(TargetRasterExtent) == 'True':
            print "Extent Overlaps"
        else:
            print "Extent does not Overlaps"