R terra получает соответствующие ячейки одного растра, которые перекрывают другой растр
У меня есть растр с высоким разрешением (highRes) и растр с более низким разрешением (lowRes), и я хотел бы вернуть для каждой ячейки с низким разрешением значения и процент покрытия ячеек с высоким разрешением в каждой ячейке с низким разрешением.
Я могу сделать это с помощью следующего кода, и он работает. Здесь я иду ячейка за ячейкой, преобразуя каждую ячейку в многоугольник, чтобы определить, какие ячейки высокого разрешения пересекаются этим многоугольником. Но это невероятно медленно с большими растрами. Для проекта, над которым я работаю, растр высокого разрешения имеет глобальное разрешение 1x1 км, а растр низкого разрешения — 25x25 км.
Мне интересно, есть ли более быстрый способ сделать это.
library(terra)
highRes <- rast()
res(highRes) <- c(0.1, 0.1)
highRes[] <- sample(1:10, ncell(highRes), replace = TRUE)
lowRes <- rast(res = c(1, 1), xmin = -170, xmax = -50, ymin = 0, ymax = 80)
lowRes[] <- 1
lowRes <- shift(lowRes, dx = 0.02, dy = 0.03) # not aligned
for (i in 1:ncell(lowRes)) {
xx <- rast(lowRes, vals = NA)
xx[i] <- 1
p <- as.polygons(xx)
p <- project(p, crs(highRes))
e <- extract(highRes, p, weights = TRUE, exact = TRUE)
}
Это дает мне значения ячеек, которые перекрываются с ячейкой lowRes, и процент покрытия этих ячеек.
Ответы 2
Создайте один полигон SpatVector и/или sf со всеми представленными ячейками, затем запустите extract() на crop()ed-версии highRes. Вы также можете использовать exactextractr::exact_extract(), хотя он не работает с объектами SpatVector и генерирует объект списка. Однако это намного быстрее, чем extract().
library(terra)
library(sf)
library(exactextractr)
# Example SpatRaster data
highRes <- rast(ext(-180, 180, -90, 90), res = 1000/111320, crs = "EPSG:4326")
set.seed(1)
highRes [] <- sample(1:10, ncell(highRes), replace = TRUE)
lowRes <- rast(ext(-170, -50, 0, 80), res = 25000/111320, crs = crs(highRes))
lowRes <- shift(lowRes, dx = 0.02, dy = 0.03)
# Add individual value to each cell
lowRes[] <- 1:ncell(lowRes)
# Extend lowRes for cropping highRes (to account for shift)
lowRes_c <- rast(extend(ext(lowRes), c(1,1)), crs = crs(highRes))
# Crop highRes
highRes_c <- crop(highRes, lowRes_c)
# Convert lowres to SpatVector for terra::extract()
v_poly <- as.polygons(lowRes)
# Convert lowres to sf for exactextractr::exact_extract()
sf_poly <- v_poly |>
st_as_sf()
# Extract values from highRes_c using sf_poly (dataframe)
st <- Sys.time()
e <- extract(highRes_c, v_poly, weights = TRUE, exact = TRUE)
Sys.time() - st
# Time difference of 1.303376 hours
nrow(e)
# 128510304
head(e)
# ID lyr.1 weight
# 1 1 8 0.08010511
# 2 1 3 0.13965326
# 3 1 7 0.13965326
# 4 1 9 0.13965326
# 5 1 2 0.13965326
# 6 1 7 0.13965326
# Extract values from highRes_c using sf_poly (list object)
st <- Sys.time()
e1 <- exact_extract(highRes_c, sf_poly, weights = "area")
Sys.time() - st
# Time difference of 4.259271 mins
head(e1[1][[1]])
# value weight coverage_fraction
# 1 8 173939.9 0.08007456
# 2 3 173939.9 0.13959999
# 3 7 173939.9 0.13959999
# 4 9 173939.9 0.13959999
# 5 2 173939.9 0.13959999
# 6 7 173939.9 0.13959999
Обратите внимание, что значения различаются между двумя методами. Разработчики exactextractr заявляют, что этот метод более точен, чем другие подходы, подробности см. в разделе Точность экстрактора GitHub.
Это здорово - спасибо!
Вы можете ускорить процесс, сначала векторизовав весь растр низкого разрешения, а затем извлекая его за один раз:
## set distinct values for the lowRes cells, otherwise only one polygon
## (the bounding box) will be returned:
values(lowRes) <- 1:ncell(lowRes)
## vectorize:
p_low <- as.polygons(lowRes)
## extract:
res <- extract(highRes, p_low, weights = TRUE, exact = TRUE)
> head(res)
ID lyr.1 weight
1 1 10 0.2408351
2 1 2 0.3010438
3 1 1 0.3010438
4 1 9 0.3010438
5 1 3 0.3010438
6 1 10 0.3010438
На моем стандартном ноутбуке это примерно в 4,5 раза быстрее, чем при использовании поячеечного подхода (при этом для данных вашего примера все равно требуется около 45 с).
Если вас интересует только совокупное значение для каждой ячейки с низким разрешением (при сохранении веса), resampleing может быть гораздо более быстрой альтернативой.
Вы видели этот ответ? gis.stackexchange.com/a/423700/156904