Archivos necesarios

El problema

Antes de empezar leeré los archivos y cargaré los archivos necesarios

library(tidyverse)
library(raster)
library(rgdal)
TempMinJanLGMPat <- read_rds("TempMinJanLGMPat.rds")
TempMinPatNow <- read_rds("TempMinPatNow.rds")
BatPatagonia <- read_rds("BatPatagonia.rds")

La idea de este trabajo que te mando es el hacer un downscaling de las capas climáticas del PaleoView (Fordham et al. 2017), las cuales están en una resolución muy grosera, para esto, me estoy basando en el trabajo linkeado acá (Schmatz et al. 2015), a modo de ejemplo usaré un pedazo de Sudamérica que elegí al azar sin ningún sesgo alguno y que puedes ver en la figura 1.

Figura 1. Temperatura mínima del mes de Enero en un lugar elejido al azar en el tiempo actual

Figura 1. Temperatura mínima del mes de Enero en un lugar elejido al azar en el tiempo actual

Como base tomé el mismo espacio y desde paleoView saqué la información del mismo mes para el último máximo glacial en 21.000 años antes del presente, con la diferencia en grados desde el presente (figura 2).

Figura 2. Temperatura mínima del mes de Eneto en el mismo lugar aleatorio hace 21.000 años.

Figura 2. Temperatura mínima del mes de Eneto en el mismo lugar aleatorio hace 21.000 años.

La figura 2 tiene una resolución mucho mas baja, para eso usaremos el método delta acoplado con correcciones para las masas de tierra observadas en el paper ya citado (Schmatz et al. 2015).

Para esto lo primero que necesitamos es una capa de la misma área con la batimetría y altura de la zona. Para la cual obtuve de la base de datos GEBCO (Weatherall et al. 2015) (figura 3).

Figura 3. Batimetría y altitud de patagonia.

Figura 3. Batimetría y altitud de patagonia.

De acuerdo a la literatura el nivel del mar en el último máximo glacial era 120 metros menor a la actualidad (Fleming et al. 1998, @milne2005modelling), para generar, podemos usar la base de datos de batimetría para hacer una estimación de como habría sido la linea de costa en ese momento, lo cual vemos en el siguiente código y en la figura 4.

MaxGlacier <- BatPatagonia

values(MaxGlacier) <- ifelse(values(MaxGlacier) < -120, NA, values(MaxGlacier))
Figura 4. Linea de costa hace 21.000 años

Figura 4. Linea de costa hace 21.000 años

Esto implica que hay que extender el clima base (figura 1) hasta la costa que tiene la figura 4 para poder hacer el método delta. En el paper que estoy siguiendo, lo primero que hacen es hacer una regresión de ventana móvil por el raster, donde ven la Temperatura explicada por altura, latitud, y longitud:

\[Temp = \beta_0 + \beta_1 Altura + \beta_2Latitud + \beta_3logitud\] Esta regresión se hace en una ventana móvil de 25 por 25 celdas, donde se generarán un raster de cada parámetro \(\beta_0\), \(\beta_1\), \(\beta_2\) y \(\beta_3\). Para hacer esto lo primero que necesitamos es armar un stack con las tres variables, pero antes de eso crearé una mascara del área de estudio:

Mask <- TempMinPatNow

values(Mask) <- ifelse(is.na(values(Mask)), NA, 1)

A continuación en la figura 5 vemos la máscara que usaremos para

Figura 5. Máscara del área a trabajar

Figura 5. Máscara del área a trabajar

Altura

El siguiente código es el de la altura.

Elev <- BatPatagonia

values(Elev) <- ifelse(values(Elev) < 0, NA, values(Elev))

Latitud

El siguiente es el código de la latitud.

LatitudTot <- BatPatagonia

data_matrix <- raster::xyFromCell(BatPatagonia, 1:ncell(BatPatagonia))

values(LatitudTot) <- data_matrix[, 2]
Latitud <- LatitudTot * Mask

Latitud

El siguiente es el código de la longitud.

LongitudTot <- BatPatagonia

values(LongitudTot) <- data_matrix[, 1]
Longitud <- LongitudTot * Mask

Stack de predictores

Desde esto creamos el stack de predictores y la temperatura mínima como vemos en la figura 6.

Preds <- stack(Elev, Longitud, Latitud, TempMinPatNow)

