2

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?

4
  • 1
    I agree that scale (or radius) is important. Just tried to read QGIS documentation and gdaldem tpi. I felt the same as you probably did, that gdal was simply looking at its neighbor cells. You might have to change the size of DEM pixels... probably not a preferred option. Instead, SAGA will give you more control. Processing | SAGA | Terrain Analysis | Morphometry | Topographic position index (tpi).
    – Kazuhito
    Dec 1, 2016 at 16:18
  • @Kazuhito, yes that's what I thought it was considering, unfortunately... I didn't try with SAGA though. And probably GRASS will have more hypotheses as well! Thanks for your suggestion!
    – mto23
    Dec 1, 2016 at 22:20
  • 1
    I would just use raster algebra to obtain the desired result. The TIP is just the focal cell subtracted from the focal mean window. You can calculate a focal mean, using any window desired, then subtract the DEM from the resulting focal raster and you have the TPI. Dec 2, 2016 at 17:03
  • 1
    Raster algebra solutions are detailed at gis.stackexchange.com/questions/6056.
    – whuber
    Dec 3, 2016 at 18:38

2 Answers 2

7

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.

enter image description here

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.

2
  • Very helpful, thanks Kazuhito. Still there is a problem and I'm not sure where: Following your instructions I get an error when running the tool: "invalid syntax (, line 17) See log for more details" (well actually this is the log...). I would appreciate your further help.
    – esther
    Jul 24, 2017 at 14:33
  • @esther Apologies... I cannot test it myself right now. It can be from various causes; most likely my faulty script, or recent change in raster package, or configuration in the QGIS processing. Will you post it as new question? If the focus is on the error message, I guess it will not be marked as dupulicate.
    – Kazuhito
    Jul 25, 2017 at 8:08
0

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

the final model with 3 inputs, annulus mask, r.neighbors, raster calculator and TPI as outputs

"""
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()


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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