Finding locations of highest values in raster using ArcGIS Desktop?

Question

Using ArcGIS 10, I have a raster where I would like to find the pixel with the maximum value in the raster and return its location (center of pixel) in decimal degrees. I would like to iterate through this process returning the location of the second highest value of the raster, then third, and so on so that in the end I have a list of N locations that have the highest values in the raster in order.

I imagine that this might be most easily done using a Python script but I'm open to other ideas if there is a better way.

Answer 1

I did this some time ago, although my solution is using GDAL (so, this is not just for ArcGIS). I think you can get a NumPy array from a raster in ArcGIS 10, but I don't know for sure. NumPy supplies simple and powerful array indexing, like argsort and others. This example doesn't handle NODATA or transform coordinates from projected to lat/long (but this isn't difficult to do with osgeo.osr, provided with GDAL)

import numpy as np
from osgeo import gdal

# Open raster file, and get GeoTransform
rast_src = gdal.Open(rast_fname)
rast_gt = rast_src.GetGeoTransform()

def get_xy(r, c):
    '''Get (x, y) raster centre coordinate at row, column'''
    x0, dx, rx, y0, ry, dy = rast_gt
    return(x0 + r*dx + dx/2.0, y0 + c*dy + dy/2.0)

# Get first raster band
rast_band = rast_src.GetRasterBand(1)

# Retrieve as NumPy array to do the serious work
rast = rast_band.ReadAsArray()

# Sort raster pixels from highest to lowest
sorted_ind = rast.argsort(axis=None)[::-1]

# Show highest top 10 values
for ind in sorted_ind[:10]:
    # Get row, column for index
    r, c = np.unravel_index(ind, rast.shape)
    # Get [projected] X and Y coordinates
    x, y = get_xy(r, c)
    print('[%3i, %3i] (%.3f, %.3f) = %.3f'%
          (r, c, x, y, rast[r, c]))

Shows the following for my test raster file:

[467, 169] (2813700.000, 6353100.000) = 844.538
[467, 168] (2813700.000, 6353200.000) = 841.067
[469, 168] (2813900.000, 6353200.000) = 840.705
[468, 168] (2813800.000, 6353200.000) = 840.192
[470, 167] (2814000.000, 6353300.000) = 837.063
[468, 169] (2813800.000, 6353100.000) = 837.063
[482, 166] (2815200.000, 6353400.000) = 833.038
[469, 167] (2813900.000, 6353300.000) = 832.825
[451, 181] (2812100.000, 6351900.000) = 828.064
[469, 169] (2813900.000, 6353100.000) = 827.514

Answer 2

The answer can be obtained by combining an indicator grid of the top 1% of values with grids of latitude and longitude. The trick lies in creating this indicator grid, because ArcGIS (still! after 40 years!) does not have a procedure to rank raster data.

One solution for floating-point rasters is iterative, but mercifully quick. Let n be the number of data cells. The empirical cumulative distribution of values consists of all pairs (z, n(z)) where z is a value in the grid and n(z) is the number of cells in the grid with values less than or equal to z. We get a curve connecting (-infinity, 0) to (+infinity, n) out of the sequence of these vertices ordered by z. It thereby defines a function f, where (z, f(z)) always lies on the curve. You want to find a point (z0, 0.99*n) on this curve.

In other words, the task is to find a zero of f(z) - (1-0.01)*n. Do this with any zero-finding routine (that can handle arbitrary functions: this one is not differentiable). The simplest, which is often efficient, is guess-and-check: initially you know z0 lies between the minimum value zMin and the maximum zMax. Guess any reasonable value strictly between these two. If the guess is too low, set zMin = z0; otherwise set zMax = z0. Now repeat. You will quickly converge to the solution; you're close enough when zMax and zMin are sufficiently close. To be conservative, choose the final value of zMin as the solution: it might pick up a few extra points that you can discard later. For more sophisticated approaches, see Chapter 9 of Numerical Recipes (the link goes to an older free version).

Looking back at this algorithm reveals you need to perform only two kinds of raster operations: (1) select all cells less than or equal to some target value and (2) count selected cells. Those are among the simplest and fastest operations around. (2) can be obtained as a zonal count or by reading one record from the attribute table of the selection grid.

Answer 3

If you are happy to use R, there is a package called raster. You can read in a raster using the following command:

install.packages('raster')
library(raster)
test <- raster('F:/myraster')

Then, when you go to look at it (by typing test), you can see the following info:

class       : RasterLayer 
dimensions  : 494, 427, 210938  (nrow, ncol, ncell)
resolution  : 200, 200  (x, y)
extent      : 1022155, 1107555, 1220237, 1319037  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 
values      : F:/myraster 
min value   : 0 
max value   : 1 

There may be better ways of manipulating the raster, but one way of finding the info you want could be to find the highest value, and get it's matrix location, and then add that to the lower extents.