Ejecución de un script de Python para extraer datos raster en QGIS

Question

Estoy intentando utilizar QGIS para extraer estadísticas rasterizadas de los shapefiles (tengo shapefiles de distribución de especies y quiero extraer datos ambientales del área de distribución de cada especie). Esto es posible utilizando el plugin Zonal Statistics, pero tengo unos 300 shapefiles de los que necesito obtener datos. He intentado fusionar los shapefiles en uno solo, pero esto da lugar a archivos masivos que hacen que QGIS se bloquee (el uso de shapefiles de más de 15mB causa problemas).

¿Hay alguna manera de automatizar el procedimiento utilizando un script de Python? He encontrado un post que detalla cómo ejecutar ZonalStats usando Python ( ¿Cómo se calculan las estadísticas ráster de los polígonos? segunda respuesta) pero no consigo que funcione. No sé absolutamente nada acerca de Python y por lo tanto no estoy seguro de si hay algo que hay que hacer antes de intentar ejecutar el script o si el problema es otra cosa. Un script similar que se ejecute por sí mismo sin la entrada del usuario sería ideal.

Estoy ejecutando QGIS 1.8 en una máquina Linux.

Cualquier ayuda será muy apreciada.

Answer 1

Sí se puede hacer esto en un script. Para hacer esto de manera eficiente, usted necesitará SciPy . SciPy es un primo de NumPy y juntos amplían Python con un montón de capacidades de manejo de arrays extremadamente potentes. En SciPy se utilizará la función ' ndimage Módulo '.

Su procedimiento es hacer un raster de sus zonas. Luego importa tu raster original y el raster de las zonas en un array NumPy. Haz una tercera matriz numpy con los valores únicos de tu matriz de zonas para usarla como índice para la función ndimage. Por último, utiliza una de las estadísticas de ndimage para extraer el resultado, por ejemplo

myIndex = numpy.unique(zonesArray)    
minResult = ndimage.minimum(myRasterArray, labels=zonesArray, index=indexArray)

(obviamente no he incluido todo el código para abrir un archivo raster o vectorial, etc., pero vea más abajo la ayuda para ello).

A menudo escribo mis resultados en un simple archivo de texto csv y luego puedo hacer una unión de tablas para extraer las estadísticas de nuevo en mis polígonos de zona. Si usted tiene muchos rasters, entonces usted puede automatizar el proceso en un proceso por lotes.

Aunque se puede hacer esto con GDAL y, de hecho, se necesitará GDAL para leer los datos originales de ráster y polígono, no recomendaría que se intentara sin SciPy porque el proceso tendrá una presteza glacial. Ver esta excelente serie sobre la codificación de GDAL con Python. Observe el tutorial 6 de la serie, ya que es especialmente relevante. Sin embargo, no menciona SciPy ya que es un poco antiguo (y donde menciona 'Numeric' entienda que éste ha sido reemplazado por 'NumPy'). ndimage es muy rápido y mucho mejor que iterar manualmente sobre un ráster de forma bit a bit o incluso iterar sobre un array celda a celda. La diferencia de velocidad al usar SciPy frente a NumPy puede ser de segundos comparado con horas para los mismos datos (por experiencia personal). Además, se necesita menos código :)

Este mensaje contiene muchas cosas y puede parecer desalentador si no sabes nada de Python. En ese caso, comience con la serie de tutoriales que mencioné anteriormente, que asume cero experiencia en codificación. Luego, después de las primeras lecciones, salta al material de rasterización. Una vez que sepas cómo abrir un raster y ponerlo en un array Numpy (también conocido como 'numérico'), lee la documentación de SciPy sobre ndimage .

Answer 2

En general: Mi QGIS (1.9 - Master ( ;-) ) en Debian Linux) funciona bastante bien con shapefiles grandes (diablos, tengo algunos que son más de 300mb de tamaño). Tal vez esto es un problema con su proceso de conversión / fusión y si es así usted debe cambiar a otras herramientas (buscar en este sitio web).

En cuanto a sus datos:

• No has especificado la naturaleza de tus "shapefiles". Si tiene una capa de puntos, podría utilizar la herramienta de muestreo de puntos en QGIS (búsquela a través del descargador de plugins)
• Si tienes un shapefile con polígonos y de alguna manera tienes experiencia con los análisis estadísticos en R, entonces mira este pequeño ejemplo de mina . Utilice el R-Code para extraer los valores de su trama.
• Si lo que quieres es estrictamente un código python, deberías buscar en el api gdal y codificarlo.
• También te animamos a que utilices otras herramientas muy buenas que hay por ahí (incluso para linux). Tanto GRASS como SAGA tienen una función que hace exactamente lo que quieres. Y si usted es inteligente, incluso podría utilizar esas funciones directamente en QGIS a través de la caja de herramientas Sextante.

Answer 3

La solución que buscaba ya ha sido respondida en otra pregunta que publiqué ( ¿Cómo utilizar el plugin QGIS Zonal Stats desde la consola de Python? ). Sin embargo, la solución real que utilicé fue en R, como sugirió aquí Curlew.

El código que utilicé fue el siguiente:

library(raster)
library(rgdal)
library(foreign)

# extract raster data from overlying polygon
rast <- raster('raster.tif')
poly <- readOGR(dsn="polygons", layer="polygon")
ext.poly <- extract(rast, poly, fun = mean, na.rm=TRUE, df=TRUE)

# and then append the result to the .dbf file of the original polygon
poly.dbf <- read.dbf("polygons/polygon.dbf", as.is = TRUE)
poly.dbf$result <- ext.poly$extraction
write.dbf(poly.dbf, "polygon_with_extraction.dbf")

Gracias a todos por la ayuda.