Растр в R: создание зонального подсчета определенных значений ячеек без переклассификации

Я хотел бы знать, есть ли способ создать zonal statistics для RasterLayerObjects, в частности count заданного значения ячейки (например, класса землепользования) в R без необходимости переклассифицировать весь растр. Решение должно эффективно использовать память для работы с большими растровыми файлами, т. е. извлечение значений в матрицу в R нежелательно.

Ниже пример того, как я справляюсь с этим до сих пор. В этом случае я реклассифицирую исходный растр, чтобы он содержал только 1 для интересующего значения и отсутствовал для всех остальных значений.

Мое предлагаемое решение создает как избыточные данные, так и дополнительные этапы обработки, чтобы приблизить меня к моей первоначальной цели. Я думал, что что-то вроде zonal(r1[r1==6],r2,"count") будет работать, но, очевидно, это не так (см. ниже).

# generate reproducible Raster 
library("raster")

## RASTER 1 (e.g. land-use classes)
r1 <- raster( crs="+proj=utm +zone=31")
extent(r1) <- extent(0, 100, 0, 100)
res(r1) <- c(5, 5)
values(r1) <- sample(10, ncell(r1), replace=TRUE)
plot(r1)

## RASTER 2 (containing zones of interest)
r2 <- raster( crs="+proj=utm +zone=31")
extent(r2) <- extent(0, 100, 0, 100)
res(r2) <- c(5, 5)
values(r2) <- c(rep(1,100),rep(2,100),rep(3,100),rep(4,100))
plot(r2)

# (1) ZONAL STATISTICS
# a. how many cells per zone (independent of specific cell value)
zonal(r1,r2,"count")

# b. how many cells per zone of specific value 6
zonal(r1[r1==6],r2,"count")
# -> fails

# with reclassification
r1.reclass<-
  reclassify(r1,
             matrix(c(1,5,NA,
                    5.5,6.5,1, #class of interest 
                    6.5,10,NA),
                    ncol=3,
                    byrow = T),
             include.lowest=T # include the lowest value from the table.
           )
zonal(r1.reclass,r2,"count")

r geospatial raster
person joaoal    schedule 07.12.2017    source источник

Ответы (1)


arrow_upward
2
arrow_downward

вы можете использовать raster::match.

zonal(match(r1, 6),r2, "count")

Как видно из plot(match(r1, 6)), он возвращает только те растровые ячейки, которые содержат нужные значения. Все остальные ячейки NA.

r1==6, использованный в вашей попытке, к сожалению, возвращает вектор и, следовательно, больше не может использоваться в фокусе.

person loki    schedule 07.12.2017
comment
именно то, что мне нужно. Спасибо!! - person joaoal; 07.12.2017
comment
Я должен упомянуть, что растровый пакет создаст временные файлы для выполнения этого расчета, и они огромны. если место на диске заполнено, используйте options(rasterTmpDir='/some_external_drive'), чтобы временно поместить файлы на внешний диск. - person joaoal; 07.12.2017
comment
Я также должен упомянуть, что этот подход на самом деле не быстрее, чем ручная реклассификация. Вероятно, потому что растровый пакет не очень эффективен при создании файлов tmp. Мой темп. растры имеют размер 25 ГБ по сравнению с моим вручную переклассифицированным бинарным растром, который потребляет всего 100 МБ... - person joaoal; 07.12.2017