Topographic position index (TPI) is a method of terrain classification where the altitude of each data point is evaluated against its neighbours. If a point is higher than its surroundings, the index will be positive, as for example on ridges and hilltops, while the figure will be negative for sunken features such as valleys. Note that the size of the examined neighbourhood matters; what may globally be a large valley can locally become flat terrain or even elevated relief. Weiss (2001) has proposed a complex terrain classification on the basis of the amount of difference (hilltop, incised valley, shallow valley etc.).

TPI_diagram.png Topographic position index depends on neighbourhood size (Salinas-Melgoza et al. 2018, Fig. 2)

Technically, TPI functions as simple average filter or, more precisely, a kernel filter where each cell is assigned the average value of its neighbourhood (see on wikipedia). The neighbourhood size and shape is typically fixed, it does not vary from cell to cell. This average value is then subtracted from the actual height value to evaluate the difference (negative values will, then, represent concave surfaces and positive ones convexities).

tpi Terrain position index (concavities are in blue and convexities in red). These are animal tracks and some human paths on a Lidar derived DEM.

Today, the TPI is included in most of GIS software, such as QGIS, GRASS, ArcGis or Saga – so why should we bother? The answer is in the simplicity of the method: it can easily be customised and experimented with in order to obtain that original result which will distinguish your work form other push-button solutions. We can play with window size or specific weight factors, for instance when window periphery should count less than the window centre. Such freedom can be very useful when fine tuning Lidar visualisations, for example.

In what follows, I will only provide the basic architecture for a TPI filter, based on the sliding window routine developed for numpy. Window size and weights can be fine-tuned in relation to specific problems (see below). To test the script you can use QGIS (open the Python console and just copy-paste the code) or run in pure Python.

console Running the script in QGIS Python console

"""
Topographic position index for elevation models, 
a mock script to be tuned according to you needs.

All comments are welcome at LandscapeArchaeology.org/2019/tpi
by Zoran Čučković
"""

import gdal
import numpy as np

# -------------- INPUT -----------------
elevation_model ="__path__to__my__dem__"
output_model = "__path__to__my__output__file__"


# ----------  create the moving window  ------------
r= 5 #radius in pixels
win = np.ones((2* r +1, 2* r +1))
# ----------   or, copy paste your window matrix -------------
# win = np.array( [    [0, 1, 1, 1, 0]
#                      [1, 1, 1, 1, 1],
#                      [1, 1, 0, 1, 1],
#                      [1, 1, 1, 1, 1],
#                      [0, 1, 1, 1, 0]  ])

# window radius is needed for the function,
# deduce from window size (can be different for height and width…)
r_y, r_x  = win.shape[0]//2, win.shape[1]//2
win[r_y, r_x  ]=0  # let's remove the central cell 

def view (offset_y, offset_x, shape, step=1):
    """
    Function returning two matching numpy views for moving window routines.
    - 'offset_y' and 'offset_x' refer to the shift in relation to the analysed (central) cell 
    - 'shape' are 2 dimensions of the data matrix (not of the window!)
    - 'view_in' is the shifted view and 'view_out' is the position of central cells
    (see on LandscapeArchaeology.org/2018/numpy-loops/)
    """
    size_y, size_x = shape
    x, y = abs(offset_x), abs(offset_y)
    
    x_in = slice(x , size_x, step) 
    x_out = slice(0, size_x - x, step)

    y_in = slice(y, size_y, step)
    y_out = slice(0, size_y - y, step)
 
    # the swapping trick    
    if offset_x < 0: x_in, x_out = x_out, x_in                                 
    if offset_y < 0: y_in, y_out = y_out, y_in
 
    # return window view (in) and main view (out)
    return np.s_[y_in, x_in], np.s_[y_out, x_out]

# ----  main routine  -------

dem = gdal.Open(elevation_model)
mx_z = dem.ReadAsArray()

#matrices for temporary data
mx_temp = np.zeros(mx_z.shape)
mx_count = np.zeros(mx_z.shape)

# loop through window and accumulate values
for (y,x), weight in np.ndenumerate(win):
    
    if weight == 0 : continue  #skip zero values !
    # determine views to extract data 
    view_in, view_out = view(y - r_y, x - r_x, mx_z.shape)
    # using window weights (eg. for a Gaussian function)
    mx_temp[view_out] += mx_z[view_in]  * weight
    
   # track the number of neighbours 
   # (this is used for weighted mean : Σ weights*val / Σ weights)
    mx_count[view_out] += weight

# this is TPI (spot height – average neighbourhood height)
out = mx_z - mx_temp / mx_count

# writing output 
driver = gdal.GetDriverByName('GTiff')
ds = driver.Create(output_model, mx_z.shape[1], mx_z.shape[0], 1, gdal.GDT_Float32)
ds.SetProjection(dem.GetProjection())
ds.SetGeoTransform(dem.GetGeoTransform())
ds.GetRasterBand(1).WriteArray(out)
ds = None

Sliding window

The script is using both, window geometry (square, circle etc.) and window weights which are applied to individual values. Geometry is modelled by setting unwanted window parts to zero, while weights should be non-zero values. Note that zero values are simply skipped as if non-existent, which will not happen with, say, weight = 0.00001.

It is not difficult to create a window on the fly:

r = 5 #radius in pixels
win = np.ones((2*r +1, 2*r +1)) #square window, all values are 1
win[r, r] = 0  #skipping the central cell

For Gaussian window weights: see numpy documentation.

However, I find more useful (and intuitive) to hard-code window definitions. This is very handy when comparing results on different datasets – there can be no error when the window is specified cell by cell. We can also keep a diary with window definitions that work best, like:

  • Fine 5x5 filter :

                  [0, 0, 0.5, 0, 0],
                  [0, 0.5, 1, 0.5, 0],
                  [0.5, 1, 0, 1, 0.5],
                  [0, 0.5, 1, 0.5, 0],
                  [0, 0, 0.5, 0, 0] 
    
  • Aggressive 5x5 filter (best on such and such DEMs…)

                  [0, 1, 1, 1, 0],
                  [1, 1, 1, 1, 1],
                  [1, 1, 0, 1, 1],
                  [1, 1, 1, 1, 1],
                  [0, 1, 1, 1, 0]
    
  • A totally skewed filter for detecting linear patterns:

                  [0, 0, 0, 0.5, 1],
                  [0, 0, 0.5, 1, 0.5],
                  [0, 0.5, 0, 0.5, 0],
                  [0.5, 1, 0.5, 0, 0],
                  [1, 0.5, 0, 0, 0 ]
    

Danger in this approach is in the assumption that the central data point is in the geometric centre of the kernel window. For this to be true, dimensions of the window should be in odd numbers. For instance, in 3x3 window, the central pixel is at (1,1) position. In a 4x4 window, the centre falls between cells, which doesn’t work! Therefore, be sure to count your pixels.

Larger windows may become unwildely; in that case you can always use the on-the-fly method, as above.

Happy mapping !

Bibliography

Kernel processing on Wikipedia

A. Weiss (2001) : Topographic Position and Landforms Analysis, Poster Presentation, ESRI Users Conference, San Diego, CA

J. Jeness (2006): Topographic Position Index (TPI) v. 1.2

M. A. Salinas-Melgoza, M. Skutsch, and J. C. Lovett (2018): Predicting aboveground forest biomass with topographic variables in human-impacted tropical dry forest landscapes. Ecosphere 9(1):e02063. 10.1002/ecs2.2063