#Cargaremos el CENSO 2010 de México, datos de las alcaldías de CDMX y datos sobre viviendas dañadas en el terremoto del 19 de septiembre de 2019.

#Cargamos el archivo “Diccionario”, que contiene definiciones de variables del censo 2010 de México #Cargamos el archivo del Censo 2010 de México, que posee todas las variables censadas #Cargamos el archivo shapefile con los radios censales (AGEB: Area Geo-estadística Básica) #Cargamos el archivo Shapefile con las Alcaldías de CDMX (16 alcaldías)

CDMX_Diccionario <- read.csv("/Users/martinsingla/Documents/UTDT - MEU/Ciencia de Datos para Ciudades II/INFO SELECCIONADA/CDMX_censo10_AGEB_diccionario.csv",
                             encoding= "UTF-08")
CDMX_Censo <- read.csv("/Users/martinsingla/Documents/UTDT - MEU/Ciencia de Datos para Ciudades II/INFO SELECCIONADA/CDMX_censo10_AGEB.csv", 
                       encoding = "UTF-08",
                       stringsAsFactors = TRUE)
Mexico_AGEB <- read_sf("/Users/martinsingla/Documents/UTDT - MEU/Ciencia de Datos para Ciudades II/INFO SELECCIONADA/Mexico_AGEB/ageb_urb_2009.shp")
CDMX_Alcaldias <- read_sf("/Users/martinsingla/Documents/UTDT - MEU/Ciencia de Datos para Ciudades II/INFO SELECCIONADA/CDMX_Alcaldias/alcaldias.shp")

#Acotamos el Censo Nacional a los radios censales de Ciudad de México, y retiramos todas las filas que denotan valores agregados (alcaldías, municipios y total CDMX), nos quedamos solo con los datos por AGEB. Hacemos las transformaciones necesarias

CDMX_Censo_AGEB <- CDMX_Censo %>%
  filter(nom_loc == "Total AGEB urbana") %>%
  arrange(ageb)

CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 1] <- "001"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 2] <- "002"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 3] <- "003"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 4] <- "004"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 5] <- "005"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 6] <- "006"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 7] <- "007"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 8] <- "008"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 9] <- "009"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 10] <- "010"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 11] <- "011"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 12] <- "012"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 13] <- "013"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 14] <- "014"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 15] <- "015"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 16] <- "016"
CDMX_Censo_AGEB$mun[CDMX_Censo_AGEB$mun == 17] <- "017"

CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 1] <- "0001"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 11] <- "0011"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 15] <- "0015"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 17] <- "0017"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 19] <- "0019"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 20] <- "0020"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 21] <- "0021"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 24] <- "0024"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 26] <- "0026"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 27] <- "0027"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 29] <- "0029"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 33] <- "0033"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 36] <- "0036"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 152] <- "0152"
CDMX_Censo_AGEB$loc[CDMX_Censo_AGEB$loc == 300] <- "0300"

CDMX_Censo_AGEB$AGEB <- paste(CDMX_Censo_AGEB$loc, CDMX_Censo_AGEB$ageb, sep = "")
CDMX_Censo_AGEB$AGEB <- paste(CDMX_Censo_AGEB$mun, CDMX_Censo_AGEB$AGEB, sep = "")

CDMX_Censo_AGEB <- CDMX_Censo_AGEB %>%
  arrange(AGEB)

#Hacemos las modificaciones necesarias en el archvio con los radios censales (AGEB), como filtrar todos los radios censales y quedarnos con los de CDMX, así como extraer las coordenadas de radio censal.

CDMX_AGEB <- Mexico_AGEB %>%
  filter(str_detect(CVEGEO, "^09"))

CDMX_AGEB <- CDMX_AGEB %>%
  mutate(AGEB = str_sub(CVEGEO, -11, -1)) %>%
  arrange(AGEB)

#Unificamos el archivo Shapefile de AGEB de CDMX, con la info censal por AGEB de CDMX

CDMX <- left_join(CDMX_AGEB, CDMX_Censo_AGEB)
## Joining, by = "AGEB"

#Calculamos el area en KM2 de cada polígono

CDMX$area_km <- as.numeric(st_area(CDMX))

#Mapeamos densidad de población por AGEB; superponemos mapa de las alcaldías de CDMX. En el sur de la ciudad, parte de las alcadías estás vacías de radios censales, esto es porque es zona de montaña deshabitada. Por supuesto, la urbanización crece a modo de sprawl por sobre la montaña conforme pasan las décadas. Así mismo, la ciudad continúa hacia el norte (Estado de México, ya no DF).

ggplot() +
  geom_sf(data= CDMX, aes(fill= pobtot /area_km), color= NA) +
  scale_fill_viridis_c() +
  geom_sf(data= CDMX_Alcaldias, fill= NA, color= "grey") +
  labs(title= "CDMX: Densidad poblacional",
       subtitle= "Densidad poblacional por AGEB y delimitación de alcaldías",
       caption= "Fuente: Censo 2010 - INEGI")

#El 19 de Septiembre de 2017 hubo un terremoto de magnitud 7.1 Mw que afectó toda la región central de México. El mismo provocó destrozos, derrumbes y decenas de víctimas fatales en toda la CDMX. En los siguientes Shapefiles se muestran las viviendas (unifamiliares y edificios), afectados por el terremoto en toda la ciudad.

