5

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'
3

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!

1

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

-1

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

We're looking for long answers that provide some explanation and context. Don't just give a one-line answer; explain why your answer is right, ideally with citations. Answers that don't include explanations may be removed.

  • 2
    Could you please share your answer in full detail! The link to the console reference tells nothing about your solution of the python problem regarding the interpolation. Thanks! – root676 Mar 21 '17 at 21:28
  • I figured out the solution myself and posted a detailed answer - ready to use for other users who encounter the same problem as in the question. – root676 Mar 23 '17 at 21:40
-2

InterpolationAttribute -> interpolationAttribute

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

We're looking for long answers that provide some explanation and context. Don't just give a one-line answer; explain why your answer is right, ideally with citations. Answers that don't include explanations may be removed.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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