How to compute an interpolation raster from the python console in QGIS?

Question

I am trying to add an interpolation layer based on an existing vector layer by using the Python Console (initially in the menu Raster/Interpolation), but can't find anything related in the documentation.

I saw the QgsIDWInterpolator in the API, but can't get it working. It seems it expects a list of LayerData but it doesn't work neither. Here my current code:

ld = qgis.analysis.QgsIDWInterpolator.LayerData
ld.vectorLayer = vlayer
ld.zCoordInterpolation = 0
ld.interpolationAttribute = 5
ld.mInputType = 0
ldlist = [ld]
itp = qgis.analysis.QgsIDWInterpolator(ldlist)

I get the following error for the last command:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
TypeError: arguments did not match any overloaded call:
  QgsIDWInterpolator(list-of-QgsInterpolator.LayerData): argument 1 has unexpected type     'list'
  QgsIDWInterpolator(QgsIDWInterpolator): argument 1 has unexpected type       'list'

Answer 1

The documentation on pyqgis is not very self-explanatory, but i figured out how to properly call the associated interpolation classes (QgsInterpolator, QgsTINInterpolator, QgsIDWInterpolator, QgsGridFileWriter) from python. I am going to describe every step of the script in great detail:

Step 1:

Import the core and analysis module and get the desired vector layer for interpolation by selecting it with a mouseclick in the layer tab.

import qgis.core
import qgis.analysis
layer = qgis.utils.iface.activeLayer()

Step 2:

Prepare the interpolation classes with the necessary Parameters. The exact parameters for initialization of the LayerData struct can be found in the QGIS API docs (searchterm: QgsInterpolator).

layer_data = QgsInterpolator.LayerData()
layer_data.vectorLayer = layer
layer_data.zCoordInterpolation=False
layer_data.interpolationAttribute =0
layer_data.mInputType = 1

Please notice that I don't use the z Coordinate, I get the first available field (index = 0) as interpolation attribute, and use POINTS as input type.

Step 3:

Choose your interpolation engine. Here you can choose between the TIN-Interpolation method (QgsTINInterpolator) and IDW-Interpolation (QgsIDWInterpolator). I took the QgsTINInterpolator in my code.

tin_interpolator = QgsTINInterpolator([layer_data])

Keep in mind that you have to pass a python list of layer_data to the interpolation engine! This also allows you to add multiple layer_data scenarios.

Step 4:

Setup the parameters that are needed for the export of the interpolation-output (see documentation of QgsGridFileWriter). Those include similar information as the interpolation gui (filepath, extent, resolution, number of colums and rows).

export_path ="C:/SomeFolder/output.asc"
rect = layer.extent()
res = 10
ncol = int( ( rect.xMaximum() - rect.xMinimum() ) / res )
nrows = int( (rect.yMaximum() - rect.yMinimum() ) / res)
output = QgsGridFileWriter(tin_interpolator,export_path,rect,ncol, nrows,res,res)
output.writeFile(True)
iface.addRasterLayer(export_path, "interpolation_output") 

Be aware of the file extension of your output-raster as QgsGridFileWriter only writes ASCII-grids (.asc). The data gets written to disk by calling the writeFile() method. After export you can add the grid-file as raster to the canvas.

Full script for reference:

import qgis.analysis
import qgis.core
layer = qgis.utils.iface.activeLayer() 
layer_data = QgsInterpolator.LayerData()
layer_data.vectorLayer = layer
layer_data.zCoordInterpolation=False
layer_data.interpolationAttribute =0
layer_data.mInputType = 1
tin_interpolator = QgsTINInterpolator([layer_data])
export_path = "E:/GIS_Workbench/script_output/test.asc"
rect = layer.extent()
res = 10
ncol = int( ( rect.xMaximum() - rect.xMinimum() ) / res )
nrows = int( (rect.yMaximum() - rect.yMinimum() ) / res)
output = QgsGridFileWriter(tin_interpolator,export_path,rect,ncol,nrows,res,res)
output.writeFile(True)

Keep in mind that the QGIS-API is currently rewritten to version 3.0 and the used interpolation-classes are moved from qgis.analysis to qgis.core! This will have a huge impact on the functionality of this script so that it must be rewritten for version 3.0!

Answer 2

To reiterate what Kyungdog Kim said, there's a typo in root676's answer. The line layer_data.InterpolationAttribute =0 should be replaced with layer_data.interpolationAttribute =0. The first "i" in ".interpolationAttribute" should be lowercase.

In root676's example code, I believe it just happened to be that the required interpolation attribute (the 0-index) coincided with the method's default behavior (and I'm not sure why an exception wasn't raised).

If you were to try putting that block of code into a loop and iterating over the first 5 columns of the attribute table you would be interpolating the first column 5 times.

for i in np.arange(0, 5):
  <other relevant code here>
  layer_data.InterpolationAttribute =i

However, if you correct the typo and again try iterating over the first 5 columns of the attribute table, you now correctly traverse the columns.

for i in np.arange(0, 5):
  <other relevant code here>
  layer_data.interpolationAttribute =i

See documentation for QgsInterpolator::LayerData::interpolationAttribute method here: https://qgis.org/api/structQgsInterpolator_1_1LayerData.html#ad7fca7c06b5787b81fa6e10441045853

Answer 3

Nevermind, found it in the documentation : http://docs.qgis.org/testing/en/docs/user_manual/processing/console.html

Answer 4

InterpolationAttribute -> interpolationAttribute

It should be lowercase than uppercase. Otherwise, the attribute index is automatically set to the default value of 0