Obtenga valores ráster de una superposición de polígonos en las soluciones de SIG de Opensource

Tengo dos capas Una capa con forma de polígono con muchos mosaicos y una capa ráster que contiene la cobertura terrestre CORINE 2006 con muchas categorías en un mapa de colores. Quiero obtener para cada polígono en la capa de forma una suma de cada categoría de cobertura del suelo de la capa ráster.

Por ejemplo, hay un polígono con id '2' y quiero atributos como este para este polígono (en porcentaje o metros cuadrados):

• Tierra arable: 15%
• Bosque: 11%
• Calles: 2% (... y uno)

Traté de hacerlo en hierba, qgis (sin función), saga (solo suma cada uno a un valor total) r (suma total), pero todavía no encontré ninguna solución. La mayoría de los complementos (estadísticas zonales en qgis) solo admiten 0-1 capas ráster. v.rast.stats tampoco ayudó. ¡Estoy abierto a cualquier solución buena e inteligente! Tal vez incluso utilicé un enfoque incorrecto o cometí errores.

En Arcgis esta tarea es bastante fácil, si recuerdo bien, pero todavía me falta una buena solución para su usuario de Linux diario.

Estoy ejecutando un sistema Debian Linux y es por eso que solo puedo usar programas para este sistema operativo.

EDITAR: Debido a que esta pregunta todavía tiene muchas vistas y visitantes: escribí un complemento QGIS, que también es capaz de calcular la cobertura de la capa ráster. Todavía no he codificado una superposición de polígonos, pero definitivamente está planificada. Encuentre el complemento aquí e instale primero la biblioteca Scipy.

Use 'extraer' para superponer características de polígono de un SpatialPolygonsDataFrame (que se puede leer desde un archivo de forma usando maptools: readShapeSpatial) en un ráster, luego use 'tabla' para resumir. Ejemplo:

> c=raster("cumbria.tif") # this is my CORINE land use raster
> summary(spd)
Object of class SpatialPolygonsDataFrame
[...]
> nrow(spd)  # how many polygons:
[1] 2
> ovR = extract(c,spd)
> str(ovR)
List of 2
 $ : num [1:542] 26 26 26 26 26 26 26 26 26 26 ...
 $ : num [1:958] 18 18 18 18 18 18 18 18 18 18 ...

Entonces, mi primer polígono cubre 542 píxeles, y mi segundo cubre 958. Puedo resumir cada uno de ellos:

> lapply(ovR,table)
[[1]]

 26  27 
287 255 

[[2]]

  2  11  18 
 67  99 792 

Mi primer polígono es 287 píxeles de la clase 26 y 255 píxeles de la clase 27. Suficientemente fácil de sumar y dividir y multiplicar por 100 para obtener porcentajes.

Quería informar de nuevo y aquí estoy. La solución de Spacedman funcionó muy bien y pude exportar toda la información para cada polígono en mi forma. En caso de que alguien tenga el mismo problema, así es como lo precedí:

...
tab <- apply(ovR,table)
# Calculate percentage of landcover types for each polygon-field.
# landcover is a datastream with the names of every polygon
for(i in 1:length(tab)){
 s <- sum(tab[[i]])
 mat <- as.matrix(tab[[i]])
 landcover[i,paste("X",row.names(mat),sep="")] <- as.numeric(tab[[i]]/s)
}

si entiendo correctamente lo que quiere, y suponiendo que ya tenga la capa vectorial 'mypolygonlayer' y la capa ráster 'corina' en su base de datos GRASS GIS:

Primero convertiría el vector en ráster. El gato se asegurará de que tengas un identificador único por polígono. Si tiene una columna con un identificador numérico único, puede usar esa columna en su lugar. La columna de etiqueta es opcional:

v.to.rast input = mypolygonlayer layer = 1 output = mypolygons use = cat labelcolumn = NameMappingUnit

Luego ejecute r.stats para obtener sus estadísticas:

r.stats -a -l input = mypolygons, corina separator =; salida = / home / paulo / corinastats.csv

