Я пытаюсь уменьшить климатические условия, используя методологию из этой статьи, используя программное обеспечение R.Я почти у цели, но мне не хватает пары шагов
Необходимые пакеты и данные
Для этого примера я загрузил некоторые данные на сайт archive.org для загрузки необходимых пакетов и используемых данныхв этом примере используйте следующий код:
library(raster)
library(rgdal)
download.file("https://archive.org/download/Downscaling/BatPatagonia.rds", "Bat.rds")
download.file("https://archive.org/download/Downscaling/TempMinPatNow.rds", "Tmin.rds")
BatPatagonia <- readRDS("Bat.rds")
TempMinPatNow <- readRDS("Tmin.rds")
BatPatagonia - это растровый файл с батиметрией и высотой области, извлеченной и преобразованной из набора данных GEBCO, тогда как TempMinPatNow является минимальной температурой той же области дляЯнварь извлечен из WorldClim.Графики наборов данных показаны ниже:
Цель этого вопроса
Чтобы уменьшить масштаб данных за последниеледниковый максимум Мне нужно смоделировать, каким будет нынешний климат, если бы уровень моря был таким же, как в прошлом.Для того, чтобы сделать это, я использую данные GEBCO и выясняю, где более или менее было побережье.В соответствии с методологией, приведенной в статье, приведенной выше, это первые три шага, которые необходимо выполнить:
- Создание матрицы высот для земли выше 20 метров над уровнем моря
- Расчет множественной линейной регрессиидвижущееся окно
- Экстраполировать коэффициенты на океан
Точка 3 - это то, что я изо всех сил пытался сделать, я покажу, как я сделал первые 2 пунктаи покажите, что я искал, пытаясь решить пункт 3
1.Создайте DEM для земли на высоте 20 метров над уровнем моря
. Для этого я взял растр BatPatagonia и заменил все значения ниже 20 метров на значения NA, используя следующий код:
Elev20 <- BatPatagonia
values(Elev20) <- ifelse(values(Elev20) <= 20, NA, values(Elev20))
Результирующий растр показан на следующем рисунке
2.Вычисление множественной линейной регрессии в движущемся окне
Согласно рукописи на стр. 2591, следующим шагом является создание множественной линейной регрессии в движущемся окне с использованием следующей формулы для высот более 20 метров:
У нас уже есть данные высот, но нам также нужны растры для широты и долготы, для этого мы используем следующий код, где мы сначала создаем Локатори растровые долготы:
Latitud <- BatPatagonia
Longitud <- BatPatagonia
data_matrix <- raster::xyFromCell(BatPatagonia, 1:ncell(BatPatagonia))
values(Latitud) <- data_matrix[, 2]
values(Longitud) <- data_matrix[, 1]
Мы умножим это на растровую маску областей, которые имеют высоты более 20 метров, так что мы получим только те значения, которые нам нужны:
Elev20Mask <- BatPatagonia
values(Elev20Mask) <- ifelse(values(Elev20Mask) <= 20, NA, 1)
Longitud <- Elev20Mask*Longitud
Latitud <- Elev20Mask*Latitud
Теперь я создам стек с переменными ответа и переменными предиктора:
Preds <- stack(Elev20, Longitud, Latitud, TempMinPatNow)
names(Preds) <- c("Elev", "Longitud", "Latitud", "Tmin")
Полученный стек показан на следующем рисунке:
Как указано в документе, движущееся окно должно иметь размер 25 на 25 ячеек, в результате чего получается 625 ячеек, однако они утверждают, что если у движущегося окна меньшечем 170 ячеек с данными, регрессия не должна быть выполнена, и она должна иметь максимум 624 ячеек, чтобы гарантировать, что мы только моделируем области близко к побережью.Результатом этой множественной регрессии с движущимся окном должен быть стек с локальным перехватом и локальной оценкой каждой из бета-версий, представленных в уравнении, показанном выше.Я узнал, как сделать это, используя следующий код с использованием функции getValuesFocal (этот цикл занимает некоторое время):
# First we establish the 25 by 25 window
w <- c(25, 25)
# Then we create the empty layers for the resulting stack
intercept <- Preds[[1]]
intercept[] <- NA
elevationEst <- intercept
latitudeEst <- intercept
longitudeEst <- intercept
Теперь мы запускаем код:
for (rl in 1:nrow(Preds)) {
v <- getValuesFocal(Preds[[1:4]], row = rl, nrows = 1, ngb = w, array = FALSE)
int <- rep(NA, nrow(v[[1]]))
x1 <- rep(NA, nrow(v[[1]]))
x2 <- rep(NA, nrow(v[[1]]))
x3 <- rep(NA, nrow(v[[1]]))
x4 <- rep(NA, nrow(v[[1]]))
for (i in 1:nrow(v[[1]])) {
xy <- na.omit(data.frame(x1 = v[[1]][i, ], x2 = v[[2]][i, ], x3 = v[[3]][i,
], y = v[[4]][i, ]))
if (nrow(xy) > 170 & nrow(xy) <= 624) {
coefs <- coefficients(lm(as.numeric(xy$y) ~ as.numeric(xy$x1) +
as.numeric(xy$x2) + as.numeric(xy$x3)))
int[i] <- coefs[1]
x1[i] <- coefs[2]
x2[i] <- coefs[3]
x3[i] <- coefs[4]
} else {
int[i] <- NA
x1[i] <- NA
x2[i] <- NA
x3[i] <- NA
}
}
intercept[rl, ] <- int
elevationEst[rl, ] <- x1
longitudeEst[rl, ] <- x2
latitudeEst[rl, ] <- x3
message(paste(rl, "of", nrow(Preds), "ready"))
}
Coeffs <- stack(intercept, elevationEst, latitudeEst, longitudeEst, (intercept + Preds$Elev * elevationEst + Preds$Longitud * longitudeEst + Preds$Latitud *latitudeEst), Preds$Tmin)
names(Coeffs) <- c("intercept", "elevationEst", "longitudeEst", "latitudeEst", "fitted", "Observed")
результат этого цикла - стек coeffs, показанный ниже:
Вот где я застрял:
Экстраполируем коэффициентык океану
Теперь цель состоит в том, чтобы экстраполировать первые 4 растра из стека Coeffs (intercept, elevationEst, longitudeEst и latitudeEst) туда, где должно быть побережье в соответствии с последним ледниковым максимумом, который был на 120 метров меньше
MaxGlacier <- BatPatagonia
values(MaxGlacier) <- ifelse(values(MaxGlacier) < -120, NA,1)
Спроектированная береговая линия показана на следующей карте:
Способ, которым авторы проецировали коэффициенты на побережье, заключался в заполнении пробелов с помощью решения уравнения Пуассона с использованием poisson_grid_fill языка NCL из NCAR .Но я хотел бы сделать это простым и попытаться сделать все на одном языке.Я также нашел похожую функцию в python .
Я был бы рад любому процессу экстраполяции, который работает хорошо, я не ограничиваю свой поиск этим алгоритмом.
Iнашел несколько r пакетов, которые заполняют пробелы, такие как Gapfill , и даже нашел обзор методов для заполнения пробелов, но большинство из них предназначены для интерполяции и в основном для слоев NDVI, которые могут бытьоснованный на других слоях, где заполнен пробел.
Есть идеи о том, как сдвинуться с места?
Спасибо