Если растровое значение NA, поиск и извлечение ближайшего пикселя, отличного от NA

При извлечении значений растра в точки я обнаружил, что у меня есть несколько NA, и вместо того, чтобы использовать аргументы buffer и fun функции extract, вместо этого я хотел бы извлечь ближайший не-NA пиксель к точке, которая перекрывает NA.

Я использую базовую функцию извлечения:

data.extr<-extract(loc.thr, data[,11:10])

person Joe    schedule 19.12.2014    source источник


Ответы (3)


Вот решение без использования буфера. Однако он вычисляет карту расстояний отдельно для каждой точки в вашем наборе данных, поэтому он может быть неэффективным, если ваш набор данных большой.

set.seed(2)

# create a 10x10 raster
r <- raster(ncol=10,nrow=10, xmn=0, xmx=10, ymn=0,ymx=10)
r[] <- 1:10
r[sample(1:ncell(r), size = 25)] <- NA

# plot the raster
plot(r, axes=F, box=F)
segments(x0 = 0, y0 = 0:10, x1 = 10, y1 = 0:10, lty=2)
segments(y0 = 0, x0 = 0:10, y1 = 10, x1 = 0:10, lty=2)

# create sample points and add them to the plot
xy = data.frame(x=runif(10,1,10), y=runif(10,1,10))
points(xy, pch=3)
text(x = xy$x, y = xy$y, labels = as.character(1:nrow(xy)), pos=4, cex=0.7, xpd=NA)

# use normal extract function to show that NAs are extracted for some points
extracted = extract(x = r, y = xy)

# then take the raster value with lowest distance to point AND non-NA value in the raster
sampled = apply(X = xy, MARGIN = 1, FUN = function(xy) r@data@values[which.min(replace(distanceFromPoints(r, xy), is.na(r), NA))])

# show output of both procedures
print(data.frame(xy, extracted, sampled))

#          x        y extracted sampled
#1  5.398959 6.644767         6       6
#2  2.343222 8.599861        NA       3
#3  4.213563 3.563835         5       5
#4  9.663796 7.005031        10      10
#5  2.191348 2.354228        NA       2
#6  1.093731 9.835551         2       2
#7  2.481780 3.673097         3       3
#8  8.291729 2.035757         9       9
#9  8.819749 2.468808         9       9
#10 5.628536 9.496376         6       6
person koekenbakker    schedule 19.12.2014
comment
Ваше решение успешно справилось с задачей с моим растром, который был не очень большим. - person Joe; 19.12.2014
comment
Спасибо koekenbakker, очень удобно. - person J. Win.; 17.08.2015
comment
Хм - но если вы не в памяти .... тогда это больше не работает. Мысли? - person jebyrnes; 16.09.2016
comment
Хорошее решение, но оно заставляет R зависать на большом количестве точек, 30 тыс. - person Herman Toothrot; 08.01.2018

Это растровое решение, при котором пиксели NA сначала заполняются ближайшим значением пикселя, не относящегося к NA. Однако обратите внимание, что при этом не учитывается положение точки в пикселе. Вместо этого он вычисляет расстояния между центрами пикселей, чтобы определить ближайший пиксель, не относящийся к числовой апертуре.

Во-первых, он вычисляет для каждого растрового пикселя NA расстояние и направление до ближайшего пикселя, не относящегося к NA. Следующим шагом является вычисление координат этой ячейки, не относящейся к Северной Америке (предполагается проекция CRS), извлечение ее значения и сохранение этого значения в местоположении Северной Северной Америки.

Исходные данные: спроецированный растр с идентичными значениями, как в ответе от koekenbakker:

set.seed(2)
# set projected CRS
r <- raster(ncol=10,nrow=10, xmn=0, xmx=10, ymn=0,ymx=10, crs='+proj=utm +zone=1') 
r[] <- 1:10
r[sample(1:ncell(r), size = 25)] <- NA

