extract () данные из растра с маленькими полигонами - округленные веса слишком малы

Используя R, я пытаюсь извлечь данные из растрового слоя, используя полигональный слой. Полигоны намного меньше ячеек растра:

введите описание изображения здесь

Теперь вызываю extract() из raster библиотеки:

a <- extract(raster, polygons, weights = TRUE, small = TRUE)
a
# ...
# [[1551]]
# value weight
#   209   0.03 # top left cell - more than 50% of the polygon area

Есть две проблемы: вес - это пропорция площади ячейки, покрытой многоугольником, и веса округляются до 1/100. В моем случае в выходных данных присутствует только верхняя левая ячейка (значение 209) - вес трех других ячеек был округлен до нуля, и они были исключены. Однако нижняя левая ячейка покрывает значительную часть многоугольника и тоже должна быть включена!

Мне нужно правильное средневзвешенное значение. Можно ли это как-то еще сделать с помощью extract()? Или любым другим способом?

PS: примечание в сторону: я думаю, что веса в extract() не очень хорошо спроектированы - вес должен быть пропорцией площади многоугольника, покрытой конкретной ячейкой, а не наоборот. Тогда будет легче вычислить взвешенное среднее для многоугольника (просто умножьте два числа в каждой строке и просуммируйте), и округление до 1/100 не будет большой проблемой.

Воспроизводимый пример - (скачать файлы - упрощенная версия, фактических данных намного больше):

require(raster)
rast <- raster("my.tif")
poly <- readOGR(".", "socc_buff_Rx")
a <- extract(rast, poly, weights = TRUE, small = TRUE)
a

По теме: Извлечь в R не удается для небольших многоугольников и растров


r gis raster r-raster spatial
person Tomas    schedule 20.07.2013    source источник
comment
Разве вы не можете расширить растровое разрешение, чтобы захватить (или приблизить) пропорцию покрытой области?   -  person Paulo E. Cardoso    schedule 21.07.2013
comment
@PauloCardoso, спасибо - на самом деле это была одна из моих 4 идей, как решить проблему: 1) убедить экстракт в разумном поведении, 2) преобразовать растровые ячейки в полигоны, 3) прочитать растр в R ячейке за ячейкой, фактически реализовав мой собственный extract, 4) увеличить разрешение растра за счет дальнейшего разбиения ячеек (ваша идея). Так что это зависит от того, какое решение на самом деле является самым простым, а также более эффективным с точки зрения вычислений. И я пока не знаю, как реализовать эти решения ...   -  person Tomas    schedule 21.07.2013
comment
Воспроизводимый пример сделает этот вопрос полезным.   -  person mdsumner    schedule 21.07.2013
comment
@mdsumner конечно. См. Обновленный вопрос.   -  person Tomas    schedule 21.07.2013
comment
@Tomas: У меня такой же вопрос, вы нашли ответ?   -  person Herman Toothrot    schedule 04.10.2013
comment
Я поддерживаю необходимость ответа на этот вопрос. Есть ли у @mdsumner какие-нибудь предложения?   -  person ego_    schedule 19.11.2013
comment
Я не нашел никакого решения для этого (я бы опубликовал его, если бы нашел). Мое обходное решение заключалось в увеличении полигонов (буферов). Реальным решением, вероятно, было бы отправить запрос на изменение авторам функции extract со ссылкой на эту страницу. Любой желающий может это сделать.   -  person Tomas    schedule 19.11.2013
comment
Вы предлагаете какой-нибудь систематический способ сделать это? Связано ли это с размером многоугольника или разрешением многоугольников? Я работаю над растром с разрешением 8х8 км.   -  person ego_    schedule 19.11.2013

Ответы (1)


arrow_upward
4
arrow_downward

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

library(raster)
r <- raster("my.tif")
pu <- shapefile("socc_buff_Rx.shp")
p <- spTransform(pu, crs(r))

extract(r, p, weights = TRUE, small = TRUE)
#[[1]]
# value weight
#   209   0.03

rr <- disaggregate(r, 10)
e <- extract(rr, p, weights = TRUE, small = TRUE)
lapply(e, function(x) { aggregate(x[,2,drop=F], list(value=x[,1]), sum ) } )

#[[1]]
#  value weight
#1   197   0.95
#2   209   3.44
#3   256   0.31
#4   293   0.04

plot(r, legend=F)
plot(p, add=T)
text(r)

растровые ячейки со значениями

person Robert Hijmans    schedule 22.12.2013