Cómo obtener el recuento de celdas ráster no NA dentro del polígono usando R

8

Me he encontrado con todo tipo de problemas usando ArcGIS ZonalStats y pensé que R podría ser una excelente manera. Diciendo que soy bastante nuevo en R, pero tengo un fondo de codificación.

La situación es que tengo varios rásteres y un archivo de forma de polígono con muchas características de diferentes tamaños (aunque todas las características son más grandes que una celda de trama y las características de polígono están alineadas con el ráster). He descubierto cómo obtener el valor medio para cada entidad poligonal usando la biblioteca ráster con extracto:

#load packages required
require(rgdal)
require(sp)
require(raster)
# ---Set the working directory-------
datdir <- "/test_data/"

#Read in grid of water depth
ras <- raster("test_data/raster/pl_sm_rp1000/w001001.adf")

#read in polygon shape file
proxNA <- shapefile("test_data/proxy/PL_proxy_WD_NA_test") 

#calc mean depth per polygon feature
#unweighted - only assigns grid to district if centroid is in that district
proxNA$RP1000 <- extract(ras, proxNA, fun = mean, na.rm = TRUE, weights = FALSE)

#plot depth values 
spplot(proxNA[,'RP1000'])

El problema que tengo es que también necesito una relación basada en el área entre el área del polígono y todas las celdas que no son de NA en el mismo polígono. Sé cuál es el tamaño de celda del ráster y puedo obtener el área para cada polígono, pero el enlace que falta es el recuento de todas las celdas que no son NA en cada entidad. Logré obtener el número de celda de todas las celdas en el polígono proxNA@data$Cnumb1000 <- cellFromPolygon(ras, proxNA)y estoy seguro de que hay una manera de obtener el valor real de la celda ráster, que luego requiere un bucle para obtener el número de todas las celdas que no son NA combinadas con un recuento, etc. PERO, estoy seguro de que hay una manera mucho mejor y más rápida de hacerlo. Si alguno de ustedes tiene una idea o puede señalarme en la dirección correcta, ¡estaría muy agradecido!

Hubert
fuente
Si realmente hubiera terminado de depurar el enfoque zonalstats (que probablemente sería la forma ideal), vería numpy antes de R. Dicho esto, ¿se están rassosteniendo banderas de NA con un valor legítimo? Parece que podría filtrar ese valor u obtener un recuento de esos valores después del hecho.
Roland
@Roland: ¡Gracias! NA es NA en ras y no tiene un valor específico. Entonces, lo que está diciendo es que podría filtrar por NA (o valor de reemplazo) y clasificar para cada polígono para obtener el recuento, luego restar del número de celda general. Interesante, pero todavía un poco largo. Esperaba una función de conteo o algo por el estilo.
Hubert
Puede obtener un aumento de semillas al no utilizar el formato ráster. La ventaja de la trama es que es segura para la memoria. Como está creando un objeto sp y luego coaccionando a ráster, pierde la ventaja. Mantenerlo como un objeto sp y usar "over" será mucho más rápido que usar "extract". También procesará todo en la memoria.
Jeffrey Evans el

Respuestas:

5

Los datos de ejemplo de Jeffrey

library(raster)
r <- raster(ncols=10, nrows=10)
set.seed(0)
x <- runif(ncell(r))
x[round(runif(25,1,100),digits=0)] <- NA
r[] <- x
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1),  Polygons(list(Polygon(cds2)), 2)))
polys <- SpatialPolygonsDataFrame(polys, data.frame(ID=sapply(slot(polys, "polygons"), function(x) slot(x, "ID"))))

Ahora usa extracto

extract(r, polys, fun=function(x, ...) length(na.omit(x))/length(x))
#[1] 0.8333333 0.6666667

Si tiene muchos rásteres, primero use la pila para combinarlos (si tienen la misma extensión y resolución)

Para obtener el área del polígono real que debe no utilizar la ranura (i, 'zona') enfoque. Para datos planos puede usar rgeos :: gArea (polys, byid = TRUE) Para datos esféricos (lon / lat) puede usar geosphere :: areaPolygon

Robert Hijmans
fuente
3

No estoy seguro de si desea que la relación se base en el "valor real" de las áreas de polígono (s) o las áreas de las celdas que las intersectan. Aquí hay un código de ejemplo que usa todas las celdas que se cruzan con los polígonos (básicamente, relación de celdas NA a celdas que no son NA). Es un ejemplo ficticio y deberá escribir su propia función.

    # Create some example data
    require(raster)
    require(sp)

    r <- raster(ncols=10, nrows=10)
      x <- runif(ncell(r))
        x[round(runif(25,1,100),digits=0)] <- NA
          r[] <- x
      cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
        cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
          polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), 
                                   Polygons(list(Polygon(cds2)), 2)))
            polys <- SpatialPolygonsDataFrame(polys, data.frame(ID=sapply(slot(polys, "polygons"), 
                                              function(x) slot(x, "ID"))))
plot(r)
  plot(polys, add=TRUE)

Puede usar este fragmento de código para agregar una columna de área a los datos de su polígono mediante la extracción de la ranura de área. Esto podría usarse si desea hacer una relación usando el área poligonal "real".

# Add area of polygon(s)
polys@data <- data.frame(polys@data, Area=sapply(slot(polys, 'polygons'), 
                         function(i) slot(i, 'area')))  

La alternativa más eficiente y considerablemente más rápida a los bucles for son las funciones "aplicar". Hay varios de estos disponibles en R que se utilizan para diferentes clases de objetos o estructuras de datos. En este caso, dado que extract devuelve una lista, usaremos lapply (aplica la lista). Esta es una manera de aplicar una función base o personalizada a un objeto de lista. La clase de objeto almacenada en la lista es un vector, la función es bastante sencilla. Si usa el extracto en un objeto raster de ladrillo o pila, los objetos resultantes almacenados en la lista serían objetos data.frame.

# On a single raster object, extract returns list object with stored vectors.                           
( vList <- extract(r, polys, na.rm=FALSE) )
  class(vList)

# Use lapply to apply function that calculates ratio of NA to non-NA values
#   wrapping lapply in unlist() collapses result into a vector  
aRatio <- function(x) { if(length(x[is.na(x)]) > 0) (length(x[is.na(x)]) / length(x[!is.na(x)])) else 0 }  
  ( vArea <- unlist( lapply(vList, FUN=aRatio ) ) )

# Assign ordered vector back to polygons
polys@data <- data.frame(polys@data, NAratio=vArea)
  str(polys@data)         
Jeffrey Evans
fuente
Gracias Jeffrey! He aprendido bastante de tu respuesta. Pero creo que no me expliqué lo suficientemente bien. La relación que busco es Área de células no AN dentro de Poly1 a Área de Poly1. Algunos polígonos no están completamente cubiertos por celdas ráster. Escribir el valor medio de todas las celdas dentro de un polígono en vList es genial. Ahora solo necesito obtener el número de células NonNA a partir de las cuales se deriva la media, ya que conozco el área de cada célula. La relación se puede derivar fácilmente por (número de celdas * área de celda) / área de polígono. ¡Muchas gracias!
Hubert
1

No tengo acceso a sus archivos, pero según lo que describió, esto debería funcionar:

library(raster)
mask_layer=shapefile(paste0(shapedir,"AOI.shp"))
original_raster=raster(paste0(template_raster_dir,"temp_raster_DecDeg250.tif"))
nonNA_raster=!is.na(original_raster)
masked_img=mask(nonNA_raster,mask_layer) #based on centroid location of cells
nonNA_count=cellStats(masked_img, sum)
Lucas Fortini
fuente