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

14

Tengo dos capas. Una capa de forma de polígono con muchos mosaicos y una capa ráster que contiene Cobertura de suelo CORINE 2006 con muchas categorías en un mapa de colores. Quiero obtener para cada polígono en la capa de formas una suma de cada categoría de cobertura terrestre 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 en metros cuadrados):

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

Intenté hacerlo en grass, qgis (sin función), saga (solo resume todo a un valor total) r (suma total), pero todavía no encontré 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 equivocado o cometí errores.

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

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

EDITAR: Porque esta pregunta todavía tiene tantas visitas y visitantes: Escribí un complemento de QGIS, que también es capaz de calcular la cobertura terrestre de la capa ráster. No he codificado una superposición de polígonos todavía, pero definitivamente planeado. Encuentre el complemento aquí e instale la biblioteca Scipy primero.

    
pregunta Curlew 17.04.2012 - 13:08

7 respuestas

12

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.

    
respondido por el Spacedman 17.04.2012 - 15:59
5

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)
}
    
respondido por el Curlew 19.04.2012 - 16:52
3

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

    
respondido por el Ecodiv 08.10.2012 - 20:45
3

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.

respondido por el Emily 15.03.2016 - 11:51
2

¿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.

    
respondido por el Andy Harfoot 17.04.2012 - 16:05
0

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.

    
respondido por el Willy 18.04.2012 - 11:54
0

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.

    
respondido por el joaoal 13.06.2017 - 18:17

Lea otras preguntas en las etiquetas