Álgebra de mapas con operadores relacionales y estructuras condicionales en GRASS con QGIS. Mapa topográfico a partir del DEM sin interpolación.

En el foro de Gabriel Ortiz un usuario señalaba que disponía de dos ráster, donde al primero lo llamaremos simplemente mareas y el otro es un DEM, y requería obtener otro ráster en el cual se volcaran los valores de mareas si se cumplía píxel a píxel la condición mareas < DEM y ceros o NULL en caso contrario. Aunque la complejidad del procedimiento no es muy grande, es interesante considerarlo porque permite poner en práctica el álgebra de mapas con operadores relacionales y estructuras condicionales; además de vislumbrar otro posible procedimiento para obtener curvas de nivel sin recurrir a métodos de interpolación. Es obvio que no se dispone del ráster mareas pero eso no impide que se pruebe el procedimiento con un ráster de valor constante en nuestro ya conocido dataset tiznados. En cónsola de GRASS se comenzará por crear uno de valor 500 (que representaría la cota correspondiente):

r.mapcalc "area_corte500 = 500"

Ahora se va a analizar el requerimiento. Si se establece la relación de mareas < DEM entonces el mapa de salida corresponde a una imagen binaria con píxeles iguales a uno donde se cumple la relación y ceros en caso contrario. Esta salida puede ser dirigida hacia un if donde si se cumple la condición (valores diferentes de cero) tome los valores de mareas y si son iguales a cero se establezcan a NULL (es lo más conveniente a la hora de vectorizar un ráster). En nuestro caso se va a sustituir la relación por “menor o igual que” (<=) y a cambiar los nombres de los ráster (mareas = area_corte500 = … = area_corte1000). El comando en cónsola de GRASS para el "corte", con el módulo r.mapcalc, sería:

r.mapcalc "corte_500 = (if((area_corte500<=demN09W068_UTM19N_canoa),area_corte500,null()))"

En la imagen siguiente se observa el resultado de corte_500, sobrepuesta sobre el DEM, que permite verificar que los píxeles que no cumplieron la relación y condición fueron convertidos a NULL.

Si se vectoriza el ráster con r.to.vect.area, se le elimina el relleno y se repite todo el proceso, cada 100 m, para las cotas entre 600 y 1000, se obtiene el siguiente resultado:

Es el mapa topográfico con curvas de nivel cada 100 m entre las cotas 500 y 1000 m. Aunque esto pudiese ser considerado más engorroso que usar directamente la herramienta Curvas de Nivel de QGIS podría ser un proceso de regeneración más exacto que aquellos que emplean interpolación (sobre todo con un DEM con una resolución menor a la del aquí usado que fue aproximadamente de 30×30). Por otra parte, hay que recordar que todo esto puede automatizarse mediante un script de bash para reducir el tiempo de regeneración de las curvas de nivel a través de líneas de comando.

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

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