3

I am trying to iteratively interpolate a shapefile using the QGIS IDW interpolator. What I need to do is take a point shapefile with time series data in the columns associated with each point, interpolate that column of data to raster, then move to the next column, creating a new raster for each column.

The process I've generally been following is outlined here: How to compute an interpolation raster from the python console in QGIS?

I've modified a few things and have been working on a test case (since the input data hasn't been provided yet - format should be similar except ~300 iterations), but I can't get the interpolator to move past the first, 0-indexed, column.

Code runs without error, outputs the 4 test raster interpolation .asc files, and adds them to the QGIS layers panel, but they are all the same. You can see in the layers panel snippet below that the raster cell bounds are all the same and visually if I toggle the layers on/off, they are all the same.

Code below (only change is that I modified actual folder paths to dummy variables):

import qgis.core
import qgis.analysis
import os


pathToFile = r"input_folder\interp_test.shp"

interp_layer = QgsVectorLayer(pathToFile, 'interp_test','ogr')

#not used yet
field_names = [field.name() for field in interp_layer.pendingFields()]

#test first 4 columns
for i in np.arange(0, 4):
    layer_data = qgis.analysis.QgsInterpolator.LayerData()
    layer_data.vectorLayer = interp_layer
    layer_data.zCoordInterpolation=False

    #this line should change the interpolation attribute/column, but it doesn't
    layer_data.InterpolationAttribute = i
    layer_data.mInputType = 1 #points

    #pick interpolator engine
    idw_interpolator = qgis.analysis.QgsIDWInterpolator([layer_data])

    export_path =r"output_folder\rasterout_{}.asc".format(i)
    rect = interp_layer.extent()
    res = 500
    ncol = int( ( rect.xMaximum() - rect.xMinimum() ) / res )
    nrows = int( (rect.yMaximum() - rect.yMinimum() ) / res)

    output = qgis.analysis.QgsGridFileWriter(idw_interpolator,export_path,rect,ncol, nrows,res,res)
    output.writeFile(True)  
    iface.addRasterLayer(export_path, "interpolation_output_{}.asc".format(i))

Input test attribute table: Input attribute table w/data to be interpolated

Output to QGIS console: Output interpolated raster, duplicated 4x

Running 64-bit QGIS 2.18.3 on windows, using the built-in python console.

1
  • Still interested in figuring out why the above code doesn't work, but for now working around this issue by setting layer_data.InterpolationAttribute = 0 and deleting the first column of data in the attribute table. Seems a bit slower than the previous version that was interpolating the first column repeatedly, but it's working.
    – nix006
    Commented Aug 7, 2017 at 18:23

1 Answer 1

4

I think you have a slight typo for the QgsInterpolator::LayerData::interpolationAttribute method.

You used:

layer_data.InterpolationAttribute = i

where "Interpolation" should start with a lower case:

layer_data.interpolationAttribute = i

Your code worked for me on QGIS 2.18.2 Win7 64-bit.

2
  • 1
    That works for me too. Thank you very much for the help! So it seems like the default behavior for the method is to use the 0-indexed column as the interpolation attribute, would you happen to know why an exception wasn't raised for the typo? Still pretty new to QGIS and python coding in general.
    – nix006
    Commented Aug 9, 2017 at 15:27
  • @nix006 - Most welcome, glad it helped! I'm not sure why it doesn't raise an exception, I have seen a couple of situations like this. Guess it's just the way the API works ;)
    – Joseph
    Commented Aug 10, 2017 at 11:32

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.