[GIS] SAGA IDW interpolation algorithm has been changed

interpolationinverse-distance-weightedqgissaga

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?

Best Answer

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.

Related Question