Creación de mapas en 3D

## definimos nuestra dirección wd dónde trabajaremos con los archivos del INEGI

setwd("C:/mapas_dina_mex")
## librerías

library(tmap) #Dibujar el mapa
library(sf) #Para leer el shapefile y reducir el tamaño del archivo
library(pryr) #Calcular el tamaño del archivo
library(readr) #para cargar csv
library(base) # para merge
library(reshape2) # para hacer dcast
library(dplyr) # para inner join

base con .shp file

Aquí es importante poner la capa ó layer y sobre todo el encoding. De otra forma no leera bien los acentos de los nombres de los municipios. Antes se usaba la librería “rgdal”, ahora con sf es posible cargar el .shp file

# importamos shp de municipios 
# el shp lo descargamos del marco geoestadistico del INEGI https://www.inegi.org.mx/app/biblioteca/ficha.html?upc=889463776079
shp_municipios <- st_read("mg_sep2019_integrado/conjunto_de_datos/00mun.shp", layer = "00mun", options = "ENCODING=WINDOWS-1252")
## options:        ENCODING=WINDOWS-1252 
## Reading layer `00mun' from data source 
##   `C:\mapas_dina_mex\mg_sep2019_integrado\conjunto_de_datos\00mun.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2465 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 911292 ymin: 319149.1 xmax: 4082997 ymax: 2349615
## Projected CRS: MEXICO_ITRF_2008_LCC

base con algunos datos sociodemográficos

# podemos descargar esta carpeta de <https://www.inegi.org.mx/programas/ccpv/2020/#Datos_abiertos>

iter_00_cpv2020 <- read_csv("C:/mapas_dina_mex/conjunto_de_datos_iter_00CSV20.csv")


# Volvemos numericas las columnas que nos interesan
# las primeras 10 no son interesantes aquí porqué son datos nacionales...
iter_00_cpv2020[, 10:230] <- sapply(iter_00_cpv2020[, 10:230],
as.numeric)
#creamos claves de estado y municipio mediante los filtros con ifelse

iter_00_cpv2020$ENTIDAD <- ifelse(iter_00_cpv2020$ENTIDAD<10,
                                  paste(0,iter_00_cpv2020$ENTIDAD,sep=""), 
                                  iter_00_cpv2020$ENTIDAD)

iter_00_cpv2020$MUN <- ifelse(iter_00_cpv2020$MUN<10,
                              paste("00",iter_00_cpv2020$MUN,sep=""),
                              ifelse(iter_00_cpv2020$MUN>=10
                                     & iter_00_cpv2020$MUN<100 ,
                                     paste("0",iter_00_cpv2020$MUN,sep=""),
                                     iter_00_cpv2020$MUN ))

iter_00_cpv2020$CVEGEO <- paste(iter_00_cpv2020$ENTIDAD,iter_00_cpv2020$MUN,
                                sep="")

# base dejando solo municipios

iter_00_cpv2020_mun <- iter_00_cpv2020[iter_00_cpv2020$NOM_LOC=="Total del Municipio",]

La Magia: Merge del shp y el df

Aquí es importante explorar y entender los merge… Éste quizá no sea el más apropiado pero funcionó para mapear las coordenadas que usaremos, por eso fue muuuuuuy importante cargar el encoding correcto, pues hice el merge con los nombres de los municipios.

shp_municipios_ <-merge(shp_municipios,iter_00_cpv2020_mun,by.x="NOMGEO",
                        by.y="NOM_MUN")
object_size(shp_municipios_) 
## 66.60 MB
#View(shp_municipios_)

Ahora creamos nuestro gráfico interactivo

#mapeamos población total

tm_shape(shp_municipios_) +
  tm_fill("POBTOT", id="NOMGEO", palette ="RdPu", 
          style="quantile", title="Pob. total") +
  tm_borders("cyan",alpha=.07) +
  tm_layout("Población total por municipio en 2020",
            main.title.position = "center") +
  tm_view(view.legend.position =c("left","bottom")) #esto sirve para bajar la 

#leyenda y que no estorbe con el título 

tmap_last <- tmap_last()

# cambiamos la vista a una donde se muestre un mapa dinámico

tmap_mode("view")
## tmap mode set to interactive viewing

Extra: gráfica de pobtotal por Edo.

#################################################################### 

################# hacer gráfica de pobtotal por edo

###################################################################### 

library("ggplot2")

pop_AGS_2020 <- iter_00_cpv2020$POBTOT[4] 
pop_EDOS_2020 <-iter_00_cpv2020 %>% filter(NOM_LOC == 'Total de la Entidad') %>%
select(POBTOT, NOM_ENT)

pop_AGS_df <- as.data.frame(pop_EDOS_2020) 
typeof(pop_AGS_df[,2])
## [1] "character"
typeof(pop_AGS_df[1,]) 
## [1] "list"
#View(pop_AGS_df)
pop_AGS_df$POBTOT_Millon <- pop_AGS_df$POBTOT / 1e6 
# Graficar con ggplot2 
pob_Ent <- ggplot(pop_AGS_df, aes(x = POBTOT_Millon, y =
NOM_ENT, colours(class))) + geom_point() + labs(title = "Población Total
por Entidad en 2020", x = "Población Total (en millones)", y ="Entidad") +
  theme_bw() + 
  geom_text(data = subset(pop_AGS_df, NOM_ENT %in% c ("Aguascalientes",
                                                      "Quintana Roo", 
                                                      "Yucatán")), 
            aes(label = paste(POBTOT_Millon)), hjust = -0.1, vjust = 0.5,
            size = 3.5, color = "cyan")
pob_Ent

######### salvar la gráfica

#ggsave("pon_ent_mex.png", plot = pob_Ent, width = 6, height = 4, dpi =300)

Ahora sí nuestro Gráfico 3D!!!!!!

############################################ 

######## RAYSHADER

########################################### 

library("ggplot2") 
#install.packages("viridis") 
library("viridis")
## Loading required package: viridisLite
#install.packages("viridisLite") 
#install.packages("tidyverse")
#install.packages("sf")
#install.packages("rayshader")
#install.packages("magick") 
library("magick") 
## Linking to ImageMagick 6.9.12.96
## Enabled features: cairo, freetype, fftw, ghostscript, heic, lcms, pango, raw, rsvg, webp
## Disabled features: fontconfig, x11
#install.packages("av")
library("av") 
library("sf")

### Transformar a 3D

library("rayshader")

# Primero hacemos un ggplot con nuestros datos sobre municipios

ggshp <- ggplot(data = shp_municipios_) + geom_sf(aes(fill =
POBTOT)) + scale_fill_viridis() + # scale_fill_viridis() en lugar de scale_color_viridis()
  ggtitle("Población México") + 
  theme_bw()

# Imprimir el ggplot

print(ggshp)

# Transformar a 3D

plot_gg(ggshp, multicore = TRUE, width = 7, height = 7, scale = 600,
        windowsize=c(1280,720),
        zoom = 0.65, phi = 50)
render_snapshot()

### Crear un vídeo
#Video <- file.path("C:/mapas_dina_mex", "pop_m.mp4")
#render_movie(filename = Video)