在R中计算栅格堆栈中的点数答案

【问题标题】：Computing number of points in a raster stack in R在R中计算栅格堆栈中的点数
【发布时间】：2023-04-04 13:08:01
【问题描述】：

我想计算栅格堆栈中的点数，包括零。为此我做了：

#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,] ?

而且我的功能不起作用。我预计res 输出中有 4、4 和 4。我的错误是什么？

【问题讨论】：

标签： r raster r-raster

【解决方案1】：

不清楚计算“栅格中的点”是什么意思。但以下是一些建议。这就是 R——当你开始为这类问题编写循环时，你几乎肯定会忽略一种更直接的方法

示例数据

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)

您想知道这些点是否与栅格相交吗？然后就可以了

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

或者

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

我想你不想要这个，因为RasterStack 中的所有层的答案都是相同的（你可以使用r 或s）

或者您想知道哪些点位于非NA 的单元格处？

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

或者计算每个单元格的点数？

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

【讨论】：

【解决方案2】：

为您的循环添加了一些修复：(1) 预先分配 res 作为长度为 3 的整数向量；（2）让pointcount返回计数的总和，这就是你说的你想要的； (3) 用每次迭代的总和结果填充res。

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

【讨论】：