names(Preds) <- c("Elev", "Longitud", "Latitud", "Tmin")
plot(Preds, colNA = "black")
Figura 6. Predictores par el modelo

Figura 6. Predictores par el modelo

Análisis de regresiones con ventana móvil

Siguiendo el manuscrito generaremos una regresión de ventana móvil que pase por cada una de las celdas con un cuadrado de 25 por 25 celdas, con la equación 1 que vemos más arriba, para eso lo primero es establecer w que es la ventana

# Establecemos la ventana de 25 por 25

w <- c(25, 25)

# Luego generamos las capas para los estimadores

intercept <- Preds[[1]]
intercept[] <- NA

elevationEst <- intercept

latitudeEst <- intercept

longitudeEst <- intercept

la función que pongo acá abajo es la que me permite crear la regresión movil de cada 25 por 25 celdas. En la figura 7 vemos el resultado de las ecuaciones, donde vemos primero el raster de los interceptos, seguido del de la pendiente de la elevación, de la pendiente de la longitud, de la latiitud, y finalmente la predicción del modelo y el observado. Todo esto calculado en el siguiente código, el cuál asegura tal como pide el paper que solo haga los calculos cuando para cada selga se tenga en la ventana movil al menos 175 celdas en la ventana móvil.

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, ]))
        # Esto nos dice que hay almenos 170 celdas no vacias
        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")
Figura 7. Stack con las estimaciones de el intercepto y las pendientes locales, mas la predicción y el observado

Figura 7. Stack con las estimaciones de el intercepto y las pendientes locales, mas la predicción y el observado

Limitando el área de la ventana movil.

Según lo que dice el texto del manuscrito en particular en la página 2952 dice que la regresión de ventana movil hay que realizarla solo para sectores mayores a 20 m.s.n.m, para esto crearé una nueva mascara basada en la elevación para excluir lugares con alturas menores a 20 metros y multiplicare lo predictores por esta mascara, el resultado lo vemos en la figura 8:

Elev20Mask <- BatPatagonia

values(Elev20Mask) <- ifelse(values(Elev20Mask) <= 20, NA, 1)

NewPreds <- Preds * Elev20Mask
names(NewPreds) <- c("Elev", "Longitud", "Latitud", "Tmin")
Figura 8. Predictrores sobre 20 msnm

Figura 8. Predictrores sobre 20 msnm

La segunda restricción para generar la capa de regresiones es que la ventana movil tenía que tener un mínimo de 170 celdas que contegan datos, lo cual ya esta incorporado en el código

Problema 1 a resolver

Acá tengo resuelto la regresión de moving window, pero posteriormente, hay que extrapolar los valores de la figura 9 hasta los límites de la figura 4, en el paper usan la solución de la ecuación de poisson via relajación, que no se bien que es.

for (rl in 1:nrow(NewPreds)) {
    v <- getValuesFocal(NewPreds[[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, ]))
        # Esto nos dice que hay almenos 170 celdas no vacias
        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(NewPreds), "ready"))
}

Coeffs <- stack(intercept, elevationEst, latitudeEst, longitudeEst, (intercept + 
    NewPreds$Elev * elevationEst + NewPreds$Longitud * longitudeEst + NewPreds$Latitud * 
    latitudeEst), NewPreds$Tmin)
names(Coeffs) <- c("intercept", "elevationEst", "longitudeEst", "latitudeEst", 
    "fitted", "Observed")
Figura 9. Stack con las estimaciones de el intercepto y las pendientes locales, mas la predicción y el observado, tomando en cuenta el quitar el borde costero

Figura 9. Stack con las estimaciones de el intercepto y las pendientes locales, mas la predicción y el observado, tomando en cuenta el quitar el borde costero

con respuesta de Stackoverflow modificada

