How can I get the number of pixels with NA value in a raster that is cliped from a large raster by many polygons?

0

Issue

I applied cloud mask to a raster image in R, and want to check how many pixels are masked out. But what I really need are only the images within some polygons (400+ of them), so I only want to get the number of pixels with no value within the polygons.

Here is what I have done:

library(raster)
library(rgdal)

##Read the raster files

tb = raster('D:/HLS/NDVI_Month_2018_TB.tif', band = 6)

##Read the polygon (400 polygons)

crops = readOGR('D:/HLS/shapefile/tb/tb.shp')

##reproject the vector

new_crops = spTransform(crops, crs(tb))

##Clip the raster with polygons

cliped = crop(tb, extent(new_crops))
output = mask(cliped, new_crops)

##Check the NA value

freq(output, value = NA)

However what I got from the freq() function seems to be all the pixels within the area (not only the polygons but the area from crop() function).

The result of freq():
1

How can I get the NA value within the polygons?

Solution

Here is a minimal, self-contained, reproducible example (taken mostly from ?raster::extract)

Example raster and polygons

library(raster)
r <- raster(ncol=90, nrow=45)
values(r) <- 1:ncell(r)
r[seq(1,ncell(r),3)] <- NA

p1 <- rbind(c(-180,-20), c(-140,55), c(0, 0), c(-140,-60), c(-180,-20))
p2 <- rbind(c(10,0), c(140,60), c(160,0), c(140,-55), c(10,0))

pols <- spPolygons(p1, p2)

Solution 1

extract(r, pols, fun=function(i, ...) sum(is.na(i)))
#     [,1]
#[1,]  215
#[2,]  178

Solution 2

z <- rasterize(pols, r)
zonal(is.na(r), z, "sum")
#     zone sum
#[1,]    1 215
#[2,]    2 178

Answered By – Robert Hijmans

This Answer collected from stackoverflow, is licensed under cc by-sa 2.5 , cc by-sa 3.0 and cc by-sa 4.0

Leave A Reply

Your email address will not be published.

This website uses cookies to improve your experience. We'll assume you're ok with this, but you can opt-out if you wish. Accept Read More