El último paso es abrir corinastats.csv en, por ejemplo, LibreOffice y crear una tabla dinámica o usar R para crear su tabla cruzada

Sé que esta publicación es bastante antigua, pero recientemente me gustaría llevar a cabo el mismo tipo de análisis, pero descargar programas como R es un poco complicado en mi computadora de trabajo y necesita aprobación. Después de muchas horas de investigar un método que podía usar solo con QGis y Excel, descubrí que este método funcionaba mejor para mí y quería ofrecerlo a personas en el mismo tipo de situación.

1. Recorte el polígono a la capa ráster (Ráster → extracción → clipper: archivo de entrada = capa ráster, elija su nombre y ubicación de salida, haga clic en la capa de máscara, elija su polígono → ok)

2. Poligonalice la capa de recortadora (Ráster → Conversión → poligonalización: archivo de entrada = su capa de clip, guarde salida → ok)

3. Cálculo del número de píxeles (haga clic en el archivo de forma que acaba de crear → calculadora de campo abierto: marque “crear nuevo campo” y agregue el nombre del campo, Función = geometría → área → ok). Ahora debería tener una nueva columna en su tabla de atributos que muestre la cantidad de píxeles.

4. Guardar capa de poligonalización (haga clic con el botón derecho en la capa de poligonalización, guardar como: formato = archivo DBF, guardar como → ok)

5. Resumiendo el número de píxeles para cada hábitat (comience Excel, abra el archivo, si no tiene títulos agregue uno ahora para cada columna, haga clic en una celda vacía, vaya a la pestaña DATOS, consolide, asegúrese de que esté en la suma, haga clic en el flecha roja al lado de "navegar ..." y seleccione dos columnas (títulos incluidos), haga clic en "agregar" y marque los cuadros "Fila superior" y "columna izquierda" → ok)

6. Si, como yo, tiene muchos polígonos para analizar y necesita compararlos en la misma tabla, este paso será útil. En un nuevo libro de Excel, enumere los números de sus hábitats en la columna A (para mí 1 a 48) y coloque las dos columnas que acaba de consolidar en las columnas B y C (hábitat en B y número de píxeles en C). En la celda D1 escriba la siguiente fórmula: = IFNA (ÍNDICE (C: C; COINCIDIR (A2; B: B; 0)); "") y arrastre (o haga doble clic en la esquina inferior derecha) hasta su último valor (por lo que si tienes 48 hábitats hasta la celda D48). El número de píxeles está ahora en las celdas correspondientes a su hábitat y puede repetir este proceso para todos sus polígonos.

¿Qué tal convertir los datos CORINE en un conjunto de datos de polígonos vectoriales usando QGIS ( Ráster> Conversión> Poligonalizar ) y luego usar la función Unión ( Vector> Herramientas de geoprocesamiento> Unión ) para combinar con los polígonos. El conjunto de datos del vector resultante contendría las áreas de cada clase CORINE en cada polígono.

QGIS.

En el tronco QGIS, hay otra versión de ZonalStats disponible, que se llama Zonal Statistics.

Esto lleva a cabo la función que necesita.

En cuanto al flujo de trabajo, no tengo claro cuántos rásteres tiene o son solo bandas en un ráster.

Opuesto a la mayoría de las respuestas anteriores, diría que la mejor opción es rasterizar sus polígonos y trabajar con dos conjuntos de datos ráster en lugar de dos conjuntos de datos de polígonos. Esto requiere mucho menos procesamiento y, en consecuencia, es la única solución fácil de implementar para procesar grandes rásteres y grandes archivos de polígonos en R.

Después de rasterizar sus polígonos exactamente con la misma extensión y resolución de los datos ráster, puede tabular estadísticas de resumen como se explica aquí , lo cual es apropiado si su ráster cabe en la memoria (capas ráster pequeñas / medianas) o puede binarizar cada categoría con la reclassfunción y luego calcular zonalestadísticas para cada clase. Aquí hay una solución que incorpora la rasterización y las estadísticas zonales en una función y funciona bien con conjuntos de datos muy grandes.