Como obter contagem de células raster não NA dentro do polígono usando R

Eu tenho encontrado todos os tipos de problemas usando o ArcGIS ZonalStats e achei que o R poderia ser uma ótima maneira. Dizendo que sou bastante novo no R, mas tenho um histórico de codificação.

A situação é que eu tenho vários rasters e um shapefile de polígono com muitos recursos de tamanhos diferentes (embora todos os recursos sejam maiores que uma célula raster e os recursos do polígono estejam alinhados ao raster). Eu descobri como obter o valor médio para cada recurso de polígono usando a biblioteca raster com extração:

#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'])

O problema que tenho é que também preciso de uma proporção baseada em área entre a área do polígono e todas as células não NA no mesmo polígono. Sei qual é o tamanho da célula da varredura e posso obter a área para cada polígono, mas o link ausente é a contagem de todas as células não NA em cada recurso. Consegui obter o número de células de todas as células no polígono proxNA@data$Cnumb1000 <- cellFromPolygon(ras, proxNA)e tenho certeza de que existe uma maneira de obter o valor real da célula de varredura, que exige um loop para obter o número de todas as células não NA combinadas com uma contagem, etc. MAS, tenho certeza que existe uma maneira muito melhor e mais rápida de fazer isso! Se algum de vocês tiver uma idéia ou puder me apontar na direção certa, ficaria muito grato!

Os dados de exemplo 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"))))

Agora use extrato

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

Se você tiver muitos rasters, primeiro use a pilha para combiná-los (se eles tiverem a mesma extensão e resolução)

Para obter a área real do polígono, você não deve usar a abordagem de slot (i, 'área'). Para dados planares, você pode usar rgeos :: gArea (polys, byid = TRUE) Para dados esféricos (lon / lat), você pode usar geosphere :: areaPolygon

Não tenho certeza se você deseja a proporção com base no "valor real" das áreas de polígono (s) ou nas áreas das células que as interceptam. Aqui está um exemplo de código que usa todas as células que interceptam os polígonos (basicamente, proporção de células NA para células não NA). É um exemplo fictício e você precisará escrever sua própria função.

    # 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)

Você pode usar esse trecho de código para adicionar uma coluna de área aos dados do polígono, extraindo do slot de área. Isso pode ser usado se você deseja propor uma proporção usando a área de polígono "real".

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

A alternativa mais eficiente e consideravelmente mais rápida aos loops for "aplicar" como funções. Há vários disponíveis no R que são utilizados para diferentes classes de objetos ou estruturas de dados. Nesse caso, como extract retorna uma lista, usaremos lapply (lista aplicável). Essa é uma maneira de aplicar uma função básica ou personalizada a um objeto de lista. A classe de objeto armazenada na lista é um vetor, a função é bastante direta. Se você usar extração em um objeto de varredura de tijolo ou pilha, os objetos resultantes armazenados na lista seriam 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)         

Não tenho acesso aos seus arquivos, mas com base no que você descreveu, isso deve 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)