3

Yesterday, I try after a long time to make an IDW interpolation on a point layer in Qgis 2.18 with SAGA tools. I use same parameters as before on Qgis 2.14, but the result is quite different. Raster isn't smooth and all points are very visible on it.

IDW interpolation

Has SAGA IDW interpolation has been changed between these 2 versions?

Improve this question
4
  • How about changing the Resampling method (in Advanced parameters) to Nearest Neighbour?
    – Kazuhito
    Oct 25, 2017 at 9:43
  • Which parameter is it for resampling ? because parameter names aren't explicit in IDW interpolation window. I have Search Range with Maximum number of nearest points and Number of points with Maximum number of nearest points too. I already use these 2 settings...
    – Tim C.
    Oct 25, 2017 at 11:39
  • which version of 2.18 are you on? If you're on 2.18.10 or later QGIS may be using a different version of SAGA. Do you know what the Inverse Distance Power and Inverse Distance weighting values before and after?
    – Steven Kay
    Oct 25, 2017 at 18:16
  • I'm on QGIS 2.18.9 with SAGA 2.2.3 and interpolation "do not work". But I don't know Qgis and Saga versions when interpolation was working well. I only know values of parameters for old version. For Inverse Distance Power : 2 and I haven't Inverse Distance weighting but Distance Weighting on inverse distance to a power.
    – Tim C.
    Oct 26, 2017 at 7:50

1 Answer 1

Reset to default
4

You can see distance-based interpolation formula in following link:

http://www.gitta.info/ContiSpatVar/en/html/Interpolatio_learningObject2.xhtml

where p is distance weighting exponent and it probably was the parameter that is different between these 2 SAGA versions.

So, it's relatively easy to code this algorithm in PyQGIS. As follow, it has inside three arbitrary points to be interpolated and it was used a five points layer as base of interpolation.

from math import sqrt

def idw(pt, points, z, beta):

    n = len(points)

    sum_inv_dist_z = 0
    sum_inv_dist = 0

    for i in range(n):

        inv_dist = 1/(sqrt(pt.sqrDist(points[i])))**beta
        inv_dist_z = z[i]*inv_dist

        sum_inv_dist_z += inv_dist_z
        sum_inv_dist += inv_dist

    z_j = sum_inv_dist_z/sum_inv_dist

    return z_j

#Code starts here

layer = iface.activeLayer()

features = layer.getFeatures()

#get features for geometry
points = [feature.geometry().asPoint() for feature in features]

#get features for attributes
features = layer.getFeatures()

attr = [feature.attributes() for feature in features]

#selecting elevation field
z = [item[4] for item in attr]

print z

#set beta value
beta = 2

#defining arbitrary points
pt = [QgsPoint(398574.065172, 4452409.334),
      QgsPoint(412049.303287, 4447317.71399),
      QgsPoint(434741.145557, 4435557.70814)]

#calculating z value for each point
for item in pt:
    print idw(item, points, z, beta)

Used layers look as at following image where green points will be interpolated and blue points have z values.

enter image description here

After running code at Python Console of QGIS I got:

[1694.0, 1700.0, 1368.0, 1439.0, 1368.0, 1943.0]
1709.79209934
1622.47616106
1436.15697073

where it can be observed that obtained values are similar to neighbors.

However, it's more useful when above algorithm was coded in my own IDW plugin to get a complete raster (perfectly aligned with base raster). In following composition I used 25 points and p different values (0.1, 0.5 and 2). It can be observed that p parameter is a smoothed factor.

enter image description here

By the way, my results are identical to obtained with IDW QGIS tool; except my interpolated rasters were perfectly aligned with base raster.

Improve this answer

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.