¿Cómo encontrar el valor promedio de ráster de un área definida por un shapefile usando R?

19

Tengo un conjunto de imágenes ráster que representan un mes específico a lo largo de los años, y quiero hacer una línea de tiempo de los valores promedio de un área usando un archivo de forma.

¿Cómo extraigo los valores de los rásteres y los importo en R de manera que pueda usarlos?

nickves
fuente

Respuestas:

23

Aquí hay un código de ejemplo. Es bastante sencillo adaptar este código para que funcione en un bucle para procesar todos sus rásteres. Si sus rásteres comparten una extensión y resolución comunes, puede crear una pila de ráster y recorrer las bandas de la pila. Para crear un vector que contenga todos los rásteres en un directorio, en un formato específico, puede usar "list.files" y luego pasar este vector a la pila.

Ejemplo:

rlist=list.files(getwd(), pattern="img$", full.names=TRUE) 

r <- stack(rlist)   


    # Add required libraries
    require(raster)
    require(sp)
    require(rgdal)

    # Set working directory, raster, in and out shapefiles
    setwd("C:/test")
    inshp="MyPolys"
    outshp="PolyMeans"
    rdata <- "Year2012.img"

    # Read polygon feature class shapefile
    sdata <- readOGR(dsn=getwd(), layer=inshp)

    # Read raster
    r <- raster(rdata)

    # Extract raster values to list object
    r.vals <- extract(r, sdata)

    # Use list apply to calculate mean for each polygon
    r.mean <- lapply(r.vals, FUN=mean)

    # Join mean values to polygon data
    sdata@data <- data.frame(sdata@data, m2012=r.mean)

    # Write results
    writeOGR(sdata, getwd(), outshp, driver="ESRI Shapefile", check_exists=TRUE, 
         overwrite_layer=TRUE)
Jeffrey Evans
fuente
1
+1 - excelente solución, bien diseñada y respuesta completa.
Simbamangu
Si tengo un data.frame(sdata@data, m2012=r.mean)archivo multipolígono: ¿cómo sabe qué polígono asignar a qué valor?
Stophface
para obtener r.mean agregado correctamente a sdata, primero tuve que anular la lista de r.mean: r.mean <- unlist (lapply (r.vals, FUN = mean))
cmbarbu el
6

Lea el archivo de forma en una SpatialPolygonsDataFrame( readOGRfunción del paquete rgdal)

Leer el ráster en un Rasterobjeto ( rasterfunción del paquete raster)

Use extract(raster, spdf)para obtener las celdas de la cuadrícula debajo de cada polígono. Entonces corre meansobre ellos.

Repita sobre su conjunto de imágenes ráster ...

Hombre espacial
fuente
¿Cómo escribir el archivo de forma (con los valores extraídos de la imagen ráster)?
Stophface