Interpolación IDW (ponderada por el inverso de la distancia) en PyQGIS

El modelo de interpolación IDW estima los valores de una propiedad, que se asume depende de su localización, con base en la asignación de pesos a los datos del entorno en función inversa a la distancia que los separa del punto en cuestión. De esta forma se acepta que los puntos más próximos al centroide “z” intervienen de manera más relevante en la construcción del valor definitivo de la propiedad para ese punto. Por otra parte, si las distancias están a su vez afectadas por un exponente de ponderación, cuanto mayor es éste más contribuyen los puntos próximos.

La interpolación ponderada con base al inverso de la distancia se basa en esta formula:

pantallazo2

donde zj es el punto problema, zi es un punto del entorno, β es el exponente de ponderación y dij es la distancia entre los puntos.

Aprovechando parte del código de la determinación de la matriz de distancia entre puntos, se propuso el siguiente para interpolar los valores de “elevación” para tres puntos circunscritos (en verde) dentro de otros 6 (en azul); todos escogidos de manera arbitraria y tal como se observan en la imagen más adelante.

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)

Los puntos en verde de la imagen siguiente son los que se encuentran en el código cómo QgsPoint y los azules tienen en su tabla atributiva los valores de elevación que han sido incorporados allí utilizando un plugin para muestreo ráster. Estos últimos son los que se emplearán como base para la interpolación.

idw1

Cuando se ejecuta el código en la Python Console se obtiene:

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

donde la lista son los valores de elevación en el vectorial y los tres últimos valores las elevaciones interpoladas para los puntos en verde. Estos valores, tal como se esperaba, son similares a sus vecinos más cercanos.

El código anterior se programó en un plugin de QGIS:

idw3

y para prever un mayor control en la verificación de los valores de salida se redujo la zona y se seleccionaron más puntos (50) como base para la interpolación en la generación del ráster correspondiente. El resultado fue el siguiente:

idw2

y es idéntico, aunque ligeramente desplazado, al que se obtiene con la herramienta de interpolación de QGIS. La ventaja que se obtiene programando el plugin es que estará alineado y tendrá la misma resolución que el ráster base.

En nuestro curso de PyQGIS puedes aprender a programar éste y otros plugins mediante el empleo de ‘Plugin Builder’.

PRÓXIMA SESIÓN EL 20 DE OCTUBRE

Esta entrada fue publicada en PyQGIS, QGIS, SIG, Software Libre. Guarda el enlace permanente.

Una respuesta a Interpolación IDW (ponderada por el inverso de la distancia) en PyQGIS

  1. Pingback: Newsletter Octubre 2015 | El Blog de José Guerrero

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s