Comment obtenir le nombre de cellules raster non NA dans un polygone à l'aide de R

J'ai rencontré toutes sortes de problèmes en utilisant ArcGIS ZonalStats et j'ai pensé que R pourrait être un excellent moyen. Dire que je suis assez nouveau pour R, mais j'ai un fond de codage.

La situation est que j'ai plusieurs rasters et un fichier de formes polygonales avec de nombreuses entités de tailles différentes (bien que toutes les entités soient plus grandes qu'une cellule raster et que les entités polygonales soient alignées sur le raster). J'ai trouvé comment obtenir la valeur moyenne de chaque entité surfacique en utilisant la bibliothèque raster avec extrait:

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

Le problème que j'ai est que j'ai également besoin d'un rapport basé sur l'aire entre l'aire du polygone et toutes les cellules non NA du même polygone. Je connais la taille des cellules du raster et je peux obtenir l'aire de chaque polygone, mais le lien manquant est le nombre de toutes les cellules non NA de chaque entité. J'ai réussi à obtenir le numéro de cellule de toutes les cellules du polygone proxNA@data$Cnumb1000 <- cellFromPolygon(ras, proxNA)et je suis sûr qu'il existe un moyen d'obtenir la valeur réelle de la cellule raster, qui nécessite ensuite une boucle pour obtenir le nombre de toutes les cellules non NA combinées avec un décompte, etc. MAIS, je suis sûr qu'il existe un moyen bien meilleur et plus rapide de le faire! Si l'un d'entre vous a une idée ou peut me diriger dans la bonne direction, je lui en serais très reconnaissant!

Les exemples de données 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"))))

Utilisez maintenant l'extrait

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

Si vous avez plusieurs rasters, utilisez d'abord la pile pour les combiner (s'ils ont la même étendue et la même résolution)

Pour obtenir la zone de polygone réelle, vous ne devez pas utiliser l'approche de la fente (i, «zone»). Pour les données planes, vous pouvez utiliser rgeos :: gArea (polys, byid = TRUE) Pour les données sphériques (lon / lat), vous pouvez utiliser geosphere :: areaPolygon

Je ne sais pas si vous voulez que le rapport soit basé sur la "valeur réelle" des zones de polygone (s) ou des zones des cellules qui les coupent. Voici un exemple de code qui utilise toutes les cellules coupant les polygones (essentiellement, le rapport des cellules NA aux cellules non NA). Il s'agit d'un exemple factice et vous devrez écrire votre propre fonction.

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

Vous pouvez utiliser cet extrait de code pour ajouter une colonne de zone à vos données de polygone en extrayant de l'emplacement de zone. Cela peut être utilisé si vous souhaitez effectuer un ratio à l'aide de la zone de polygone "réelle".

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

L'alternative la plus efficace et considérablement plus rapide aux boucles for sont des fonctions de type "appliquer". Il existe un certain nombre de ceux-ci disponibles dans R qui sont utilisés pour différentes classes d'objets ou structures de données. Dans ce cas, puisque extract renvoie une liste, nous utiliserons lapply (list apply). Il s'agit d'un moyen d'appliquer une fonction de base ou personnalisée à un objet de liste. La classe d'objets stockée dans la liste est un vecteur, la fonction est assez simple. Si vous utilisez l'extraction sur un objet raster de brique ou de pile, les objets résultants stockés dans la liste sont des objets 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)         

Je n'ai pas accès à vos fichiers, mais d'après ce que vous avez décrit, cela devrait fonctionner:

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)