## Eduard F. Martinez-Gonzalez
## Workshop Cienfi: Introduccion a Datos Raster en R
## lee: data/night_light_201301.tif, data/night_light_202301.tif, data/manzanas_cali.rds
## produce: output/cali_recorte.tif
## nota: abra el proyecto workshop-2.Rproj para que las rutas relativas funcionen

## Instalar/llamar las librerías de la clase
require(pacman)
p_load(tidyverse,
       sf,
       osmdata,
       terra,
       mapview)

## Leer un archivo raster
nl_r <- rast("data/night_light_202301.tif")

## Información general
nl_r          

## Visualización
plot(nl_r)

## Obtener los límites administrativos de Cali desde OpenStreetMap
cali <- getbb(place_name = "Cali, Colombia", 
              featuretype = "boundary:administrative", 
              format_out = "sf_polygon")

## Seleccionar la segunda geometría
cali <- cali[2, ]

## Visualizar la geometría de Cali sobre un mapa interactivo
mapview(cali)

## Recortar el raster original (nl_r) usando el polígono de Cali
cali_r <- crop(nl_r , cali , mask=T)

## Visualizar el raster recortado de Cali
mapview(cali_r)

## Guardar el raster recortado como archivo GeoTIFF en el disco
writeRaster(cali_r, "output/cali_recorte.tif", overwrite=T)

## Convertir el raster recortado en un data.frame, incluyendo las coordenadas (x, y)
df <- as.data.frame(cali_r, xy = TRUE)

# Leer archivo RDS que contiene los polígonos de manzanas censales de Cali
mnz <- read_rds("data/manzanas_cali.rds") %>%
       st_as_sf()

## Visualizar las primeras filas del objeto de manzanas
head(mnz)

## raster to sf
cali_sf <- as.polygons(cali_r) %>% st_as_sf()
st_crs(cali_sf) == st_crs(mnz)

## transformar el CRS
cali_sf <- st_transform(cali_sf , st_crs(mnz))
st_crs(cali_sf) == st_crs(mnz)

## unir los shapefiles
mnz_nl <- st_join(x=mnz , y=cali_sf , largest=T)

## Visualizar resultado
mapview(mnz_nl, zcol = "night_light_202301")

## Convertir el objeto sf en un data.frame plano
df_nl <- as.data.frame(mnz_nl, geometry = NULL)
head(df_nl)

## Graficar relación entre población y luminosidad
ggplot(data = df_nl, 
       aes(x = log(personas), y = log(night_light_202301))) +
geom_point() +
geom_smooth(method="lm", se =T , color = "blue") +
theme_bw()

## regresion
lm(asinh(night_light_202301) ~ asinh(personas), data = df_nl)

## luminosidad 2023
nl_r <- c(rast("data/night_light_201301.tif"),
          rast("data/night_light_202301.tif")) %>%
        crop(cali,mask=T) %>%
        as.polygons() %>%
        st_as_sf() %>%
        st_transform(crs=st_crs(mnz))
data <- st_join(nl_r , mnz , largest = T) %>%
        as.data.frame(xy=T)