[GIS] Clip_management to raster extent outputs larger raster than raster extent

alignmentarcgis-10.0arcpyraster

Using Python and ArcMap 10.0 SP5, I'm trying to reproject and clip a large input raster to exactly the same grid as a smaller target raster (same extent, pixel sizes and projection). The problem is that the script produces an output raster that has 1 more row than the target raster (e.g. 1324 rows instead of 1323 — columns match perfectly, pixels appear to be on the same grid). There is an extra row of pixels that overhangs the target image – everything else appears to be okay.

The following is my code for projecting and clipping to a target raster. What's needed to ensure that the output extent is exactly the same as in target raster? I'm new to ArcMap so it could be something obvious.

import arcpy
import tempfile
import os
import shutil

sourceImage = "C:/data/WorldClim_alt/WorldClim_Alt.tif"
inputImage = "E:/subset/nCom_2000/nCom2000101to107.tif"
outputFilename = "E:/subset/test/nComTest_clipped.tif"

# Snap to source image and source image extent
arcpy.env.extent = sourceImage
arcpy.env.snapRaster = sourceImage

# Get cell size
cellx = arcpy.GetRasterProperties_management(sourceImage, "CELLSIZEX")

# Get temporary directory and filename
tempDir = tempfile.mkdtemp()
tempFilename = 'E:/subset/test/' + basename(outputFilename) +'_proj.tif'

# Reproject
arcpy.ProjectRaster_management(inputImage,tempFilename,sourceImage, "#", cellx)

#
# Clip to source file
#

# Get envelope for clipping
top = arcpy.GetRasterProperties_management(sourceImage, "TOP")
left = arcpy.GetRasterProperties_management(sourceImage, "LEFT")
bottom = arcpy.GetRasterProperties_management(sourceImage, "BOTTOM")
right = arcpy.GetRasterProperties_management(sourceImage, "RIGHT")
envelope = str(left) + ' ' + str(bottom) + ' '+ str(right) + ' ' + str(top)

# Clip 
arcpy.Clip_management(tempFilename, envelope, outputFilename, sourceImage,"#","NONE")

Edit in Response to Paul's comment on Clip help:

From the help on Clip, it mentions

If the clip extent specified is not aligned with the input raster
dataset, this tool makes sure that the proper alignment is used. This
may cause the output to have a slightly different extent than
specified in the tool.

I think that the input to clip (the output from ProjectRaster) may be aligned with the in_template_dataset. The snapRaster environment variable was set to the in_template_dataset before ProjectRaster was called. Also, the output cell size of ProjectRaster was set to the in_template_dataset cell size. Visual inspection of the grids seems to indicate that they are aligned (except for an extra row and column in the ProjectRaster output).

The only effect of Clip seems to be the removal of a single column.

Best Answer

In my experience with clip management, these types of issues are often due to rounding errors. Basically it seems like actual image coordinates are often saved with greater precision than what is returned from, for instance, GetRasterProperties_management. In some tools, the output will be resampled based on the exact coordinates, but here, I believe it's just selecting input cells that are within the extent (the goodness or badness of this depends on what you're trying to do). For example, if the actual image coordinate of 'top' is

345109.3011115

And your clip extent specified is

345109.301112

an additional row will be selected, because the extent extends ever so slightly beyond that top row. To remedy this, I often select an extent that is a fraction of a cell inside of the actual clipping extent (for top and right, that's coordinates slightly smaller, while for bottom and left it's coordinates slightly larger). In your example (assuming the cells are greater than 1 unit resolution) I would do something like

top = float(arcpy.GetRasterProperties_management(sourceImage, "TOP")) - 1.0
left = float(arcpy.GetRasterProperties_management(sourceImage, "LEFT")) + 1.0
bottom = float(arcpy.GetRasterProperties_management(sourceImage, "BOTTOM")) + 1.0
right = float(arcpy.GetRasterProperties_management(sourceImage, "RIGHT")) - 1.0

