¿Cómo se mide el volumen en la trama?

Question

Hay muchas preguntas similares como esta. Pero aún no he encontrado una respuesta útil. Así que espero que algunos de ustedes tengan una sugerencia para esto.

He hecho algunos sobrevuelos con el dron en un área de construcción, y luego saqué un DSM de estos, usando Pix4D. Ahora tengo dos DSM diferentes, y luego los importé a QGIS. A través de la rastercalculadora, los he retirado uno del otro para obtener la diferencia de altura. Ver abajo:

Entonces me gustaría saber el volumen de donde ha habido una diferencia de altura.. Así que parte del área que he marcado en rojo (aumento de volumen) y en azul (disminución de volumen). También usando la rastercalculadora y la poligonización. Ver abajo.

¿Hay una forma inteligente de encontrar el volumen de las marcas azules y rojas de estas capas de trama en QGIS, o alguno de ustedes tiene un programa alternativo que haría el trabajo más fácil y rápido?

Answer 1

La mejor manera que he encontrado para hacerlo es usar una combinación de la calculadora de trama y las estadísticas zonales.

Primero hay que calcular la diferencia entre los dos rasters, usando un filtro para obtener la diferencia positiva o negativa, y ningún valor en otra parte. En la calculadora de Qgis esto daría algo como:

((A-B)>0)*(A-B) Para conseguir el cambio positivo

((A-B)<0)*(A-B)*-1 Para obtener el cambio negativo (multiplicar por -1 si te interesa el valor absoluto del cambio)

Una vez hecho esto, cree una nueva capa de vectores, mostrando la zona de interés (azul + rojo, sólo evite los lugares donde no le interesa. Como esto es una construcción, probablemente pueda seleccionar sólo las zonas donde sabe que ha habido cambios).

Entonces usa la herramienta en Raster -> Estadísticas zonales para leer la suma de los valores de todos los píxeles de la capa vectorial. Hágalo dos veces, tanto para el cambio positivo como para el negativo.

Multiplica entonces el valor obtenido para cada caso con el tamaño de tu píxel al cuadrado, para obtener el volumen total.

Answer 2

Una forma muy rápida de obtener el volumen con un código de pitón y gdal2xyz . Tienes que calcular tu diferencia sólo raster. No se preocupe por el valor positivo o negativo, lo distinguiremos más tarde.

1. Convierte tu archivo raster de diferencias en un archivo csv con el OSGeo4W Shell:

   gdal2xyz your_raster.tif your_raster.csv

2. Crear una matriz numérica a partir de este archivo csv y crear dos subconjuntos:

   import numpy as np
   x = np.genfromtxt (r'your_raster.csv', delimiter=" ", usecols=range(2,3))
   pos = x[x>0.2]
   neg = x[x<-0.2]
   pos.sum(axis=0)
   neg.sum(axis=0)

Explicación del código Python:

El archivo csv normalmente tiene 3 columnas x,y,z. Con usecols= range(2,3) se importa sólo la 3ª columna (valores de diferencia).

Con pos = x[x>0.2] guardas todos los valores positivos >0.2 en un nuevo arreglo.

Después puedes sumar los dos subconjuntos. Utiliza pos.sum(axis=0) para obtener la suma de todos los valores positivos >0,2.

Si tu mapa de bits tiene valores NODATA (-99999 por ejemplo) puedes eliminarlos de tu matriz con:

index = np.argwhere(x==-99999)
x = np.delete(x, index)

Entonces puedes ir más allá con la creación de los subconjuntos.

También hay una forma de crear una matriz numérica a partir de archivos raster directamente. Obtengo resultados ligeramente diferentes con este enfoque:

import numpy as np
from osgeo import gdal
ds = gdal.Open("your_raster.tif")
x = np.array(ds.GetRasterBand(1).ReadAsArray())
pos = x[x>0]
neg = x[x<0]
pos.sum(axis=0)
neg.sum(axis=0)

Finalmente tienes que multiplicar tu suma con el tamaño de la celda que tienes.