从 Large spatialPixelsDataFrame 计算面积
Calculating area from Large spatialPixelsDataFrame
我有一个名为 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 公顷)。
谁能把我推向正确的方向?
我真的很困惑。您说您正在处理 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
您应该提供简单的代码生成示例数据,而不是文件。例如
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")
我有一个名为 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 公顷)。
谁能把我推向正确的方向?
我真的很困惑。您说您正在处理 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
您应该提供简单的代码生成示例数据,而不是文件。例如
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")