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.

  • have you tried converting the grid to points then adding X,Y fields and sorting? Commented Aug 9, 2011 at 15:15
  • Are the raster values floats or integers?
    – whuber
    Commented Aug 9, 2011 at 15:49
  • @Jakub - No I haven't. I'm probably only going to be interested in the top 1% or so of points so I don't know if it's worth adding x,y fields for all points then sorting. Maybe if there's not a more efficient option?
    – mga
    Commented Aug 9, 2011 at 16:09
  • @whuber - The raster values are floats.
    – mga
    Commented Aug 9, 2011 at 16:13
  • @mga, its worth a try. The conversion is quite fast and adding the XY is also a default tool. Deleting unwanted records is a straight forward operation and all could be bundled into a single model. Just an idea. Commented Aug 9, 2011 at 16:58

3 Answers 3


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
  • +1 Thanks for sharing this. I don't see an inability to handle NoData as a limitation: just convert all NoData to extremely negative values before proceeding. Note, too, that if any reprojection of the grid occurs, the answers will likely change due to resampling of the grid, so normally one does not want reprojection to happen automatically during a calculation like this. Instead, the reported coordinates can be reprojected afterwards. Thus your solution is perfectly general.
    – whuber
    Commented Aug 10, 2011 at 15:08
  • The NODATA handling can be implemented by first getting the value from the raster NODATA = rast_band.GetNoDataValue(), then by using either a NaN value (rast[rast == NODATA] = np.nan) or by using a masked array (rast = np.ma.array(rast, mask=(rast == NODATA))). The more complicated trick is getting argsort to somehow remove the NODATA values from analysis, or simply just skip past them in the for-loop if they are NaN/masked.
    – Mike T
    Commented Aug 10, 2011 at 23:17

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.


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

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.

  • 1
    +1 For this reference. Once you have read the raster into R, you can use standard R functions or the getValues method to access cell values. From there it's straightforward to identify the highest values and their locations.
    – whuber
    Commented Aug 18, 2011 at 16:52
  • 1
    Thanks to your recommendation this is what I have ended up doing. Using the raster package in R has been a breeze compared to trying this in ArcGIS. I've also continued on using other spatial analyses in R and have been very happy with the results. Great advice!
    – mga
    Commented Aug 19, 2011 at 14:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.