Unidad4.qmd

---
title: "Unidad 4: Interpolación espacial y Geoestadística"
---

- Interpolación basada en métodos determinísticos
  - Vecino más cercano (Voronoi)
  - Inverse distance weighting (IDW)

- Predicción:
  - Covariables
  - Modelos regresivos
    - Regresión lineal con desface espacial
    - Regresión ponderada geográficamente
  - Métodos geoestadísticos
    - Kriging


## Paquetes utilizados en esta lección

```{r}
#| warning: false
#| message: false
library(sf)
library(tidyverse)
library(tmap)
library(raster)
library(stars)
library(gstat)
library(spgwr)
library(dismo)
library(fields)
library(automap)
```


## Actividad 1: Interpolación basada en métodos determinísticos

> Estos ejercicios se realizaron junto con Pablo Tenelanda Y Damián Tuba.

**Cargue las capas Temperatura_2.tif y recórtela hasta 200 metros más que la extensión de la capa muestreo.gpkg:: Puntos de muestreo. Y a partir de esos datos genere:**

```{r}
puntos <- st_read("data/gpkg/muestreo.gpkg", "Puntos de muestreo", quiet = TRUE)

puntos <- puntos |> 
  mutate(x = st_coordinates(puntos)[, "X"], y = st_coordinates(puntos)[, "Y"]) |> 
  as_Spatial()

temp2 <- raster("data/tif/Temperatura_2.tif")
```

Creamos el buffer de al rededor de los puntos:

```{r}
mascara <- st_bbox(puntos) |> 
  st_as_sfc(crs = st_crs(puntos)) |> #convertir el bbox a polígono
  st_buffer(200) |>  # Crear un buffer de 200 m alrededor del bbox
  as_Spatial() 
```

Recortamos la capa de temperatura con los límites de los puntos de muestreo:

::: {.column width="48%"}

```{r}
template <- temp2 %>% 
  crop(mascara) %>% 
  as("SpatialPixelsDataFrame")
```

Revisamos el resultado

```{r}
#| eval: false
template |> plot()
puntos |> plot(add = T)
```

:::

::: {.column width="48%"}

```{r}
#| echo: false
template |> plot()
puntos |> plot(add = T)
```

:::


- Modelo de tendencia espacial de los datos mediante un polinomio de segundo orden (modelo cuadrático). Puede usar la fórmula: variable ~ x + y + I(x*y) + I(x^2) + I(y^2)

```{r}
spatial_model <- krige(
  ta_media ~ x + y + I(x*y)+ I(x^2)+I(y^2), puntos, template
)
spplot(spatial_model ,"var1.pred",
       main= "Tendencia (2do orden)")
```

- Diagrama de Voronoi.

A continuación se presenta los diagramas de voronoi con el paquete dismo y con una interpolación de kriging simple.

```{r}
#| layout-ncol: 2
dismo::voronoi(coordinates(puntos)) %>% 
  plot(main = "Polígonos de Voronoi")

krige(ta_media ~ 1, puntos, template, nmax=1, set=list(idp = 0), debug.level = 0) |>
  spplot("var1.pred", main= "1 vecino próximo")
```

- Interpolación ponderada por el inverso de la distancia

```{r}
idw(
  ta_media ~ 1, puntos, template, nmax = 5, idp = 2, debug.level = 0
) |> spplot("var1.pred", main= "Inverse Distance Weighting")
```


- Interpolación Thin Plate

Con la función **Tps** y las coordenadas de la capa puntos como la matriz de variable independientes se crea el modelo, que luego es ingresado junto al raster **template** en la función **interpolate**.

```{r}

tps_model <- fields::Tps(coordinates(puntos), puntos$ta_media)
tps_grid <- interpolate(raster(template), tps_model)
spplot(tps_grid, main = "Thin plate spline")
```


-----------------------------------------------------------------------

## Actividad 2: Predicción Modelos regresivos

**Cargue las capas modis_perc_bosque_2010_clip.tif y dem200s.tif. Recórtelas tal como lo hizo con la capa de temperatura.**

- Ajuste un modelo de regresión geográficamente ponderada (GWR) con las variables "ta_media" como dependiente y "Altitud" + "Forest_P_2010" como independientes


```{r}
alt <- raster("data/tif/dem200s.tif") |>
  crop(mascara) 
#Llamar al archivo modis
forest2010 <- raster("data/tif/modis_perc_bosque_2010_clip.tif") |>
  crop(mascara) 
```


