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

En el artículo anterior se calculó el área 3D de una montaña o elevación sumando al valor que determina el módulo r.surf.area lo que éste subestima al sólo considerar la mitad del área en la zona del perímetro. Cómo el procedimiento es relativamente largo, se automatizó mediante el siguiente script de bash:

#!/bin/bash
clear
read -p "Nombre del raster para calcular area 3D = " input_raster
temp1=`r.surf.area input=$input_raster vscale=1 | grep -i 'estimated' | grep -oE '[0-9.+e]*'`
temp2=`echo $temp1 | awk '{print $5}'`
temp3=`/usr/bin/printf "%f\n" $temp2`
temp4=`echo ${temp3/,/.}`
r.mapcalc "temporal=if($input_raster,1)"
r.to.vect --quiet input=temporal output=temporal_vector feature=area -s 
g.remove rast=temporal
v.type --quiet input=temporal_vector output=temporal_vector_linea type=boundary,line
g.remove vect=temporal_vector
v.category --quiet input=temporal_vector_linea type=line layer=1 output=temporal_vector_linea_cat
g.remove vect=temporal_vector_linea
v.to.rast --quiet input=temporal_vector_linea_cat value=1 output=temporal_linea_raster use=val
g.remove vect=temporal_vector_linea_cat
temp5=`r.report map=temporal_linea_raster units=me | grep -i 'value 1' | grep -oE '[0-9,]*'`
temp6=`echo $temp5 | awk '{print $3}'`
temp7=`echo ${temp6//,/""}`
suma=$(echo "scale=5; ($temp4 + $temp7 / 2)"| bc)
g.remove rast=temporal_linea_raster
echo "El area 3D es = " $suma "m2"

El área 3D total (834.935,8 m2) que refleja la siguiente salida producto de la ejecución del script:

Nombre del raster para calcular area 3D = circulo2
 100%
Removing raster <temporal>
Removing vector <temporal_vector>
Removing vector <temporal_vector_linea>
Removing vector <temporal_vector_linea_cat>
 100%
Removing raster <temporal_linea_raster>
El area 3D es =  834935.800000 m2

real    0m8.611s
user    0m4.488s
sys     0m0.076s

es idéntica a la obtenida en el artículo anterior. La ejecución del script apenas tomó poco más de 8 segundos (incluyendo los 4 segundos empleados en la escritura del nombre del ráster, circulo2, objeto del cálculo del área 3D).

Nota 1: Se agradece la colaboración del usuario neurus (Espacio Linux) quien sugirió como cambiar el formato numérico para que finalmente funcionara la calculadora de precisión arbitraria bc dentro del script.

Nota 2: Se detectó un bug en el script, ejecutándolo para elevacion_500_dem, el cual ya fue corregido. La salida tenía que ser esta:

Nombre del raster para calcular area 3D = elevacion_500_dem
 100%
Removing raster <temporal>
Removing vector <temporal_vector>
Removing vector <temporal_vector_linea>
Removing vector <temporal_vector_linea_cat>
 100%
Removing raster <temporal_linea_raster>
El area 3D es =  18094391.000000 m2

real    0m15.741s
user    0m4.316s
sys     0m0.088s
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