从大空间像素数据帧计算面积

Question

我有一个名为 cr1 的对象，它是湖的大型 SpatialPixelsDataFrame。 这是该文件的链接： https://www.dropbox.com/s/uuvlmxmri144hp2/macrosmall.rdata?dl=0

我认为每个像素都有一个 1m x 1m 的单元格大小，但是我认为没有指定此属性。 “宏”是湖中沉水植物的测量高度。 结构是这样的。

    Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots

  ..@ data       :'data.frame': 252234 obs. of  1 variable:
  .. ..$ macro: num [1:252234] 0.0468 0.0518 0.0445 0.046 0.0477 ...

  ..@ coords.nrs : num(0) 

  ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots

  .. .. ..@ cellcentre.offset: Named num [1:2] 3404494 5872334

  .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
  .. .. ..@ cellsize         : Named num [1:2] 1 1

  .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
  .. .. ..@ cells.dim        : Named int [1:2] 776 536

  .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
  ..@ grid.index : int [1:252234] 415333 415334 415335 415336 415337 415338 
415339 414554 414555 414556 ...

  ..@ coords     : num [1:252234, 1:2] 3404666 3404667 3404668 3404669 3404670 ...
  .. ..- attr(*, "dimnames")=List of 2

  .. .. ..$ : chr [1:252234] "949" "950" "951" "952" ...

  .. .. ..$ : chr [1:2] "x" "y"

  ..@ bbox       : num [1:2, 1:2] 3404493 5872333 3405269 5872869

  .. ..- attr(*, "dimnames")=List of 2

  .. .. ..$ : chr [1:2] "x" "y"

  .. .. ..$ : chr [1:2] "min" "max"

  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

  .. .. ..@ projargs: chr NA

我想计算某些大型植物高度间隔所覆盖的面积（即“宏观”间隔所覆盖的面积）。

如何指定每个单元格的分辨率或大小 (=1m x 1m)？ 哪个包和函数处理 SpatialPixelsDataFrame 的面积估计？

其实我目前只加载了地图

library(sp)
library(raster)

load("macrosmall.rdata")

并尝试了几件事：

area(cr1)

这将是我想要什么以及我想如何计算它的一个例子，但是数据框的规范不允许它

intervals <- list(c(0.1,0.2), 
              c(0.2,0.3),
              c(0.3,0.4))

sapply(intervals, function(x) { 
  sum(cr1[] > x[1] & cr1s[] <= x[2])
})

但我基本上总是以同样的警告信息告终

警告信息： 在 .local(x, ...) 中： 此功能仅对具有经度/纬度坐标的 Raster* 对象有用

请注意，该区域非常小（25 公顷）。

谁能把我推向正确的方向？

Answer 1

我真的很困惑。您说您正在处理spatialPixelsDataFrame，但您提供的数据cr1 是一个栅格对象。

无论如何，这里我展示了一个计算面积的可能解决方案。

library(sp)
library(raster)

# Load the RData
load("macrosmall.RData")

# View the raster layer
cr1
# class       : RasterLayer 
# dimensions  : 200, 269, 53800  (nrow, ncol, ncell)
# resolution  : 1, 1  (x, y)
# extent      : 3405000, 3405269, 5872500, 5872700  (xmin, xmax, ymin, ymax)
# coord. ref. : NA 
# data source : in memory
# names       : macro 
# values      : 0, 1.896009  (min, max)

# Plot the raster layer
plot(cr1)

我不确定哪些栅格值表示“湖”。假设所有非 NA 单元格都是湖，那么我们可以执行以下操作。

# Get all cell values
cells <- values(cr1)

# Remove NA
cells_nonNA <- cells[!is.na(cells)]

# Count how many cells
length(cells_nonNA)
# [1] 36143

由于 1 个单元格是一个 1 m^2，因此总湖面积为36143 m^2。

假设湖泊的栅格值大于 1，我们可以再次对像元值进行子集化。

cells_above1 <- cells_nonNA[cells_nonNA > 1]
length(cells_above1)
# [1] 77

Answer 2

您应该提供简单的代码生成的示例数据，而不是文件。例如

library(raster)
r <- raster(nrow=20, ncol=26, ext=extent(3405000, 3405269, 5872500, 5872700))
values(r) <- 1:ncell(r) / 280
set.seed(-1)
r[sample(ncell(r), 100)] <- 0

解决方法是先制作你感兴趣的类。你可以使用cut 或reclassify。这里有reclassify：

m <- matrix(c(0, 0.1, 1,
              0.1, 0.5, 2,
              0.5, 1, 3,
              1, 2, 4), ncol=3, byrow=TRUE)

rc <- reclassify(r, m)

然后统计每个类的cell个数：

f <- freq(rc)

只要您的 CRS 不是经度/纬度，您就可以将单元格计数乘以单元格面积来获得面积（但在您的情况下面积为 1，因此没有必要）。

f <- data.frame(f)
f$area <- f$count * prod(res(rc))

如果数据在经度/纬度上，您需要这样做

a <- area(rc)
ff <- zonal(a, rc, "sum")