Вычисление площади от Large пространственногоПиксели DataFrame
У меня есть объект с именем cr1, который является большой SpatialPixelsDataFrame озера. Вот ссылка на файл: https://www.dropbox.com/s/uuvlmxmri144hp2/macrosmall.rdata?dl=0
Я думаю, что каждый пиксель имеет размер ячейки 1 х 1 м, однако я думаю, что этот атрибут нигде не указан. "макро" - это измеренная высота погруженных макрофитов в озере. Структура выглядит следующим образом.
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
Я хотел бы вычислить область, которая покрыта для определенных интервалов высоты макрофита (то есть область, покрытая для интервалов "макроса").
Как я могу указать разрешение или размер каждой ячейки (=1 м х 1 м)? Какой пакет и функция обрабатывает оценку области 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])
})
Но я в основном всегда получаю одни и те же предупреждающие сообщения
Предупреждающее сообщение: In .local(x, ...): эта функция полезна только для растровых * объектов с координатами долготы / широты
Обратите внимание, что рассматриваемая площадь довольно мала (25 га).
Кто-нибудь может подтолкнуть меня в правильном направлении?
2 ответа
Я действительно смущен. Вы сказали, что работаете над 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 м ^2, общая площадь озера составляет 36143
м ^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
Решение состоит в том, чтобы сначала сделать интересующие вас занятия. Вы можете использовать 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)
А затем посчитайте количество клеток в каждом классе:
f <- freq(rc)
Пока ваш CRS не является долготой / широтой, вы можете умножить количество ячеек на площадь ячейки, чтобы получить область (но в вашем случае эта область равна 1, поэтому в этом нет необходимости).
f <- data.frame(f)
f$area <- f$count * prod(res(rc))
Если данные на lon/lat, вам нужно сделать
a <- area(rc)
ff <- zonal(a, rc, "sum")