viviendas_uni_afectadas <- read_sf("/Users/martinsingla/Documents/UTDT - MEU/Ciencia de Datos para Ciudades II/INFO SELECCIONADA/reconstruccion-unifamiliares/reconstruccion-viviendas-unifamiliares.shp")
viviendas_multi_afectadas <- read_sf("/Users/martinsingla/Documents/UTDT - MEU/Ciencia de Datos para Ciudades II/INFO SELECCIONADA/reconstruccion-multifamiliares/reconstruccion-multifamiliares-cdmx.shp")

#Realizamos toda una serie de operaciones para seleccionar columnas y contenido de los shapefiles, para poder unirlos en uno.

viviendas_uni_afectadas_1 <- viviendas_uni_afectadas %>%
  select(geometry,
         alcaldia = alcaldia_in,
         tipo_riesgo)
viviendas_uni_afectadas_1$tipo_vivienda <- "unifamiliar/casa"

viviendas_multi_afectadas_1 <- viviendas_multi_afectadas %>%
  select(geometry,
         alcaldia,
         tipo_riesgo)
viviendas_multi_afectadas_1$tipo_vivienda <- "multifamiliar/edificio"

viviendas_afectadas <- rbind(viviendas_multi_afectadas_1, viviendas_uni_afectadas_1)
viviendas_afectadas <- viviendas_afectadas %>%
  arrange(alcaldia)
viviendas_afectadas$alcaldia[viviendas_afectadas$alcaldia == "Álvaro Obregón"] <- "ALVARO OBREGON"

viviendas_afectadas_group <- viviendas_afectadas %>%
  group_by(alcaldia, tipo_vivienda) %>%
  summarise(total = n())

#Graficamos la cantidad de viviendas afectadas, por tipo de vivienda y alcaldía:

ggplot(viviendas_afectadas_group) +
  geom_bar(position = "stack", stat = "identity", 
           aes(x= alcaldia, y= total, fill= tipo_vivienda )) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  labs(title = "CDMX: viviendas afectadas en terremoto",
       subtitle= "Viviendas afectadas en el terremoto del 19 de septiembre de 2017, por tipo de vivienda y alcaldía",
       y= "Total viviendas afectadas",
       x= "",
       caption= "Fuente: Gob. de CDMX. Dispoible en: datos.cdmx.gob.mx")

#Vemos que los barrios más afectados han sido Iztapalapa, Tlahuac, Xochimilco y Milpa Alta, mayoritariamente en viviendas unifamiliares/casas. Así mismo, veremos que estos altos valores absolutos se mantienen cuando consideramos el total de stock de viviendas de la alcaldía.

#Ploteamos estos casos en el mapa:

ggplot() +
  geom_sf(data= CDMX_Alcaldias, fill= NA, color= "grey") +
  geom_sf(data= viviendas_afectadas, aes(color= tipo_vivienda), size = 0.1) +
  labs(title = "CDMX: Daños sísimicos (2017)",
       subtitle= "Viviendas afectadas por el terremoto del 19 de septiembre de 2017",
       caption= "Fuente: Gobierno de la CDMX. Dispoible en: datos.cdmx.gob.mx")

#Unimos los dos Shapefiles

viviendas_afectadas_group1 <- viviendas_afectadas %>%
  group_by(alcaldia) %>%
  summarise(freq = n())

viviendas_afectadas_group1 <- st_set_geometry(viviendas_afectadas_group1, NULL)

CDMX_Alcaldias <- CDMX_Alcaldias %>%
  arrange(nomgeo)
CDMX_Alcaldias$nomgeo <- viviendas_afectadas_group1$alcaldia

viviendas_afectadas_alcaldia <- left_join(CDMX_Alcaldias, viviendas_afectadas_group1, by= c("nomgeo" = "alcaldia"))

CDMX_group <- CDMX %>% 
  group_by(nom_mun) %>%
  summarise(TotPop = sum(pobtot),
            TotViv = sum(vivtot)) %>%
  arrange(nom_mun)
## Warning: Factor `nom_mun` contains implicit NA, consider using
## `forcats::fct_explicit_na`
CDMX_group <- na.omit(CDMX_group)
CDMX_group <- st_set_geometry(CDMX_group, NULL)

CDMX_group$nom_mun <- viviendas_afectadas_group1$alcaldia

viviendas_afectadas_alcaldia <- left_join(viviendas_afectadas_alcaldia, CDMX_group, by= c("nomgeo" = "nom_mun"))

#Mapeamos el gráfico anterior.

ggplot() +
  geom_sf(data= viviendas_afectadas_alcaldia, aes(fill= freq / TotViv)) +
  labs(title= "CDMX: Daños sísimocos (2017)",
       subtitle= "Viviendas dañadas en el terremoto del 19 de septiembre de 2017",
       caption= "Fuente: Gobierno de la CDMX. Dispoible en: datos.cdmx.gob.mx")

#Como vemos en el primer mapa que refleja la densidad demográfica de CDMX, la alcaldía de Iztapalapa (centro-este), posee una densidad demográfica muy alta así como una concentración alta de población en términos absolutos. Por otro lado, en el mapa y gráfico 2, vemos que en esta alcaldía se concentraron la mayor cantidad de casos de daños a viviendas durante el sismo de 2017. Por otro lado, las alcaldías de Tlahuac, Xochimilco y Milpa Alta tienen una baja densidad demográfica y cantidad de viviendas absolutas. Sin embargo, por la altísima cantidad de casos de viviendas dañadas durante el sismo, la razón entre cantidad de viviendas dañadas / stock de viviendas totales, se eleva y aparecen resaltadas en el mapa.