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.

Is there a way to determine what the difference in elevation is between two adjacent pixels? I am using an SRTM DEM to map faults. I am using the vertical error (>16m) as a criterion for the difference in elevation to accurately map these faults. There must be a more efficient way to do this than by using the i-tool.

Cheers

share|improve this question

3 Answers 3

up vote 6 down vote accepted

Compute the focal range grid using a 2 x 2 neighborhood. (Use the option where NoData cells are ignored.)

  • Any two adjacent pixels will be included within at least one neighborhood. Therefore, if any pair of cells differ by more than 16 m, they will cause at least one surrounding neighborhood to have a focal range exceeding 16 m.

  • If the focal range of a neighborhood exceeds 16 m, that means (by definition) there exist two cells differing by at least 16 m. Since the neighborhood is 2 x 2, those cells will be adjacent (perhaps diagonally).

This shows that you can identify all potential fault locations by selecting areas where the focal range exceeds 16 m.

If you don't want to include diagonally adjacent pixels, you can compute two focal range grids, one with a 1 x 2 neighborhood and the other with a 2 x 1 neighborhood, and then take their (local) maximum.

share|improve this answer
    
fantastic!!! Thanks heaps -- exactly what I was looking for! Cheers –  Kate Aug 23 '11 at 5:40

If you really need the discrete difference between adjacent row and/or column pixel values, you can use NumPy's diff function along one row or column axis:

import numpy as np
# Somehow load raster into NumPy, I use GDAL
>>> dtm = np.array([[560.78, 556.04, 559.11],
                    [559.84, 560.77, 560.24],
                    [560.03, 559.82, 560.05]])
>>> print(dtm)
[[ 560.78  556.04  559.11]
 [ 559.84  560.77  560.24]
 [ 560.03  559.82  560.05]]

>>> dzdrow = np.diff(dtm, axis=0)
>>> dzdrow
[[-0.94  4.73  1.13]
 [ 0.19 -0.95 -0.19]]

>>> dzdcol = np.diff(dtm, axis=1)
>>> dzdcol
[[-4.74  3.07]
 [ 0.93 -0.53]
 [-0.21  0.23]]

I generally would have more GDAL Python code here to write the two arrays to two raster files so that I can view them in ArcGIS or wherever.

As expected, the difference along columns output has one fewer column, and the difference along rows output has one fewer row. However, as there are two outputs from this analysis, it can be more difficult to interpret (compared to calculating slope with one output). With this method, the row and column directions are two perpendicular components, so you can't combine them into one raster output.

share|improve this answer
    
Thanks for the response but I don't know how to use python all that well (I've only ever used it for super basic applications). Plus my DEM is 1G... I don't think that would work. –  Kate Aug 22 '11 at 5:57
    
The raster size is not a problem (2GB is the limit where problems typically happen), but the Python-know-how is a barrier. I think @whuber's focal range idea is a good solution. –  Mike T Aug 22 '11 at 22:28

You could calculate the slope of the raster. Then based on the cell size, you could work out the slope angle which corresponded to a 16m difference in elevation.

share|improve this answer
1  
Thanks for the response... I have calculated slope and all the other derivative surfaces there are... Do you have a method I could use to show where exactly there is a difference of 16m+? Arc uses a 3x3 matrix to calculate the slope, with the center pixel not being used in the calculation. –  Kate Aug 21 '11 at 23:59
    
If you calculate the FLOWDIRECTION (help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//…) you'll know which cell was used to calculate the slope at that point. Based on that, you can calculate the elevation difference between the adjacent cells. Does that work? –  Stephen Lead Aug 22 '11 at 2:01
    
no...thanks though –  Kate Aug 22 '11 at 5:07
    
The slope doesn't work that way, Stephen. It's based on a least squares fit of a plane to the 3 x 3 neighborhood. Typically, it doesn't correspond to the slope between any pair of cells within the neighborhood. –  whuber Aug 22 '11 at 15:32
    
Fair enough - thanks for the clarification –  Stephen Lead Aug 22 '11 at 22:51

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.