Вычисление площади от Больших пространственныхПикселейДанныхФрейма

Question

У меня есть объект с именем 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])
})

Но я в основном всегда получаю одни и те же предупреждающие сообщения

Предупреждающее сообщение: В .local (x, ...): Эта функция полезна только для растровых * объектов с координатами долготы / широты

Обратите внимание, что рассматриваемая площадь довольно мала (25 га).

Кто-нибудь может подтолкнуть меня в правильном направлении?

Answer 1

Вы должны предоставить простые данные, сгенерированные кодом, а не файл. Например

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

Answer 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