Cómo usar el lenguage R con Python y QGIS-GRASS

Recientemente, un usuario en un foro de GIS planteaba la necesidad de corregir un ráster DEM, en una zona delimitada por un archivo vectorial, con base en que éste no reflejaba la elevación producto de una densa población de árboles de unos 30 metros de altura como valor promedio. Hacer una corrección del DEM con base a la resta de 30 m dentro de la zona del perímetro no es difícil pero el efecto que se produce es como el de un “zurcido invisible”.

Por tanto, para tener un efecto más “natural” lo que se plantea es generar una serie de puntos aleatorios en el área del vectorial, incorporarle la columna de altura de los arboles en la tabla atributiva con la variable aleatoria (altura de los árboles) y luego interpolar para “tejer” el DEM.

Para determinar la variable aleatoria es necesario partir de la suposición que la altura de los árboles sigue un patrón conocido de distribución. La manera más sencilla de abordar ésto es asumir (puede que no sea cierto) que es normal. Una distribución normal de media μ y desviación típica σ se designa por N(μ, σ) y su gráfica es la “campana de Gauss”. La distribución normal estándar es aquella que tiene por media el valor cero, μ =0, y por desviación típica la unidad, σ =1. El área de esta última es uno y está relacionada con probabilidades. Para hacer referencia a probabilidades se tiene que transformar la variable aleatoria (X) que sigue una distribución N(μ, σ) en otra variable Z que siga una distribución N(0, 1). La expresión es X = Z.σ + μ. Por tanto, lo que resta por conocer en nuestro caso es el valor de σ (se determina experimentalmente y lo asumimos como 5 m) porque los valores de Z los producimos con un generador de números pseudoaleatorios con media cero, μ =0, y con desviación típica uno, σ =1.

Una manera expedita de producir números pseudoaleatorios con media cero, μ =0, y con desviación típica uno, σ =1, es mediante el comando rnorm del lenguaje R. Para usar R con Python y el plugin de Spqr de QGIS es necesario instalar en Debian los paquetes python-rpy y python-rpy2.

Estas simples intrucciones:

from rpy2.robjects import r
x = r('x <- rnorm(100)')  
print x

a pesar de que la visualización es en formato de texto de R, internamente se interpreta como un arreglo de floats y se puede usar directamente, en un bucle, para llenar la lista de la variable aleatoria mediante la expresión X = Z.σ + μ. Luego, se copian a la tabla atributiva del vectorial de puntos mediante db.execute de GRASS. El código completo es el siguiente:

#!/usr/bin/env python
#-*- coding: utf-8

from os import system
from rpy2.robjects import r
import grass.script as grass
from math import sqrt

system('clear')

#Valores de números pseudoaleatorios, z, a generar
n = 1000

z = r('x <- rnorm(' + str(n) + ')')

#Valores asumidos para la distribución
mean = 30
sdev = 5

#Declaración de lista vacia para la variable aleatoria
X = []

sum = 0
sum2 = 0

for i in range(n):
	temp = z[i] * sdev + mean
	X.append(temp)
	sum += temp
	sum2 += temp*temp 

#Determina promedio y desviación estándar del arreglo generado
#para verificar si tiene media 30 y desviación estandar 5
promedio = sum/n
desvestand = sqrt((sum -sum2/n)/(n-1))

#Parámetros estadísticos del arreglo de 1000 valores generados
print "Parámetros estadísticos del arreglo de 1000 valores generados"
print "media = %.2f sdev = %.2f" % (promedio, desvestand)

#Se obtiene información sobre las columnas del vectorial
info= grass.read_command("v.info", 
                    flags="c", 
                    map="puntos_poligono_vector",
                    quiet="True")
  
info= info.split("\n")

#Crear columna variable si no existe
val = 0
  
for strings in info:
    if 'variable' in strings:
        val = 1
        print "La columna variable existe"
        break
  
if val == 0:
    grass.run_command("v.db.addcol", 
                        map="puntos_poligono_vector", 
                        columns="variable double precision")
  
    print "Se creó la columna variable"

print "Se está escribiendo el arreglo a la tabla..."
 
#Introduce los valores del arreglo array en la columna variable del vectorial
for i in range(n):
    input="UPDATE puntos_poligono_vector SET variable=" + str(X[i]) + " WHERE cat = " + str(i+1)
    grass.write_command("db.execute", stdin = input)

print "Listo!" 

Su ejecución, en cónsola, permite verificar que la media del arreglo de 1000 valores generados tiene como parámetros valores cercanos a los supuestos.

Parámetros estadísticos del arreglo de 1000 valores generados
media = 29.81 sdev = 5.38
La columna variable existe
Se está escribiendo el arreglo a la tabla...
Listo!

La imagen siguiente permite verificar que los valores del arreglo fueron copiados a la columna variable de la tabla atributiva y el histograma, producido mediante el plugin Spqr, refleja una distribución normal.

R

Esta entrada fue publicada en Código Python, GRASS, Lenguaje R, QGIS, SIG, Software Libre. Guarda el enlace permanente.

8 respuestas a Cómo usar el lenguage R con Python y QGIS-GRASS

  1. juan dijo:

    Hola Jose, como se podria hacer estos pasos en Windows?
    Muchas Gracias
    Salu2

  2. Mira, yo aprendí por mi cuenta, sin ningún tutorial, instalando distribuciones de Linux en discos viejos y de poca capacidad. Hoy en día, con la capacidad que tienen las pendrives y las tarjetas de memoria puedes probar con ellas el particionado. Aquí tienes un ejemplo:

    https://joseguerreroa.wordpress.com/2010/06/10/instalar-linux-en-el-disco-duro-uso-previo-de-gparted/

    Eso si. ASEGÚRATE DE CUAL ES EL NOMBRE DE TU DISPOSITIVO REMOVIBLE PARA QUE NO VAYAS A CREAR UNA PARTICIÓN SOBRE TU WINDOWS. No hay que tener miedo. Sólo hay que ser previsivo y hacer click cuando se está seguro de que es lo que se quiere hacer. De todas maneras, el sistema particionador siempre avisa que es lo que va a hacer antes de “ese click”. Si la “embarras” es por tu culpa (al no entender las instrucciones) y si tienes problemas con el inglés asegúrate que las instrucciones de particionado están en castellano.

    Por otra parte, si tienes Windows Vista, 7 o superior la mejor manera de particionar es desde el propio Windows. Estas son la instrucciones:

    http://www.taringa.net/posts/info/5524138/Particionar-disco-duro-en-Windows-Vista-Windows-7.html

    Después, con el Gparted de un CD Live de Linux, se formatea ese espacio creado por Windows con un sistema de archivos compatible a Linux (puede ser ext3 o reiserfs, por ejemplo) para que en la instalación sea ese (y no otro) el que “busque” por defecto.

  3. david arturomendez matute dijo:

    mugracias me agrado el tema

  4. Pingback: Cómo usar el lenguage R con Python y QGI...

  5. Buenas, Me podrias ayudar a instalar R en python 3.5 y los comandos para importalo a python …. el sistema operativo Windows … gracias y espero puedas ayudarme.

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