Широта в операции вычисления растра

Question

Я пытаюсь создать карту мира Коппена, используя данные http://worldclim.org. Чтобы найти правильный климат Коппена, мне нужны данные об осадках и температуре (у меня есть одна растровая карта на каждый месяц для каждой из этих двух переменных)и широта.

Я попытался сделать следующее:

prast <- list.files(path = "prec25/", pattern = glob2rx('*.tif'), full.names = T)
trast <- list.files(path = "temp25/", pattern = glob2rx('*.tif'), full.names = T)
lrast <- c(prast, trast)
climrast <- stack(lrast)

koppen_map <- calc(climrast, filename = "koppen.tif", fun = function(x) koppen(x[13:24], x[1:12], yFromCell(climrast, x[1])))

climrast - это RasterStack с 24 различными слоями (12 слоев с данными о температуре и 12 слоев с данными об осадках).Функция koppen нуждается в векторе с 12 значениями для температуры (это будет x[13:24]) и 12 значениями для температуры (x[1:12]).

yFromCell(climrast, x[1]) должен дать мне широту, но calc операция завершается неудачно, потому что yFromCell(climrast, x[1]) возвращает NA в некоторых случаях.

Если я заменю yFromCell(climrast, x[1]) на произвольное число, например 10, операция calc будет работать нормально.

Любая идея, чтоЯ делаю не так?

Answer 1

Безопасный для памяти (и простой) способ получить RasterLayer со значениями широты, вы можете сделать:

x <- init(climrast, 'y')

Рабочий пример с данными worldclim:

library(raster)
prast <- getData('worldclim', var='prec', res=10)
tmin <- getData('worldclim', var='tmin', res=10)
tmax <- getData('worldclim', var='tmin', res=10)
trast <- (tmin + tmax) / 2

lat <- init(trast, 'y')

lrast <- stack(prast, trast, lat)
climrast <- crop(lrast, extent(25,30,-5,0))

# example function
koppen <- function(temp, prec, lat) {
    (sum(temp * prec) + lat) / 1000
}

koppen_map <- calc(climrast, filename = "koppen.tif", fun = function(x) koppen(x[13:24], x[1:12], x[25]), overwrite=TRUE)

Answer 2

В вашем calc вы передаете x[1] на yFromCell.Но x[1] - это значение растровой ячейки, тогда как вам нужно передать номер ячейки в yFromCell.Я могу проиллюстрировать это на минимальном примере:

Сначала давайте создадим небольшой фиктивный растр

library(raster)
set.seed(0)
clim = raster(matrix(sample(c(1:10,NA), 100, T), 10, 10))

Теперь давайте попробуем получить его широты, используя аналогию с тем, что вы имели в примере

lat = calc(clim, function(x) yFromCell(clim, x))
plot(lat)

Как видите, это совсем не так - мы получили совершенно неправильные значения широты, потому что мы передали значение ячейки, а не номер ячейки.

Итак, давайте создадим растровый слой с правильными широтами

lat = clim
lat[] = yFromCell(clim, 1:ncell(clim))
plot(lat)

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

climrast = stack(list(clim, lat))
koppen = calc(climrast, function(x) x[1]*x[2])