gaps = which(is.na(Coeffs[[1]] & is.na(Mask))[])
intercept.ext = Coeffs[[1]]
w = matrix(c(0, 0.25, 0, 0.25, 0, 0.25, 0, 0.25, 0), nc = 3, nr = 3)
max.it = 1000
for (i in 1:max.it) intercept.ext[gaps] = focal(intercept.ext, w = w, na.rm = TRUE)[gaps]
intercept.ext = mask(intercept.ext, MaxGlacier)
saveRDS(intercept.ext, "intercept.ext.rds")
gaps = which(is.na(Coeffs[[1]] & is.na(Mask))[])
elevation.ext = Coeffs[[2]]
w = matrix(c(0, 0.25, 0, 0.25, 0, 0.25, 0, 0.25, 0), nc = 3, nr = 3)
max.it = 1000
for (i in 1:max.it) elevation.ext[gaps] = focal(elevation.ext, w = w, na.rm = TRUE)[gaps]
elevation.ext = mask(elevation.ext, MaxGlacier)
saveRDS(elevation.ext, "elevation.ext.rds")
gaps = which(is.na(Coeffs[[1]] & is.na(Mask))[])
latitude.ext = Coeffs[[3]]
w = matrix(c(0, 0.25, 0, 0.25, 0, 0.25, 0, 0.25, 0), nc = 3, nr = 3)
max.it = 1000
for (i in 1:max.it) latitude.ext[gaps] = focal(latitude.ext, w = w, na.rm = TRUE)[gaps]
latitude.ext = mask(latitude.ext, MaxGlacier)
saveRDS(latitude.ext, "latitude.ext.rds")
gaps = which(is.na(Coeffs[[1]] & is.na(Mask))[])
longitude.ext = Coeffs[[4]]
w = matrix(c(0, 0.25, 0, 0.25, 0, 0.25, 0, 0.25, 0), nc = 3, nr = 3)
max.it = 1000
for (i in 1:max.it) longitude.ext[gaps] = focal(longitude.ext, w = w, na.rm = TRUE)[gaps]
longitude.ext = mask(longitude.ext, MaxGlacier)
saveRDS(longitude.ext, "latitude.ext.rds")
Coeffs2 <- stack(intercept.ext, elevation.ext, latitude.ext, longitude.ext, 
    (intercept.ext + MaxGlacier * elevation.ext + LongitudTot * longitude.ext + 
        LatitudTot * latitude.ext), NewPreds$Tmin)
names(Coeffs2) <- c("intercept", "elevationEst", "longitudeEst", "latitudeEst", 
    "fitted", "Observed")
plot(Coeffs2, colNA = "black")

saveRDS(Coeffs2, "Coeffs2.rds")

Llenar el fitted con los valores anteriores:

TempPred <- Coeffs2$fitted
values(TempPred) <- ifelse(is.na(values(TempPred)), values(TempMinPatNow), values(TempPred))
TempPred <- TempPred/10

Sumar el delta

Primero interpolamos bilinearmente la diferencia:

TempMinJanLGMPat <- resample(TempMinJanLGMPat, Coeffs2$fitted)

lo cual se ve así:

para luego sumar esta capa de delta interpolada bilinearmente a la reconstrucción actual de acuerdo a lo predicho:

TempMinJanLGMPat <- TempPred + TempMinJanLGMPat

Comparemos la capa actual con la pasada:

Referencias

Fleming, Kevin, Paul Johnston, Dan Zwartz, Yusuke Yokoyama, Kurt Lambeck, and John Chappell. 1998. “Refining the Eustatic Sea-Level Curve Since the Last Glacial Maximum Using Far-and Intermediate-Field Sites.” Earth and Planetary Science Letters 163 (1-4). Elsevier: 327–42.

Fordham, Damien A, Frédérik Saltré, Sean Haythorne, Tom ML Wigley, Bette L Otto-Bliesner, Ka Ching Chan, and Barry W Brook. 2017. “PaleoView: A Tool for Generating Continuous Climate Projections Spanning the Last 21 000 Years at Regional and Global Scales.” Ecography 40 (11). Wiley Online Library: 1348–58.

Milne, Glenn A, Antony J Long, and Sophie E Bassett. 2005. “Modelling Holocene Relative Sea-Level Observations from the Caribbean and South America.” Quaternary Science Reviews 24 (10-11). Elsevier: 1183–1202.

Schmatz, D, J Luterbacher, N Zimmermann, and P Pearman. 2015. “Gridded Climate Data from 5 Gcms of the Last Glacial Maximum Downscaled to 30 Arc S for Europe.” Climate of the Past Discussions 11 (3): 2585–2613.

Weatherall, Pauline, Karen M Marks, Martin Jakobsson, Thierry Schmitt, Shin Tani, Jan Erik Arndt, Marzia Rovere, Dale Chayes, Vicki Ferrini, and Rochelle Wigley. 2015. “A New Digital Bathymetric Model of the World’s Oceans.” Earth and Space Science 2 (8). Wiley Online Library: 331–45.