Calculating number of pixel(count) raster files in R

whj picture whj · Jun 3, 2018 · Viewed 9.9k times · Source

I have a huge number of raster files and a polygon that is within the extent of the raster files. I want get the pixel number(count) for each raster files within the polygon. Additionally I want to create a table with the name of raster file and number of pixels(listed) for each raster file. I have tried with stacking but with that I can not keep track of the name. Is there any other way of performing this task in R?

Answer

Robert Hijmans picture Robert Hijmans · Jun 4, 2018

Always include example data, please

library(raster)
fn <-system.file("external/rlogo.grd", package="raster")
s <- stack(fn)
s[[1]][1:5000] <- NA
s[[2]][5001:ncell(s)] <- NA
names(s)
#[1] "red"   "green" "blue" 
p <- rbind(c(5,20), c(25,55), c(50, 20), c(20,6), c(5,20))
pol <- spPolygons(p)

plot(s, addfun=function() lines(pol, lwd=2))

I am not quite sure what you are after. The number of cells (pixels) would be the same for all rasters if you can stack them (which you say you can). I am assuming that you want the sum of the cells that are not NA. If you actually have rasters with a different origin/resolution, you can repeat these steps, but there no need to stack them into a RasterStack, but you would need to adjust the approach to also count the NA cells.

Simple approach for smaller objects:

m <- mask(s, pol)
cellStats(m, function(i, ...) sum(!is.na(i)))
# red green  blue 
# 600   506  1106 

If that runs out of memory, you can do:

m <- mask(s, pol)
x <- reclassify(m, cbind(-Inf, Inf, 1))
names(x) <- names(m)
cellStats(x, 'sum')
#red green  blue 
#600   506  1106 

You can also try:

extract(s, pol, fun=function(x,...)length(na.omit(x)))
#     red green blue
#[1,] 600   506 1106

If you want to count all the cells (whether NA or not), you can do something like

# example RasterLayer
r <- s[[1]]
# this step may help in speed if your polygon is small relative to the raster
r <- crop(r, pol)

x <- rasterize(pol, r, 1)
cellStats(x, 'sum')
#[1] 1106