# create sample points
xy = data.frame(x=runif(10,1,10), y=runif(10,1,10))

# use normal extract function to show that NAs are extracted for some points
extracted <- raster::extract(x = r, y = xy)

Вычислите расстояние и направление от всех пикселей NA до ближайшего пикселя, не относящегося к NA:

dist <- distance(r)  
# you can also set a maximum distance: dist[dist > maxdist] <- NA
direct <- direction(r, from=FALSE)

Получить координаты пикселей NA

# NA raster
rna <- is.na(r) # returns NA raster

# store coordinates in new raster: https://stackoverflow.com/a/35592230/3752258 
na.x <- init(rna, 'x')
na.y <- init(rna, 'y')

# calculate coordinates of the nearest Non-NA pixel
# assume that we have a orthogonal, projected CRS, so we can use (Pythagorean) calculations
co.x <- na.x + dist * sin(direct)
co.y <- na.y + dist * cos(direct)

# matrix with point coordinates of nearest non-NA pixel
co <- cbind(co.x[], co.y[]) 

Извлечь значения ближайшей не-NA ячейки с координатами 'co'

# extract values of nearest non-NA cell with coordinates co
NAVals <- raster::extract(r, co, method='simple') 
r.NAVals <- rna # initiate new raster
r.NAVals[] <- NAVals # store values in raster

Заполнить исходный растр новыми значениями

# cover nearest non-NA value at NA locations of original raster
r.filled <- cover(x=r, y= r.NAVals)

sampled <- raster::extract(x = r.filled, y = xy)

# compare old and new values
print(data.frame(xy, extracted, sampled))

#           x        y extracted sampled
# 1  5.398959 6.644767         6       6
# 2  2.343222 8.599861        NA       3
# 3  4.213563 3.563835         5       5
# 4  9.663796 7.005031        10      10  
# 5  2.191348 2.354228        NA       3
# 6  1.093731 9.835551         2       2
# 7  2.481780 3.673097         3       3
# 8  8.291729 2.035757         9       9
# 9  8.819749 2.468808         9       9 
# 10 5.628536 9.496376         6       6

Обратите внимание, что точка 5 принимает другое значение, чем ответ Koekenbakker, поскольку этот метод не учитывает положение точки в пикселе (как упоминалось выше). Если это важно, это решение может оказаться неподходящим. В других случаях, например. если растровые ячейки малы по сравнению с точечной точностью, этот растровый метод должен дать хорошие результаты.

person Lennert    schedule 23.05.2016
comment
Отличное объяснение, я рад, что вы разъяснили это - person Joe; 26.05.2016
comment
возможно, это может квалифицироваться как другой вопрос ... но было бы неплохо, чтобы новая точка xy с не-NA возвращалась вместе с соответствующими данными - person Joe; 01.11.2016

Для растрового стека используйте решение @koekenbakker выше и превратите его в функцию. Слот @layers растрового стека представляет собой список растров, поэтому наложите его внахлест и продолжайте оттуда.

#new layer
r2 <- raster(ncol=10,nrow=10, xmn=0, xmx=10, ymn=0,ymx=10)
r2[] <- 1:10
r2[sample(1:ncell(r2), size = 25)] <- NA

#make the stack
r_stack <- stack(r, r2)

#a function for sampling
sample_raster_NA <- function(r, xy){
  apply(X = xy, MARGIN = 1, 
        FUN = function(xy) r@data@values[which.min(replace(distanceFromPoints(r, xy), is.na(r), NA))])

}

#lapply to get answers
lapply(r_stack@layers, function(a_layer) sample_raster_NA(a_layer, xy))

Или, чтобы быть модным (улучшение скорости?)

purrr::map(r_stack@layers, sample_raster_NA, xy=xy)

Что заставляет меня задаться вопросом, можно ли еще больше ускорить все это с помощью dplyr...

person jebyrnes    schedule 16.09.2016