Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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. Thanks!

share|improve this question
    
have you tried converting the grid to points then adding X,Y fields and sorting? –  Jakub Aug 9 '11 at 15:15
    
Are the raster values floats or integers? –  whuber Aug 9 '11 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 Aug 9 '11 at 16:09
    
@whuber - The raster values are floats. –  mga Aug 9 '11 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. –  Jakub Aug 9 '11 at 16:58
show 1 more comment

3 Answers

up vote 5 down vote accepted

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.

share|improve this answer
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 Aug 18 '11 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 Aug 19 '11 at 14:52
add comment

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.

share|improve this answer
add comment

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
share|improve this answer
    
+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 Aug 10 '11 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 Aug 10 '11 at 23:17
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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