Calcular media aritmética por columnas en un ráster usando un script de bash

En un artículo anterior se calculó la media aritmética por columnas en un ráster usando QGIS y lenguaje R. En el presente artículo se va a calcular la media aritmética y la desviación estándar por columnas en un ráster pero usando un script de bash en ambiente de GRASS-QGIS. El script usa g.región para “cortar” la imagen ráster original en todas las posibles imágenes ráster con espesor de una columna, las vectoriza como puntos y usa la tabla atributiva de estas para, mediante la adaptación de un script anterior, determinar las medias aritméticas y desviaciones estándar por columnas en el ráster.

El script es el siguiente:

#! /bin/bash
clear
array=`g.region -g|grep -oE '[0-9.]*'`
norte=`echo $array|awk '{print $1}'`
sur=`echo $array|awk '{print $2}'`
oeste=`echo $array|awk '{print $3}'`
este=`echo $array|awk '{print $4}'`
eores=`echo $array|awk '{print $6}'`
cols=`echo $array|awk '{print $8}'`
oeste=`echo $(echo "scale=8; $oeste + 0.00000000" | bc )`
este=`echo $(echo "scale=8; $oeste + $eores" | bc )`
g.region w=$oeste e=$este
r.mapcalc raster=N09W068_UTM19N_canoa 
r.to.vect input=raster output=raster_puntos feature=point -s --quiet
#Calculo del promedio y desviacion estándar
temporal=`v.db.select map=raster_puntos col=value | grep -oE '[0-9.-]*'`
suma=0
suma2=0
iter=0
for var in $temporal
do 
	suma=$(echo "$suma + $var" | bc )
	suma2=$(echo "$suma2 + $var*$var" | bc )
	iter=$(expr $iter + 1)
done
promedio=$(echo "scale=5; $suma/$iter" |bc)
sdev=$(echo "scale=5; sqrt (($suma2-$suma*$suma/$iter)/($iter-1))" | bc -l)
echo $promedio $sdev >> estadistica
g.remove vect=raster_puntos --quiet
#for i in $(seq 1 $(expr $cols - 1 )) 
for i in $(seq 1 5) #Para probar el script con solo 6 valores pareados: tiene 3603!
   do
      oeste=`echo $(echo "scale=8; $oeste + $eores" | bc )`
      este=`echo $(echo "scale=8; $oeste + $eores" | bc )`
      g.region w=$oeste e=$este
      r.mapcalc raster=N09W068_UTM19N_canoa 
      r.to.vect input=raster output=raster_puntos feature=point -s --quiet
      #Calculo del promedio y desviacion estándar
      temporal=`v.db.select map=raster_puntos col=value | grep -oE '[0-9.-]*'`
      suma=0
      suma2=0
      iter=0
      for var in $temporal
      do 
         suma=$(echo "$suma + $var" | bc )
	 suma2=$(echo "$suma2 + $var*$var" | bc )
	 iter=$(expr $iter + 1)
      done
      promedio=$(echo "scale=5; $suma/$iter" |bc)
      sdev=$(echo "scale=5; sqrt (($suma2-$suma*$suma/$iter)/($iter-1))" | bc -l)
      echo $promedio $sdev >> estadistica
      g.remove vect=raster_puntos --quiet
done
g.region rast=N09W068_UTM19N_canoa

Se grabó como script_media_raster, se le dio permisos de ejecución con chmod +x script_media_raster y se ejecutó con ./script_media_raster. La salida está en el archivo estadistica que, por razones de ejecución, sólo contiene los primeros 6 valores pareados (de los 3603 posibles en este caso). Para producirlos todos hay que descomentar en el código a continuación (correspondiente al script anterior) la línea superior y eliminar la inferior:

.
.
.
#for i in $(seq 1 $(expr $cols - 1 )) 
for i in $(seq 1 5) #Para probar el script con solo 6 valores pareados: tiene 3603!
.
.
.

Para este caso, el archivo estadística (pareados medias y desviación estándar) tiene estos valores:

25.78429 118.11121
86.05812 218.15340
126.67603 240.70155
144.53085 236.86797
163.83057 232.39561
176.89008 225.83195

cuyas medias son idénticas a las 6 primeras producidas por el procedimiento encontrado en:

Calcular media aritmética por columnas en un ráster usando QGIS y lenguaje R

cuyo extracto del archivo (contiene sólo medias) es el siguiente:

25.7843 86.05813 126.6760 144.5309 163.8306
176.8901 188.0510 197.6449 205.6521 212.4804
218.9931 218.7366 218.5802 218.4242 218.2702
218.1824 218.2399 218.3791 218.5138 218.6488
218.8008 218.9394 219.0609 219.1474 219.1975
219.1857 219.1617 219.1245 219.0612 219.0631
219.1758 219.4176 219.6727 219.8837 220.0719
220.2551 220.4085 220.5099 220.6240 220.7601
220.9402 221.0989 221.1543 221.1738 221.1727
221.1570 221.2063 221.4386 221.7298 222.0634
222.4113 222.7441 222.9592 223.2025 223.4501
223.6766 223.9474 224.2028 224.3325 224.4050
224.4862 224.6240 224.7854 225.0179 225.2884
225.5449 225.7815 226.0242 226.2804 226.5496
226.7799 226.9364 227.0204 226.9876 226.8928
226.7884 226.727 226.5959 226.3300 226.0496
225.8534 225.7474 225.5857 225.3904 225.2168
.
.
.

Esta concordancia parece señalar que el script trabaja de manera adecuada.

Nota: Como tenía curiosidad por saber cuánto tiempo tardaría el script para procesar las 3603 columnas (para medias y desviaciones estándares de 3630 filas), lo ejecuté con time. El tiempo de ejecución fue de 15 horas, 51 minutos, 18.411 segundos. Todos los valores obtenidos (7206) concordaron de la manera adecuada.

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

2 respuestas a Calcular media aritmética por columnas en un ráster usando un script de bash

  1. Pingback: Calcular media aritmética por columnas en un ráster usando un script de bash | Geoprocessing - Geoprocessamento | Scoop.it

  2. Pingback: Calcular media aritmética por columnas en un ráster usando un script de python |

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