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.
import arcpy from arcpy import env from arcpy.sa import * # Check out the ArcGIS Spatial Analyst extension license arcpy.CheckOutExtension("Spatial") master_raster="NCALM_all.tif" master_mask="Extract_peak_box" # 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 rectmax=rectExtract.maximum 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" try: arcpy.RasterToPoint_conversion(OutRas2, "Temp_peak", "VALUE") print "finished raster to point conversion" finally: print "job complete"