5.1. Follow Along: Creating a script for the mask file

Now, to get the r.neighbors algorithm to work correctly, we need to create a mask file script.

Warning

This is a exercise. Only do this if you have extra time left. Otherwise, go directly to the solution. Doing this exercise will also help you with Follow Along: Making a script for batch rasterizing.

5.1.1. Creating a model

1. We will create a model that we will convert to a script. Click ► Create New Model…

2. Drag in the following inputs:

   • Number:

     • Description: Inner radius

     • Number type: Integer

     • Minimum value: 0

     • Default value: 1

   • Number: Outer radius, similar to Inner radius

3. Name the model Annulus mask for r.neighbors. Optionally, also give it a group name.

4. Click the Export as Script Algorithm icon.

   The following script will appear. Places where we will insert some of our own code are highlighted.

5.1.2. Add additional parameters

In the model, we have only added two inputs. However, our algorithm should also have an output. At line 11, insert the following:

from qgis.core import QgsProcessingParameterFileDestination

and within the initAlgorithm function (line 18) insert:

                "innerradius",

5.1.3. Convert parameters to workable format

There is one problem with these processing parameters, however: they are not actually values that we can work with. However, we want to be able to use them as numbers or strings (in the case of file names). For this we will use the parameterAsInt() and parameterAsFileDestination() At line 23, insert the following.:

        self.addParameter(
            QgsProcessingParameterNumber(
                "outerradius",
                "outer radius",

5.1.4. Perform calculations

What we want is a function that creates a file like below for an inner, resp. outer radius of: \(r_i=1,r_o=3\)

0 0 1 0 0
0 1 1 1 0
1 1 0 1 1
0 1 1 1 0
0 0 1 0 0

This file is a mask file with weights e.g. numbers between 0 and 1, that tell GRASS how much this cell matters for the calculation of the tpi.

Note

I did not come up with these calculations myself, but found them on stackexchange. Sadly, I forgot where.

Note that the 0 in between all the 1s is the center is the point that corresponds to the center. It is actually at coordinates \((x_0,y_0)=(3,3)\) (start counting at 0). This is what r.neighbors expects. It follows that \(x_0=y_0=r_o\). Let \(d\) be the distance to this point. Then, we want all points to be 1 for which:

\[d\geq r_i \land d \leq r_0\]

holds and 0 otherwise. The (eucludian) distance can be calculated by:

\[\begin{split}d := \sqrt{(x-x_0)^2+(y-y_0)^2}\\ = \sqrt{(x-r_o)^2+(y-r_o)^2}\end{split}\]

where \(x,y\) are the coordinates of the Currently processed point. To put this in code, we first need to import the corresponding functions:

from numpy import sqrt, fromfunction, logical_and, savetxt

and then make the calculations. Here a ** b means \(a^b\). Also note that the size of our array is \(2r_o+1\)

                minValue=2,
                defaultValue=3,

Then, we save our file to outloc in decimal ("%d") format:

        )

Your final script should look like this:

Solution

If you didn’t follow the above Follow Along:, you can use the below script.

1. In the Processig Toolkbox, click the ► Create New Script…

2. copy-paste the following code into the text editor that popped up:

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

3. Save the script. It should now show up in the toolbox:

   

5.1.5. Testing the annulus mask

Now you have made the annulus mask file, either by following the instructions or skipping to the solution, now it is time to test whether the annulus mask we just made actually works.

1. Search for annulus mask for r.neighbors in the processing pane and run it.

   Use default settings, but set annular mask to a file name with a .txt extension

   

2. If you have any errors, read the error, see where it comes from and resolve them.

3. open the file and verify if it is made correctly.