Zonal Peak

Here is a short Python recipe for locating the highest point (peak) in an elevation grid or raster for a set of zones:

import arcpy
from arcpy.sa import *

def main(zones, zone_field, depth_raster, output_points):
    """Find the highest points for each raster zone."""
    arcpy.env.overwriteOutput = True

    # Get highest (max) value for each zone.
    arcpy.env.cellSize = depth_raster
    arcpy.env.snapRaster = depth_raster
    zonal_raster = ZonalStatistics(zones, zone_field, depth_raster, "MAXIMUM")

    # Retrieve all cells in depth raster equal to zonal raster values (max values).
    con_raster = Con(Raster(depth_raster) == zonal_raster, depth_raster)

    # Convert the raster cells to point features - output the highest locations.
    arcpy.conversion.RasterToPoint(con_raster,output_points,"Value")
# End main function

if __name__ == '__main__':
    # Retrieve input parameter values and run main.
    in_zones = arcpy.GetParameterAsText(0)
    zone_id_field = arcpy.GetParameterAsText(1)
    elev_raster = arcpy.GetParameterAsText(2)
    out_points = arcpy.GetParameterAsText(3)
    main(in_zones, zone_id_field, elev_raster, out_points)

To find the lowest (deepest) point in each zone, change the keyword “MAXIMUM” to “MINIMUM”.

About these ads

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s