So ermitteln Sie mithilfe von R die Anzahl der Nicht-NA-Rasterzellen im Polygon

Ich bin mit ArcGIS ZonalStats auf alle möglichen Probleme gestoßen und dachte, R könnte eine großartige Möglichkeit sein. Ich bin ziemlich neu in R, habe aber einen Codierungshintergrund.

Die Situation ist, dass ich mehrere Raster und ein Polygon-Shapefile mit vielen Features unterschiedlicher Größe habe (obwohl alle Features größer als eine Rasterzelle sind und die Polygon-Features auf das Raster ausgerichtet sind). Ich habe herausgefunden, wie der Mittelwert für jedes Polygon-Feature mithilfe der Raster-Bibliothek mit Extrakt ermittelt werden kann:

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

Das Problem, das ich habe, ist, dass ich auch ein flächenbasiertes Verhältnis zwischen der Fläche des Polygons und allen Nicht-NA-Zellen im selben Polygon benötige. Ich weiß, wie groß die Zellen des Rasters sind, und ich kann die Fläche für jedes Polygon ermitteln, aber das fehlende Glied ist die Anzahl aller Nicht-NA-Zellen in jedem Feature. Ich habe es geschafft, die Zellennummer aller Zellen im Polygon zu ermitteln, proxNA@data$Cnumb1000 <- cellFromPolygon(ras, proxNA)und ich bin sicher, dass es eine Möglichkeit gibt, den tatsächlichen Wert der Rasterzelle zu ermitteln, für die dann eine Schleife erforderlich ist, um die Anzahl aller Nicht-NA-Zellen zu ermitteln eine Zählung usw. ABER ich bin sicher, dass es einen viel besseren und schnelleren Weg gibt, das zu tun! Wenn jemand von euch eine Idee hat oder mich in die richtige Richtung weisen kann, wäre ich sehr dankbar!

Die Beispieldaten von 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"))))

Verwenden Sie nun Extrakt

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

Wenn Sie viele Raster haben, kombinieren Sie diese zuerst mit dem Stapel (wenn sie den gleichen Umfang und die gleiche Auflösung haben).

Um den tatsächlichen Polygonbereich zu erhalten, sollten Sie nicht den Slot-Ansatz (i, 'area') verwenden. Für planare Daten können Sie rgeos :: gArea (polys, byid = TRUE) verwenden. Für sphärische Daten (lon / lat) können Sie geosphere :: areaPolygon verwenden

Ich bin mir nicht sicher, ob Sie das Verhältnis basierend auf dem "realen Wert" der Polygonbereiche oder der Bereiche der Zellen, die sie schneiden, möchten. Hier ist ein Beispielcode, der alle Zellen verwendet, die die Polygone schneiden (im Grunde das Verhältnis von NA-Zellen zu Nicht-NA-Zellen). Es ist ein Dummy-Beispiel und Sie müssen Ihre eigene Funktion schreiben.

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

Mit diesem Codeausschnitt können Sie Ihren Polygondaten eine Bereichsspalte hinzufügen, indem Sie sie aus dem Bereichsbereich extrahieren. Dies kann verwendet werden, wenn Sie das Verhältnis unter Verwendung des "echten" Polygonbereichs verwenden möchten.

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

Die effizienteste und wesentlich schnellere Alternative zu for-Schleifen sind "Apply" -ähnliche Funktionen. In R sind einige davon verfügbar, die für verschiedene Objektklassen oder Datenstrukturen verwendet werden. In diesem Fall verwenden wir lapply (Liste gilt), da der Extrakt eine Liste zurückgibt. Auf diese Weise können Sie eine Basis- oder benutzerdefinierte Funktion auf ein Listenobjekt anwenden. Die in der Liste gespeicherte Objektklasse ist ein Vektor, die Funktion ist recht einfach. Wenn Sie Extrakt für ein Baustein- oder Stapelrasterobjekt verwenden, sind die in der Liste gespeicherten resultierenden Objekte data.frame-Objekte.

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

Ich habe keinen Zugriff auf Ihre Dateien, aber basierend auf dem, was Sie beschrieben haben, sollte dies funktionieren:

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)