CODEWITHSUNDEEP
rrasterr-raster
Leprechault
Leprechault

Reputation: 1823

Computing number of points in a raster stack in R

I'd like to compute the number of points in a raster stack, including zeros. For this I make:

#Packages
library(raster)
library(sp)
library(rgdal)

## Create a simulated RBG rasters
r <- raster(nc=30, nr=30)
r <- setValues(r, round(runif(ncell(r))* 255))
g <- raster(nc=30, nr=30)
g <- setValues(r, round(runif(ncell(r))* 255))
b <- raster(nc=30, nr=30)
b <- setValues(r, round(runif(ncell(r))* 255))
rgb<-stack(r,g,b)
plotRGB(rgb,
        r = 1, g = 2, b = 3)

##Given interesting points coordinates
xd     <- c(-24.99270,45.12069,99.40321,73.64419)
yd  <- c(-45.435267,-88.369745,-7.086949,44.174530)
pts <- data.frame(xd,yd)
pts_s<- SpatialPoints(pts)
points(pts_s, col="black", pch=16)

#Count number of points inside each raster
res<-NULL
for(i in 1:3){
  pointcount = function(r, pts){
  # make a raster of zeroes like the input
  r2 = r
  r2[] = 0
  # get the cell index for each point and make a table:
  counts = table(cellFromXY(r,pts))
  # fill in the raster with the counts from the cell index:
  r2[as.numeric(names(counts))] = counts
  return(r2)
 }
r2 = pointcount(rgb[[i]], pts_s)
res<-rbind(res,c(r2))
}
#

> res
     [,1]
[1,] ?   
[2,] ?   
[3,] ?

And my function doesn't work. I expected 4, 4 and 4 in res output. What is my error?

Upvotes: 0

Views: 484

Answers (2)

Robert Hijmans
Robert Hijmans

Reputation: 47146

It is not clear what you mean with counting "points in a raster". But below are some suggestions. This is R -- when you start writing a loop for this type of question you are almost certainly overlooking a more direct approach

Example data

r <- raster(nc=30, nr=30)
values(r) <- 1:ncell(r)
s <- stack(r, r*2)
# 5 points, 1 outside raster
xd  <- c(-24,45,99,73,200)
yd  <- c(-45,-88,-7,44,100)
pts <- cbind(xd,yd)

Do you want to know if the points intersect with the raster? Then you can do

 sum(!is.na(cellFromXY(r, pts)))
 #[1] 4

Or

 length(intersect(SpatialPoints(pts), r))
 #[1] 4

I suppose you do not want this, as the answer is the same for all layers in the RasterStack (and you can use either r or s)

Or do you want to know which points are at cells that are not NA?

e <- extract(s, pts)
colSums(!is.na(e))
#layer.1 layer.2 
#      4       4

Or count the number of points per cell?

rp <- rasterize(pts, r, fun="count")
freq(rp)
#     value count
#[1,]     1     4
#[2,]    NA   896

Upvotes: 1

ThetaFC
ThetaFC

Reputation: 670

Added a few fixes to your loop: (1) pre-allocated res as an integer vector of length 3; (2) made pointcount return the sum of counts, which is what you say you want; (3) filled res with the sum result of each iteration.

res<-rep(NA_integer_, 3) #pre-allocates the result vector that will be filled in the loop
for(i in 1:3){
  pointcount = function(r, pts){
    # make a raster of zeroes like the input
    #r2 = r #not necessary
    #r2[] = 0 #not necessary
    # get the cell index for each point and make a table:
    counts = table(cellFromXY(r,pts))
    # fill in the raster with the counts from the cell index
    #r2[as.numeric(names(counts))] = counts #don't need this
    #return(r2) #don't need this
    sum(counts) #return the sum of counts; 'return' is implied so not necessary code
  }
  result = pointcount(rgb[[i]], pts_s)
  #res<-rbind(res,c(r2)) #not functional: r2 is a raster, so rbind does not make sense
  res[i] = result #fill the correct slot of res with result (i.e. sum of counts)
}

Upvotes: 0

Related Questions