{"id":29,"date":"2021-12-01T08:00:00","date_gmt":"2021-12-01T08:00:00","guid":{"rendered":"http:\/\/wcwp.rf.gd\/2021\/12\/01\/https-landscapearchaeology-org-2021-python-tpi\/"},"modified":"2021-12-01T08:00:00","modified_gmt":"2021-12-01T08:00:00","slug":"python-tpi","status":"publish","type":"post","link":"https:\/\/landscapeanalysis.org\/?p=29","title":{"rendered":"Topographic position index (TPI): a Python implementation"},"content":{"rendered":"

Topographic position index (TPI) is a method of terrain classification where the altitude of each data point is evaluated against its neighbourhood. In a nutshell, for each data point, i.e. a pixel in raster DEM, we calculate its height difference from the immediate neighbourhood. A detailed explanation can be found in the post on the TPI method<\/a>. Here, I present a Python implementation of the TPI, for all those who wish to venture into custom filter programming.<\/p>\n

\"tpi\"
\nTerrain position index (concavities are in blue and convexities in red). These are animal tracks and some human paths on a Lidar derived DEM.<\/em><\/p>\n

The TPI method 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.<\/p>\n

In what follows, I will only provide the basic architecture for a TPI filter, based on the sliding window routine developed for numpy<\/a>. In essence, TPI is just a kernel routine, which gathers information around each data point. At the same time it applies a weight factor, for instance to correct for the distance difference between diagonal pixels (touching by a corner) and orthogonal ones (touching by a side; below).<\/p>\n

\"2021-12-filter_standard.jpg\"
\nStandard TPI kernel (moving window) with distance correction<\/em><\/p>\n

Below is a script that can be used in QGIS (open the Python console and just copy-paste the code) or run in pure Python.<\/p>\n

\"console\"
\nRunning the script in QGIS Python console<\/em><\/p>\n

\n
\n
"""\nTopographic position index for elevation models, \na mock script to be tuned according to you needs.\n\nAll comments are welcome at landscapeanalysis.org\/2019\/tpi\nby Zoran \u010cu\u010dkovi\u0107\n"""<\/span>\n\nimport<\/span> gdal<\/span>\nimport<\/span> numpy<\/span> as<\/span> np<\/span>\n\n# -------------- INPUT -----------------\n<\/span>elevation_model<\/span> =<\/span>\"__path__to__my__dem__\"<\/span>\noutput_model<\/span> =<\/span> \"__path__to__my__output__file__\"<\/span>\n\n\n# ----------  create the moving window  ------------\n<\/span>r<\/span>=<\/span> 5<\/span> #radius in pixels\n<\/span>win<\/span> =<\/span> np<\/span>.<\/span>ones<\/span>((<\/span>2<\/span>*<\/span> r<\/span> +<\/span>1<\/span>,<\/span> 2<\/span>*<\/span> r<\/span> +<\/span>1<\/span>))<\/span>\n# ----------   or, copy paste your window matrix -------------\n# win = np.array( [    [0, 1, 1, 1, 0]\n#                      [1, 1, 1, 1, 1],\n#                      [1, 1, 0, 1, 1],\n#                      [1, 1, 1, 1, 1],\n#                      [0, 1, 1, 1, 0]  ])\n<\/span>\n# window radius is needed for the function,\n# deduce from window size (can be different for height and width\u2026)\n<\/span>r_y<\/span>,<\/span> r_x<\/span>  =<\/span> win<\/span>.<\/span>shape<\/span>[<\/span>0<\/span>]<\/span>\/\/<\/span>2<\/span>,<\/span> win<\/span>.<\/span>shape<\/span>[<\/span>1<\/span>]<\/span>\/\/<\/span>2<\/span>\nwin<\/span>[<\/span>r_y<\/span>,<\/span> r_x<\/span>  ]<\/span>=<\/span>0<\/span>  # let's remove the central cell \n<\/span>\ndef<\/span> view<\/span> (<\/span>offset_y<\/span>,<\/span> offset_x<\/span>,<\/span> shape<\/span>,<\/span> step<\/span>=<\/span>1<\/span>):<\/span>\n    """\n    Function returning two matching numpy views for moving window routines.\n    - 'offset_y' and 'offset_x' refer to the shift in relation to the analysed (central) cell \n    - 'shape' are 2 dimensions of the data matrix (not of the window!)\n    - 'view_in' is the shifted view and 'view_out' is the position of central cells\n    (see on landscapeanalysis.org\/2018\/numpy-loops\/\n    """<\/span>\n    size_y<\/span>,<\/span> size_x<\/span> =<\/span> shape<\/span>\n    x<\/span>,<\/span> y<\/span> =<\/span> abs<\/span>(<\/span>offset_x<\/span>),<\/span> abs<\/span>(<\/span>offset_y<\/span>)<\/span>\n    \n    x_in<\/span> =<\/span> slice<\/span>(<\/span>x<\/span> ,<\/span> size_x<\/span>,<\/span> step<\/span>)<\/span> \n    x_out<\/span> =<\/span> slice<\/span>(<\/span>0<\/span>,<\/span> size_x<\/span> -<\/span> x<\/span>,<\/span> step<\/span>)<\/span>\n\n    y_in<\/span> =<\/span> slice<\/span>(<\/span>y<\/span>,<\/span> size_y<\/span>,<\/span> step<\/span>)<\/span>\n    y_out<\/span> =<\/span> slice<\/span>(<\/span>0<\/span>,<\/span> size_y<\/span> -<\/span> y<\/span>,<\/span> step<\/span>)<\/span>\n \n    # the swapping trick    \n<\/span>    if<\/span> offset_x<\/span> <<\/span> 0<\/span>:<\/span> x_in<\/span>,<\/span> x_out<\/span> =<\/span> x_out<\/span>,<\/span> x_in<\/span>                                 \n    if<\/span> offset_y<\/span> <<\/span> 0<\/span>:<\/span> y_in<\/span>,<\/span> y_out<\/span> =<\/span> y_out<\/span>,<\/span> y_in<\/span>\n \n    # return window view (in) and main view (out)\n<\/span>    return<\/span> np<\/span>.<\/span>s_<\/span>[<\/span>y_in<\/span>,<\/span> x_in<\/span>],<\/span> np<\/span>.<\/span>s_<\/span>[<\/span>y_out<\/span>,<\/span> x_out<\/span>]<\/span>\n\n# ----  main routine  -------\n<\/span>\ndem<\/span> =<\/span> gdal<\/span>.<\/span>Open<\/span>(<\/span>elevation_model<\/span>)<\/span>\nmx_z<\/span> =<\/span> dem<\/span>.<\/span>ReadAsArray<\/span>()<\/span>\n\n#matrices for temporary data\n<\/span>mx_temp<\/span> =<\/span> np<\/span>.<\/span>zeros<\/span>(<\/span>mx_z<\/span>.<\/span>shape<\/span>)<\/span>\nmx_count<\/span> =<\/span> np<\/span>.<\/span>zeros<\/span>(<\/span>mx_z<\/span>.<\/span>shape<\/span>)<\/span>\n\n# loop through window and accumulate values\n<\/span>for<\/span> (<\/span>y<\/span>,<\/span>x<\/span>),<\/span> weight<\/span> in<\/span> np<\/span>.<\/span>ndenumerate<\/span>(<\/span>win<\/span>):<\/span>\n    \n    if<\/span> weight<\/span> ==<\/span> 0<\/span> :<\/span> continue<\/span>  #skip zero values !\n<\/span>    # determine views to extract data \n<\/span>    view_in<\/span>,<\/span> view_out<\/span> =<\/span> view<\/span>(<\/span>y<\/span> -<\/span> r_y<\/span>,<\/span> x<\/span> -<\/span> r_x<\/span>,<\/span> mx_z<\/span>.<\/span>shape<\/span>)<\/span>\n    # using window weights (eg. for a Gaussian function)\n<\/span>    mx_temp<\/span>[<\/span>view_out<\/span>]<\/span> +=<\/span> mx_z<\/span>[<\/span>view_in<\/span>]<\/span>  *<\/span> weight<\/span>\n    \n   # track the number of neighbours \n<\/span>   # (this is used for weighted mean : \u03a3 weights*val \/ \u03a3 weights)\n<\/span>    mx_count<\/span>[<\/span>view_out<\/span>]<\/span> +=<\/span> weight<\/span>\n\n# this is TPI (spot height \u2013 average neighbourhood height)\n<\/span>out<\/span> =<\/span> mx_z<\/span> -<\/span> mx_temp<\/span> \/<\/span> mx_count<\/span>\n\n# writing output \n<\/span>driver<\/span> =<\/span> gdal<\/span>.<\/span>GetDriverByName<\/span>(<\/span>'GTiff'<\/span>)<\/span>\nds<\/span> =<\/span> driver<\/span>.<\/span>Create<\/span>(<\/span>output_model<\/span>,<\/span> mx_z<\/span>.<\/span>shape<\/span>[<\/span>1<\/span>],<\/span> mx_z<\/span>.<\/span>shape<\/span>[<\/span>0<\/span>],<\/span> 1<\/span>,<\/span> gdal<\/span>.<\/span>GDT_Float32<\/span>)<\/span>\nds<\/span>.<\/span>SetProjection<\/span>(<\/span>dem<\/span>.<\/span>GetProjection<\/span>())<\/span>\nds<\/span>.<\/span>SetGeoTransform<\/span>(<\/span>dem<\/span>.<\/span>GetGeoTransform<\/span>())<\/span>\nds<\/span>.<\/span>GetRasterBand<\/span>(<\/span>1<\/span>).<\/span>WriteArray<\/span>(<\/span>out<\/span>)<\/span>\nds<\/span> =<\/span> None<\/span>\n<\/code><\/pre>\n<\/div>\n<\/div>\n
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.<\/p>\n
It is not difficult to create a window on the fly:<\/p>\n
\n
\n
r<\/span> =<\/span> 5<\/span> #radius in pixels\n<\/span>win<\/span> =<\/span> np<\/span>.<\/span>ones<\/span>((<\/span>2<\/span>*<\/span>r<\/span> +<\/span>1<\/span>,<\/span> 2<\/span>*<\/span>r<\/span> +<\/span>1<\/span>))<\/span> #square window, all values are 1\n<\/span>win<\/span>[<\/span>r<\/span>,<\/span> r<\/span>]<\/span> =<\/span> 0<\/span>  #skipping the central cell\n<\/span><\/code><\/pre>\n<\/div>\n<\/div>\n
For Gaussian window weights: see numpy documentation<\/a>.<\/p>\n
However, I find more useful (and intuitive) to hard-code window definitions. This is very handy when comparing results on different datasets \u2013 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:<\/p>\n