```{r}
bwG <- gwr.sel(ta_media ~ Altitud + Forest_P_2010,
               data = puntos,
               gweight = gwr.Gauss,
               verbose = FALSE)
bwG

# Ajuste del modelo sin predicción

gwr.ta <- gwr(ta_media ~ Altitud + Forest_P_2010,
              data = puntos,
              bandwidth = bwG,
              gweight = gwr.Gauss,
              fit.points = template)

gwr.ta

```

::: {.column width="45%"}

```{r}

template2 <- stack(list(ta_media = raster(template),
                         Altitud = alt,
                         Forest_P_2010 = forest2010)) |>
  mask(alt) |> 
  as("SpatialPixelsDataFrame")

gwr.ta2 <- gwr(ta_media ~ Altitud + Forest_P_2010,
               data = puntos,
               bandwidth = bwG,
               gweight = gwr.Gauss,
               fit.points = template2,
               predictions = TRUE)
```

:::

::: {.column width="50%"}

```{r}
spplot(gwr.ta2$SDF, "pred")
```

:::


- Describa el resultado.

> El modelo GWR, es una de las tecnicas de regresión espacial empleado ampliamente en geografia y en otras areas. El cual mediante una ecuación de regresión evalua las variables que influyen en el proceso que se intenta comprender o precedir. Obteniendo como resultado los coeficientes de correlación entre el proceso a comprender o predecir y cada variable evaluada. 
Por lo tanto, considerando a la variable temperatura como dependiente de la altitud y Forest  2010. Se obtuvo que la altitud presenta coeficeintes negativos (-0.005: -0.003), mientras que la cobertura de bosque, presenta coeficientes entre -0.001 a 0.1. 
En ambos casos, su influencia es baja con respecto a la temperatura. 


- Genere las capas de valores beta.

```{r}
gwr.ta.stk <- raster::stack(gwr.ta$SDF)
plot(gwr.ta.stk)
```

- Mediante álgebra de mapas genere la predicción para la variable de temperatura, usando las capas beta y las variables independientes .

```{r}
ta_pred <- gwr.ta.stk[[2]] + alt * gwr.ta.stk[[3]] + forest2010 * gwr.ta.stk[[4]]

spplot(ta_pred)
```


------------------------------------------------------------------------

## Actividad 3. Predicción Geoestadística

- A partir de los mismos datos de muestreo genere el semivariograma experimental.

```{r}
v <- variogram(ta_media ~ 1, locations = puntos)
plot(v)
```


- Ajuste un modelo teórico que más considere pertinente al semivariograma, y justifique por qué utilizó ese Modelo teórico.

Para ajustar el modelo vamos a revisar los disponibles

```{r}
gstat::show.vgms()
```

Entre los que mejor se ajustan al semivariograma experimental están "Sph", "Cir" y "Pen"

```{r}
#| layout-ncol: 3
plot(v, model = fit.variogram(v, model = vgm(model = "Sph")))
plot(v, model = fit.variogram(v, model = vgm(model = "Cir")))
plot(v, model = fit.variogram(v, model = vgm(model = "Pen")))
```

El mejor modelo de variograma que se ajusta, tanto para distancias cortas como para distancias largas, es el modelo "Pen"

```{r}
fit.v <- fit.variogram(v, model = vgm(model = "Pen"))
plot(v, model = fit.v)
```

- Genere los resultados de predicción con los tres tipos de Krigin (Va a necesitar ajustar un semivariograma adicional)

> Kriging simple

```{r}
simple_krige <- krige(ta_media ~ 1, puntos, template, model = fit.v, 
            beta = mean(puntos$ta_media))
```

> Kriging Ordinario

```{r}
ordinary_krige <- krige(ta_media ~ 1, locations = puntos, newdata = template,
            model = fit.v)
```

> Kriging Universal

```{r}
universal_krige <- krige(ta_media ~ x + y, puntos, template, 
            model = automap::autofitVariogram(ta_media ~ x + y, puntos)$var_model, 
            block = c(100, 100))
```



```{r}
#| layout-ncol: 3
spplot(simple_krige, "var1.pred", main = "Simple")
spplot(ordinary_krige, "var1.pred", main = "Ordinario")
spplot(universal_krige, "var1.pred", main = "Universal")
```