Obtener valores rasterizados de una superposición de polígonos en soluciones SIG de código abierto

Question

Tengo dos capas. Una capa de forma poligonal con muchos azulejos y una capa de trama que contiene Cobertura terrestre de CORINE 2006 con muchas categorías en un mapa de colores. Quiero obtener para cada polígono del shapelayer una suma de cada categoría de cobertura del suelo del raster-layer.

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

• Tierras cultivables: 15 %.
• Bosque: 11 %
• Calles:2 % (... y así sucesivamente)

Intenté hacerlo en grass, qgis (sin función), saga (sólo suma cada uno a un valor total) r(suma total), pero aún no encontré solución. La mayoría de los plugins (estadísticas zonales en qgis) sólo soportan 0-1 capas raster. v.rast.stats tampoco ayudó. Estoy abierto a cualquier solución buena e inteligente. Tal vez he utilizado un enfoque equivocado o cometido errores.

En Arcgis esta tarea es bastante fácil, si no recuerdo mal, pero sigo echando en falta una buena solución para el usuario cotidiano de linux.

Estoy usando un sistema linux debian y por eso solo puedo usar programas para este sistema operativo.

---

EDITAR: Debido a que esta pregunta todavía tiene muchos puntos de vista y los visitantes: Escribí un QGIS-plugin, que también es capaz de calcular la landcover de la capa raster. No he codificado una superposición de polígono todavía, pero definitivamente planeado. Encuentra el plugin aquí e instalar la biblioteca Scipy primero.

Answer 1

Utilice 'extract' para superponer las características de los polígonos de un SpatialPolygonsDataFrame (que puede ser leído desde un shapefile utilizando maptools:readShapeSpatial) en un raster, y luego utilice '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í, mi primer polígono cubre 542 píxeles, y el segundo 958. Puedo resumir cada uno de ellos:

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

 26  27 
287 255 

[[2]]

  2  11  18 
 67  99 792 

Así que mi primer polígono tiene 287 píxeles de clase 26, y 255 píxeles de clase 27. Es bastante fácil sumar, dividir y multiplicar por 100 para obtener los porcentajes.

Answer 2

Quería informar y aquí estoy. La solución de Spacedman funcionó muy bien y pude exportar toda la información de cada polígono de mi forma. Sólo en caso de que alguien tiene el mismo problema, aquí es cómo he precedido:

...
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)
}

Answer 3

si entiendo bien lo que quieres, y asumiendo que tienes la capa vectorial 'mypolygonlayer' y la capa raster 'corina' ya en tu base de datos de GRASS GIS:

Primero convertiría el vector a raster. El gato se asegurará de que usted tendrá un identificador único por polígono. Si tiene una columna con un identificador numérico único, puede utilizar esa columna en su lugar. El labelcolumn es opcional:

v.to.rast input=mi capa de polígonos layer=1 output=mis polígonos use=cat labelcolumn=NameMappingUnit

A continuación, ejecute r.stats para obtener sus estadísticas:

r.stats -a -l input=mispolígonos,corina separator=; output=/home/paulo/corinastats.csv

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

Answer 4

¿Qué le parece convertir los datos de CORINE en un conjunto de datos poligonales vectoriales utilizando QGIS ( Raster > Conversión > Poligonización ) y luego utilizando la función Unión ( Vectorial > Herramientas de geoprocesamiento > Unión ) para combinarlos con los polígonos. El conjunto de datos vectoriales resultante contendría las áreas de cada clase de CORINE en cada polígono.

Answer 5

QGIS.

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

Esto lleva a cabo la función que usted requiere.

En cuanto al flujo de trabajo, no tengo claro cuántas tramas tienes o son sólo bandas en una trama?