Obtener valores ráster de una superposición de polígonos en soluciones GIS de código abierto

Use 'extraer' para superponer entidades poligonales de un SpatialPolygonsDataFrame (que puede leerse desde un shapefile usando maptools: readShapeSpatial) en un raster, luego use 'table' 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 ...

Así que 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. Lo 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 de cada polígono en mi forma. En caso de que alguien tenga el mismo problema, aquí es cómo 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 tiene la capa vectorial 'mypolygonlayer' y la capa raster 'corina' ya en su base de datos de GRASS GIS:

Primero convertiría vector a raster. 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 label 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, separador de corina =; output = / 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 he llevado a cabo el mismo tipo de análisis, pero descargar programas como R es un poco complicado en la computadora de mi trabajo y necesita aprobación. Después de muchas horas de investigación sobre un método que podía usar solo con QGis y Excel, encontré que este método funcionaba mejor para mí y quería ofrecerlo a las personas en el mismo tipo de situación.

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

2. Poligonizar la capa del recortador (Raster → Conversión → polygonise: archivo de entrada = su capa de recorte, guardar salida → ok)

3. Cálculo del número de píxeles (Haga clic en el archivo de formas que acaba de crear → abrir calculadora de campo: 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 poligonización (haga clic con el botón derecho en la capa de poligonal, guardar como: formato = archivo DBF, guardar como → ok)

5. Resumiendo el número de píxeles para cada hábitat (inicie Excel, abra el archivo, si no tiene títulos, agregue uno ahora para cada columna, haga clic en una celda vacía, acceda a DATOS tab, consolidar, asegúrese de que esté en la suma, haga clic en la flecha roja junto a "navegar ..." y seleccione sus dos columnas (títulos incluidos), haga clic en "agregar" y marque las casillas "Fila superior" y "Columna izquierda" → ok)

6. Si, como yo, tienes muchos polígonos para analizar y necesitas 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í del 1 al 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 (INDICE (C: C; MATCH (A2; B: B; 0)); "") y arrastre (o haga doble clic en la esquina inferior derecha) hacia abajo a su último valor (por lo tanto, si tiene 48 hábitats hasta la celda D48). El número de píxeles se encuentra ahora en las celdas correspondientes a su hábitat y puede repetir este proceso para todos sus polígonos.

¿Qué hay de convertir los datos CORINE en un conjunto de datos de polígonos vectoriales usando QGIS ( Raster > Conversion & Polygonize ) y luego usar la función Union ( Vector > Geoprocessing Tools > Union ) para combinar con los polígonos. El conjunto de datos vectoriales resultante contendría las áreas de cada clase CORINE en cada polígono.

QGIS.

En el troncal QGIS, hay otra versión de ZonalStats disponible, se llama Estadísticas Zonales.

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, argumentarí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ígono. Esto requiere mucho menos procesamiento y, por lo tanto, es la única solución fácil de implementar para procesar grandes rásteres y archivos grandes de polígonos en R.

Después de rasterizar sus polígonos a la misma extensión y resolución de los datos ráster, puede tabular las estadísticas de resumen como se explica here , lo que es apropiado si su raster encaja en la memoria (capas raster pequeñas / medianas) o puede binarizar cada categoría con la función reclass y luego calcular zonal de estadísticas para cada clase. Aquí es 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.