Good luck, ArcGIS loves to play games with these types of things.

Related Solutions

[GIS] When a raster is clipped to a polygon with Python, why does it still take on the extent of the entire shapefile

cross-posting this solution to this thread AND user3397's other thread

@lpinner was on the right track, by stating that the extent being used is the entire extent of the 'clipto' shapefile. This is an underlying attribute of that shapefile, and doesn't change according to the current feature being processed by a cursor. So, every time you're feeding the extent object into the Clip tool, you're giving it this same overall extent, over and over.

You can manipulate the output by changing the Extent environment setting - this is done by setting the gp.Extent property. In my experience, there are two different ways of doing this:

the legacy method involved getting and setting the extent as a space-delimited string,
the newer (9.3+ I think) method involves instantiating an 'Extent' object, and getting and setting its Xmax, YMax, etc properties

Option 1: legacy method

# Set extent property with a space-delimited string
# extents should be in this order: XMin YMin XMax YMax 
gp.Extent = "-180 -90 180 90"

Option 2: modern method

# Create Extent object, and set its attributes
ext = gp.CreateObject('Extent')
ext.XMin, ext.YMin, ext.XMax, ext.YMax = -180, -90, 180, 90
gp.Extent = ext

I'm not sure which ArcGIS versions restrict you to which method (off the top of my head, I think you can do both in 9.3+, but you can't do Option 2 in 9.2-).

So, to solve your problem, you need to extract the XMin, YMin, XMax and YMax lat/longs from each polygon, and change gp.Extent every time. I.e.:

while row:
    # do your other processing and extract extent co-ords from current polygon        
    ext = gp.CreateObject('Extent')
    ext.XMin, ext.YMin, ext.XMax, ext.YMax = currentXMin, currentYMin, currentXMax, currentYMax
    gp.Extent = ext
    gp.Clip_management #...continue with Clip operation

A third option, I guess, would be to write each polygon to a new temporary feature class. This temporary FC's Extent property would then be the same as the individual polygon's extent. This would be a slow operation if you wrote these temp files to disk, but you can write them to ArcGIS's 'in_memory' workspace to keep it nippy. I.e.:

# Get the primary key field name for the polygon shapefile
oidField = gp.Describe(clipit).OIDFieldName

while row: # Iterate through each polygon
    oid = row.GetValue(oidField) # get the unique id for this polygon

    # Use Select tool to extract individual polygon to temporary in_memory file
    # use the unique ID as the SQL query to extract this record
    gp.Select_management(clipit, 'in_memory/temp', '"%s" = %s' % (oidField, oid))  

    # get the polygon's extent
    ext = gp.Describe('in_memory/temp').Extent

From here, you can either use this in_memory/temp polygon as ClippingGeometry for the raster Clip tool, or you can extract the Extent object's values and feed them into the Clip tool's 'extent' parameters as text.

I'll cut my answer off here, and wait and see if it helps. I can clarify any points some more my referring to my scripts at work - I'm writing this off the top off my head at the moment.

[GIS] How to clip a raster by extent using Processing from the Python console

You can get the correct syntax from the Processing | History and log after running the tool once from the toolbox:

processing.runalg("gdalogr:cliprasterbyextent",   
                  "C:/Users/anita/Geodata/qgis_sample_data/raster/SR_50M_alaska_nad.tif",
                  "",
                  "2001518.84322,3027531.34481,4007380.09563,4664339.34638",
                  "",
                  None)

enter image description here

This calls gdal_translate like this:

gdal_translate -of GTiff -projwin -6232946.6727 9275122.96868 6363148.43764 -735684.661767 
C:/Users/anita/Geodata/qgis_sample_data/raster/SR_50M_alaska_nad.tif "[temporary file]"

Best Answer

Related Solutions

[GIS] When a raster is clipped to a polygon with Python, why does it still take on the extent of the entire shapefile

[GIS] How to clip a raster by extent using Processing from the Python console

Related Question