R растр :: функция извлечения не работает должным образом

Question

Я пытаюсь сделать несколько полигонов / растровых работ в R. Функция извлечения не работает, как я ожидал.Вот воспроизводимый пример.

library(sf)
library(raster)
library(maptools)

raster      <-     download.file('https://raw.githubusercontent.com/JimShady/london_osm_canyons/master/test.tif', 'downloaded_raster')
raster      <- raster('test.tif')

polygons    <- read_sf('https://raw.githubusercontent.com/JimShady/london_osm_canyons/master/test2.geojson')
polygons    <- st_transform(polygons, crs(raster)@projargs)
polygons$id <- 1
polygons    <- as(polygons, 'Spatial')

new_polygon <- elide(polygons, rotate=-50)
proj4string(new_polygon) <- crs(raster)@projargs

polygons    <- rbind(polygons, new_polygon)

plot(raster)
plot(polygons, add=T)

Я хочу получить количество ячеек, с которыми пересекается каждый многоугольник, даже если они NA,и сохранить их в многоугольнике в виде столбца

polygons$cell_count     <- extract(raster, polygons, fun = function(x,...)length(x), na.rm=F)

Но это возвращает число 5 для многоугольника один, когда ясно, что многоугольник пересекает более 5 ячеек растра. Что с этим происходит?

Затем я хочу получить средневзвешенное значение.если есть какие-либо пересеченные ячейки NA, то это должно «потерпеть неудачу» и вернуть NA. Я думаю, что это работает, но поскольку вышеприведенное неверно, я не очень уверен в своем методе сейчас.

polygons$mean           <- round(extract(raster, polygons, fun = mean, weights=T, na.rm=F),2)

Я также хочу знать количество ячееккоторые пересекаются, которые меньше единицы. Я не уверен, как это сделать.У кого-нибудь есть хорошая идея?

polygons$almost_zero    <-

Спасибо за любую помощь.

Answer 1

Извлечение ведет себя так, как и предполагалось, находя ячейки, в которых центроид находится внутри многоугольника.Один из способов справиться с частичным покрытием - преобразовать полигон в растровую маску, используя rasterize.Обратите внимание, что в дальнейшем я переименовал ваш растр и полигон r1 и p1 соответственно.Поскольку не рекомендуется использовать имена функций для объектов:

r_mask = rasterize(p1, r1, getCover=TRUE)
r_mask[r_mask==0] = NA

Теперь мы можем использовать эту маску, чтобы получить нужные значения:

r2 = mask(r1, r_mask)
cellStats(!is.na(r2), sum)
# [1] 18
cellStats(r2, mean)
# [1] 1.944444

Чтобы определить, сколько ячеек меньшечем один, мы можем сделать

cellStats(r2<1, sum)
# [1] 13

Answer 2

Из extract help:

Ячейка покрыта, если ее центр находится внутри многоугольника (но см. Параметр весов для рассмотрения частично покрытых ячеек и аргумент small для получения значений для маленьких многоугольниковв любом случае).

Следовательно, вы можете использовать это, чтобы иметь количество ячеек:

data <- extract(raster, polygons,  weights = T,  na.rm = F)[[1]]
nrow(data)

#[1] 18

Где data:

      value      weight
 [1,]     8 0.003610108
 [2,]     6 0.025270758
 [3,]     0 0.010830325
 [4,]     7 0.104693141
 [5,]     0 0.066787004
 [6,]     0 0.088447653
 [7,]     0 0.001805054
 [8,]     0 0.104693141
 [9,]     0 0.010830325
[10,]     0 0.093862816
[11,]     0 0.030685921
[12,]     0 0.063176895
[13,]     0 0.057761733
[14,]     0 0.063176895
[15,]     7 0.003610108
[16,]     7 0.102888087
[17,]     0 0.093862816
[18,]     0 0.074007220

Комуиметь количество ячеек, которые пересекаются и меньше чем 1:

sum(data[, 1] < 1)

#[1] 13

---

РЕДАКТИРОВАТЬ

В случае нескольких полигонов, вы должны сначала растворить ваши полигоныиспользуя rgeos::gUnaryUnion:

library(rgeos)

polygons$id <- 1
polygons <- gUnaryUnion(polygons, id = polygons$id)
# As suggested by @RobertHijmans, you can also use `raster::aggregate`
# polygons <- raster::aggregate(polygons, by = "id")

plot(raster)
plot(polygons, add = T)