Determinación del área 3D de una montaña o elevación mediante un script de python (GRASS-QGIS)

La estimación del área 3D de una superficie en GRASS se realiza con el módulo r.surf.area. Sin embargo, además de otros posibles errores, en el manual se admite la subestimación de la mitad del área del perímetro en la determinación. Para evaluar la magnitud del error se puede crear en el SIG el DEM de un sólido geométrico hipotético, por ejemplo un cono, para el cual se pueda determinar, por consideraciones puramente geométricas, la superficie 3D. En este caso, el procedimiento involucraría la determinación de la superficie 3D mediante r.surf.area, a la cual luego se le sumaría la mitad del área correspondiente al perímetro ráster de la base del cono.

Por ejemplo, si se asume la presencia de una montaña perfectamente cónica con la cúspide en la cota 800 m y en la cota 600 m la distancia al centro del círculo base es de 500 m (radio = r), la altura del cono (h) será de 200 m y su superficie 3D (S3D) vendrá dada por π.r.g; donde g es la generatriz del cono (hipotenusa del triángulo rectángulo que conforman el radio y la altura). Por tanto S3D = 3,1416*500*raiz(5002 + 2002) = 845899,709823203 m2.

Para modelar esto en el SIG primero se obtiene un ráster de círculos concéntricos (circulo) en un punto seleccionado arbitrariamente sobre el área de trabajo del dataset mediante el módulo r.circle; donde el radio del círculo mínimo fue 0 m y el del círculo máximo 500 m. Sin embargo, para determinar el área 3D se requiere un ráster DEM (circulo_dem)y para ello hay que transformarlo mediante álgebra de mapas. Si se asume que la cúspide de la montaña está en la cota 800 m y la base en la cota 600 m, la altura de la montaña cónica es de 200 m y la disminución de la altura en cada anillo concéntrico es proporcional a la tangente del ángulo que forman la generatriz y el radio, es decir, 200/500 = 0.4. Expresado en álgebra de mapas con r.mapcalc sería:

r.mapcalc "circulo_dem=800-0.4*circulo"

Para este insumo se generó el siguiente código:

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

from os import system 

import grass.script as grass

system("clear")

print "Espere..."

#Determina el área 3D
texto=grass.read_command("r.surf.area",
					input="circulo_dem",
					quiet="true")

#Filtra la salida para obtener sólo el área 3D
texto=texto.split("\n")

texto=texto[6].split(":")

area_sub=float(texto[1])

#Se obtiene un ráster unitario a partir del área del original
#Sirve para obtener un vectorial monocategoría
grass.mapcalc("temporal=if($circulo_dem, 1)", 
				circulo_dem="circulo_dem",
				quiet="true")

#Se obtiene un vectorial monocategoría
grass.run_command("r.to.vect",
					input="temporal",
					output="circulo_dem_area",
					feature="area",
					quiet="true")

grass.run_command("g.remove",
					rast="temporal",
					quiet="true")

#Un vectorial monocategoría tipo área producirá una línea monocategoría
#pero de categoría cero
grass.run_command("v.type",
					input="circulo_dem_area",
					output="temporal",
					type="boundary,line",
					quiet="true")

grass.run_command("g.remove",
					vect="circulo_dem_area",
					quiet="true")

grass.run_command("v.category",
					input="temporal",
					output="temporal_",
					option="add",
					cat=1,
					quiet="true")

grass.run_command("v.to.rast",
					input="temporal_",
					output="temporal_raster",
					use="cat",
					quiet="true")

texto=grass.read_command("r.report",
							flags="hn",
							map="temporal_raster",
							units="me")

texto=texto.split("|")

texto=(texto[12])

texto=texto.split(",")

sum_texto=""

for i in range(len(texto)):

	sum_texto=sum_texto + texto[i]

area_perimetro_sub=float(sum_texto)/2

area_total = area_sub + area_perimetro_sub

print "\narea total 3D = %.4f" % area_total, "m2"

grass.run_command("g.remove",
					vect="temporal",
					quiet="true")

grass.run_command("g.remove",
					vect="temporal_",
					quiet="true")

grass.run_command("g.remove",
					rast="temporal_raster",
					quiet="true")

La ejecución del código anterior, en cónsola de GRASS-QGIS, produjo el siguiente resultado:

Este resultado de 842405.0600 m2 difiere en un 0,41 % del valor obtenido por consideraciones geométricas (845899,709823203 m2).

Esta entrada fue publicada en Código Python, 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