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.

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.

share|improve this question
    
I haven't figured out why Clip_management doesn't work - I suspect that the envelope has rounding errors despite being extracted from the image. However, a work-around is to use arcpy.sa.ExtractByMask(tempFilename, sourceImage). This works --no extra rows or columns! –  Peter Jul 23 '13 at 21:40
add comment

2 Answers

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.

Have you tried running the tool without using the optional parameter in_template_dataset?

Also, you might find this way easier to create the envelope. The numbers in the brackets refer to index numbers for the format function, and the values are auto casted to string. The spaces between brackets are important.

envelope = "{0} {1} {2} {3}".format(left, bottom, right, top)

I notice that you are creating a temporary folder for your projected image. You might be interested in using the in_memory workspace, which is designed to store intermediate steps of workflows. Assuming you have enough memory (RAM), this route is usually faster than saving a temp file to disk. Just remember to Delete afterwards!

share|improve this answer
    
Regarding the help on Clip: shouldn't the clip extent be aligned with the input raster because the input raster was created with snapRaster = in_template_dataset? I tried running Clip without the in_template_dataset -- it resulted in an extra column and extra row. Thanks for the tips on envelope creation and in_memory_workspace. I'll be sure to implement them. –  Peter Jul 19 '13 at 15:36
add comment

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.

share|improve this answer
add comment

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.