This is a variant of finding max elevations in a raster and converting them to point containing the max value. I have done this successfully using the full raster (39000x39000), zonal stats, and a CON expression. However, I want to use a small rectangle to define the search area, returning the max value and converting that to a point.

To speed processing, I first extract the portion of the raster delimited by master_mask (rectExtract). I then set the extent to that smaller raster and extract the maximum value within it using zonal statistics and store to a new raster (maxpixel). All proceeds fine to this point, both rasters having the same extent. (Alternatively, I skip this step and simply use rectExtract.maximum).

Then I use a CON expression to compare rectExtract with maxpixel. This is where the problem is. The output raster, instead of being the same extent as the inputs, is the extent of the full raster (39000x39000) and no match value is found, even though I explicitly set the extent to the size of the clipped raster beforehand. The convert-to-point then fails.

Any ideas?

import arcpy
from arcpy import env  
from arcpy.sa import *
# Check out the ArcGIS Spatial Analyst extension license


# Set the extent environment using a shapefile
arcpy.env.extent = master_mask
print arcpy.env.extent
arcpy.env.snapRaster = master_raster
env.workspace = "C:/sapyexamples/output"

# Execute ExtractByRectangle
rectExtract = ExtractByRectangle(master_raster, master_mask, "INSIDE")
print "Max of extracted rectangle"
print rectExtract.maximum,rectExtract.height,rectExtract.width
arcpy.env.snapRaster = rectExtract

maxpixel = ZonalStatistics(master_mask, "FID", rectExtract, "MAXIMUM", "NODATA")
print "below is the max pixel height"
print "maxpixel",maxpixel.maximum,maxpixel.height,maxpixel.width

arcpy.env.extent = master_mask
print "Extents before CON"
print arcpy.env.extent

#Problem is here, env.extent not being enforced
OutRas2 = Con(rectExtract == maxpixel,maxpixel)
print "OutRas2",OutRas2.maximum,OutRas2.height,OutRas2.width
#OutRas2 = Con(rectExtract == rectmax, rectmax)  (gives the same problem)
print "Finished peak pixel IDs"     

    arcpy.RasterToPoint_conversion(OutRas2, "Temp_peak", "VALUE")
    print "finished raster to point conversion"

    print "job complete"
  • 1
    Just an idea replace this line arcpy.env.extent = master_mask just before the print "Extents before CON" statement with arcpy.env.extent = rectExtract? – Hornbydd Jul 7 at 17:34
  • Thanks for replying, but the same problem. – William Jul 8 at 2:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.