从 Large spatialPixelsDataFrame 计算面积

Question

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

我认为每个像素都有一个 1m x 1m 的像元大小，但是我认为这个属性没有指定。 "macro" 是湖中沉水水生植物的实测高度。结构如下所示。

    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

我想计算某些大型植物高度间隔覆盖的面积（即 "macro" 间隔覆盖的面积）。

如何指定每个单元格的分辨率或大小 (=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])
})

但我基本上总是收到相同的警告消息

Warning message: In .local(x, ...) : This function is only useful for Raster* objects with a longitude/latitude coordinates

请注意，相关区域很小（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)

我不确定哪些栅格值表示 "lake"。假设所有的非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

一个解决方案是先制作你感兴趣的classes。你可以使用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)

然后统计每个单元格的个数class:

f <- freq(rc)

只要你的CRS不是longitude/latitude，你可以用单元格数乘以单元格面积得到面积（但你的面积是1，所以没有必要）。

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

如果数据在 lon/lat 上，您需要执行

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

从 Large spatialPixelsDataFrame 计算面积

Calculating area from Large spatialPixelsDataFrame

maps

r

raster

area

spatial-data-frame