Проблема с оценкой взвешенного среднего по растру для формы многоугольника в R - PullRequest
Новая
       8

Проблема с оценкой взвешенного среднего по растру для формы многоугольника в R

0 голосов
/ 19 декабря 2018

Мне нужно оценить средневзвешенное значение растровых значений для многоугольника, показанного в квадратах.Я хочу получить значение растра и его вес в каждом квадрате в форме многоугольника.(Как показано в этом сообщении: Как я могу извлечь взвешенную сумму площади из растра в многоугольник в R? )

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

Пожалуйста, смотрите прикрепленные входные данные, установленные здесь: https://bft.usu.edu/w8crs

Спасибо. 

library(raster)
library(sp)
library(rgdal)
library(rgeos)

rlist = list.files(getwd(), pattern = "tif$", full.names = TRUE)
inshp = "Test" 
rdata <- rlist[1]

r <- raster(rdata)
sdata <- readOGR(dsn=getwd(), layer=inshp)
sdata <- spTransform(sdata, crs(r))

extract(r, sdata, weights=TRUE)

enter image description here

Выход:

enter image description here

[[1]]
    value weight
 56.75139      1

[[2]]
    value weight
 61.18781      1

[[3]]
    value weight
 56.75139      1

[[4]]
    value weight
 61.18781      1

1 Ответ

0 голосов
/ 19 декабря 2018

Вот воспроизводимый пример 

library(raster)
packageVersion("raster")
#[1] ‘2.8.4’
r <- raster(xmn=0, xmx=1, ymn=0, ymx=1, nrow=2, ncol=2)
values(r) <- 1:4

m <- matrix(c(0.4, 0.6, 0.8, 0.6, 0.7, 0.2, 0.3, 0.2), ncol=2, byrow=TRUE)
s <- spPolygons(m)
plot(r); lines(s)

extract(r, s, weights=TRUE)
#[[1]]
#     value weight
#[1,]     1 0.0625
#[2,]     2 0.1875
#[3,]     3 0.3125
#[4,]     4 0.4375

Это не сработало для вас, потому что ваш полигон был очень мал по сравнению с размером растровой ячейки.Я изменил функцию, так что она повышает точность для этих случаев.Теперь я получаю это с вашими данными: 

> extract(r, sdata, weights=TRUE)
[[1]]
    value weight
 56.75139      1

[[2]]
        value    weight
[1,] 61.18781 0.6592593
[2,] 56.75139 0.3407407

[[3]]
    value weight
 56.75139      1

[[4]]
        value    weight
[1,] 61.18781 0.5522388
[2,] 56.75139 0.4477612

Чтобы сделать его воспроизводимым без загрузок, для одного из ваших полигонов: 

library(raster)
r <- raster(ncol=2, nrow=1, xmn=596959.624056728, xmx=624633.120455544, ymn=568805.230192675, ymx=582641.978392083, crs='+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m')
values(r) <- c(61.18781, 56.75139)    
g <- data.frame(matrix(c(rep(1, 18), rep(0,6), 611318.079488842,611440.751254539,610712.115334383,609842.749239201, 609703.303842618,611318.079488842,581038.816616668,579434.971927127, 579381.167042005,579315.223934334,580917.724282178,581038.816616668), ncol=6))
colnames(g) <- c('object','part','cump','hole','x','y')
p <- as(g, "SpatialPolygons")
crs(p) <- crs(r)

extract(r, p, weights=TRUE)

#[[1]]
#        value    weight
#[1,] 61.18781 0.6592593
#[2,] 56.75139 0.3407407
...