Standard scale value for TPI?

Question

The estimation of the Topographic Position Index (TPI) is scale dependent, but when estimating this index in QGIS (Raster > Analyses > DEM (Terrain models)), I don't see where I can change this scale, nor its standard value.

Do you know how do I obtain/change that(if possible)? Is it related with the pixel value used in the input DEM raster?

Answer 1

I do not know much about GRASS, but given there are variety of tools such as r.neighbors (focal statistics), I'm sure you already have found the solution with GRASS.

BTW, this is just to follow up on the other lead - comments by Jeffrey and whuber, both recommended raster algebra.

I tried to put it in R script, which can be run by copy/paste with Script Editor - via Processing Toolbox | R scripts | Tools | Create new R script.

TPI.rsx -- please copy this to Script Editor

##Raster_Tools= group
##Input_Raster_Layer= raster
##Band_number= number 1
##Outer_Window_Size= number 3
##Inner_Window_Size= optional number 1
##Compute_edges= boolean TRUE
##Topographic_Position_Index= output raster

library(raster)
ow <- Outer_Window_Size
iw  <- Inner_Window_Size
edge <- Compute_edges
Layer <- Input_Raster_Layer[[Band_number]]
Outer_window <- focal(Layer, w=matrix(1, nrow=ow, ncol=ow), 
                      fun=mean, na.rm=FALSE, NAonly=TRUE, pad=edge)
if(iw>=3){Inner_window = focal(Layer, w=matrix(1, nrow=iw, ncol=iw), 
                      fun=mean, na.rm=FALSE, NAonly=TRUE, pad=edge)}
if(iw==1){Inner_window = Layer}
Topographic_Position_Index <- Inner_window - Outer_window

If you save this file into ~/.qgis2/processing/rscripts/ as TPI.rsxthere would appear a new geoalgorithm TPI under Raster Tools group.

Note

1. Outer window size: default= 3 (cells)...odd number 3, 5, 7, ...
2. Inner window size: default= 1 (cell)...odd number 1, 3, 5, ...

If we run it with Outer=3 and Inner= 1, it produces a raster layer which is equivalent to the gdaldem tpi. Please modify Outer window size as necessary.

Usually Inner window size would always be 1, unless you need to work with very fine grid (e.g. LiDAR) which can make your TPI analysis difficult.

Answer 2

Late to the party, but you can also easily implement the TPI yourself and set all values yourself (both inner and outer radii). here is a solution (more thoroughly explained here):

1. Add the script on the bottom as a tool to your processing toolbox
2. Create a model
3. drag in 2 "number" inputs and a "Raster layer" input
4. drag in your annular mask model
   1. use the number inputs as inner, resp. outer radii
5. drag in an r.neighbors
   1. use the DEM as Input raster layer
   2. use outer radius as Neighborhood size
   3. set neighborhood operation to average
6. drag in a GDAL raster calculator to calculate DEM-Neighbors and the output is the TPI

"""
Model exported as python.
Name : annulus mask for r.neighbours
Group : 
With QGIS : 32200
"""

from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterNumber
from qgis.core import QgsProcessingParameterFileDestination
import processing

from numpy import sqrt, fromfunction, logical_and, savetxt

class AnnulusMaskForRneighbours(QgsProcessingAlgorithm):

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterNumber('innerradius', 'inner radius', type=QgsProcessingParameterNumber.Integer, minValue=1, defaultValue=1))
        self.addParameter(QgsProcessingParameterNumber('outerradius', 'outer radius', type=QgsProcessingParameterNumber.Integer, minValue=2, defaultValue=3))
        self.addParameter(QgsProcessingParameterFileDestination('outfile', 'annular mask'))

    def processAlgorithm(self, parameters, context, model_feedback):
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(0, model_feedback)
        
        r_in = self.parameterAsInt(parameters, 'innerradius', context)
        r_out = self.parameterAsInt(parameters, 'outerradius', context)
        
        outloc = self.parameterAsString(parameters, 'outfile', context) #+ ".txt"
        
        d = sqrt(fromfunction(lambda x,y: (x-r_out)**2+(y-r_out)**2,(2*r_out+1,2*r_out+1)))
        m = logical_and(d>=r_in, d<r_out)
        feedback.pushInfo(str(m))
        savetxt(outloc,m,fmt="%d")
        
        results = {}
        outputs = {}

        return results

    def name(self):
        return 'annulus mask for r.neighbours'

    def displayName(self):
        return 'annulus mask for r.neighbours'

    def group(self):
        return ''

    def groupId(self):
        return ''

    def createInstance(self):
        return AnnulusMaskForRneighbours()