Apuntes

Mapas de España con R

En esta página se mostrará los métodos que se irán aprendiendo por el alumno para la representación de mapas, de momento de España, en R. En ocasiones sólo se presentará la visualización del mapa, no la respresentación de datos sobre él.

1. Archivos shapefile

Aquí se muestran los pasos seguidos por los autores de la web R que R. Enlace a página.

Otros ejemplos:

• Rpubs

Paquetes de R

library(rgdal) # para importar archivos shapefiles
library(broom) # Para transformar los archivos shapefiles 

Mapa de las CCAA

Archivo shapefile

arcgis.com

Comprobar que la carpeta es correcta

file.exists('Comunidades_Autonomas_ETRS89_30N.shp')

## [1] TRUE

Para que readOGR cargue el path el script .rmd o .r tienen que estar en la misma carpeta.

shapefile_ccaa <- readOGR(dsn = "Comunidades_Autonomas_ETRS89_30N.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/albtorval/Documents/Aprendizaje/Trabajos/Mapas/Comunidades_Autonomas_ETRS89_30N.shp", layer: "Comunidades_Autonomas_ETRS89_30N"
## with 19 features
## It has 3 fields

# shapefile_ccaa <- readOGR(dsn = path.expand("Comunidades_Autonomas_ETRS89_30N.shp"))

Otra forma:

library(raster)
s <- shapefile("Comunidades_Autonomas_ETRS89_30N.shp")

Para convertir el archivo shapefile en un dataframe utilizamos la función tidy()

data_ccaa <- tidy(shapefile_ccaa)

Vemos que el dataframe creado no contiene los nombres de las distintas CCAA españolas pero, en cambio, dicho dataframe indica el “id” de las mismas. No obstante, el archivo shapefile original si contiene información de los nombres de las Comunidades Autónomas. Por tanto, podemos extraer los nombres en un nuevo dataframe (nombres_ccaa) y añadir el “id” de cada una de las regiones. El objetivo es, posteriormente, unir los dos dataframes por la columna común “id”.

nombres_ccaa <- tibble(shapefile_ccaa$Texto) %>% 
  mutate(id = as.character(seq(0, nrow(.)-1)))

data_ccaa_mapa <- data_ccaa %>% 
  left_join(nombres_ccaa, by = "id") %>% 
  rename(ccaa = `shapefile_ccaa$Texto`)

Mapa

data_ccaa_mapa %>%
  ggplot() +
  geom_polygon(aes( x= long, y = lat, group = group),
               fill = "black", alpha = 0.8, size = 0.05 ) +
  theme_void() +
  theme(panel.background = element_rect(size= 0.5, color = "white", fill = "white")) +
  labs(title = "CCAA", subtitle = "España") 

Mapa de las provincias

Enlace a página

shapefile_provincias <- readOGR("Provincias_ETRS89_30N.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/albtorval/Documents/Aprendizaje/Trabajos/Mapas/Provincias_ETRS89_30N.shp", layer: "Provincias_ETRS89_30N"
## with 52 features
## It has 5 fields

data_provincias <- tidy(shapefile_provincias)
nombres_provincias <- tibble(shapefile_provincias$Texto) %>% 
  mutate(id = as.character(seq(0, nrow(.)-1)))
data_provincias_mapa <- data_provincias %>% 
  left_join(nombres_provincias, by = "id")

• Mapa:

data_provincias_mapa %>%
  ggplot() +
  geom_polygon(aes( x= long, y = lat, group = group),
               fill = "black", alpha = 0.8, size = 0.05 ) +
  theme_void() +
  theme(panel.background = element_rect(size= 0.5, color = "white", fill = "white")) +
  labs(title = "CCAA", subtitle = "España") 

Mapa de los municipios

Enlace a página

La operativa es similar a la realizada en el punto anterior aunque en lugar de utilizar los nombres de los municipios utilizaremos los códigos provistos por el INE (CODIGOINE).

shapefile_municipios <- readOGR("Municipios_IGN.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/albtorval/Documents/Aprendizaje/Trabajos/Mapas/Municipios_IGN.shp", layer: "Municipios_IGN"
## with 8205 features
## It has 10 fields
## Integer64 fields read as strings:  OBJECTID

data_municipios <- tidy(shapefile_municipios)
codigo_municipios <- tibble(shapefile_municipios$CODIGOINE) %>% 
  mutate(id = as.character(seq(0, nrow(.)-1)))
data_municipios_mapa <- data_municipios %>% 
  left_join(codigo_municipios, by = "id") %>%
  rename(codigo_ine = `shapefile_municipios$CODIGOINE`) %>% 
  mutate(codigo_ine = as.character(codigo_ine)) 

• Mapa:

!!!! TARDA EN EVALUAR

data_municipios_mapa %>%
   ggplot() +
   geom_polygon(aes(x= long, y = lat, group = group),
                fill = "black", alpha = 0.8, size = 0.05 )+
   theme_void() +
   theme(panel.background = element_rect(size= 0.5, color = "white", fill = "white")) +
   labs(title = "CCAA", subtitle = "España")

2. Paquete leaflet

Enlace a página Este paquete funciona de forma muy parecida a la sintaxis de ggplot2.

if (!"leaflet" %in% installed.packages()) install.packages("leaflet")
library(leaflet)

my_map <- leaflet() %>% 
    addTiles() %>%
    addMarkers(lat=36.1848, lng=-5.49, 
               popup="Population explosion")
my_map

Mapa de España

m <- leaflet() %>% setView(lat  = 40.416775, lng = -3.703790, zoom = 5)
m %>% addTiles()

Modificaciones

• addProviderTiles()

m %>% addProviderTiles(providers$Stamen.Toner)

m %>% addProviderTiles(providers$CartoDB.Positron)

m %>% addProviderTiles(providers$Esri.NatGeoWorldMap)

… seguir aquí

3. Paquete maps

Enlace a página

if (!"maps" %in% installed.packages()) install.packages("maps")
library(maps)

Mapa del mundo

map("world")

Otros aspectos a controlar:

• interior = FALSE
• boundary = FALSE
• fill = FALSE
• col = "salmon"
• add = TRUE. Add the interior in the desired colour. As the ‘add’ argument is TRUE add this command to the current plot. If FALSE a new plot is begun.

Se pueden sobrescribir los mapas indicando filly col distintos.

map('world', interior = FALSE, fill = TRUE, col = '#B4045F')
map('world', interior = TRUE, fill= FALSE, col = 'white', add = TRUE)

Mapa de España

world <- map("world", plot = FALSE, namesonly = T)  # Assign the world's regions' names to a new object we called 'world'
canaryIslands <- world[grep("[Cc]anary.*", world)]  # Search for all the regions containing the 'canary' word
canaryIslands

## [1] "Canary Islands:Hierro"        "Canary Islands:Gran Canaria" 
## [3] "Canary Islands:Gomera"        "Canary Islands:Tenerife"     
## [5] "Canary Islands:Fuerteventura" "Canary Islands:La Palma"     
## [7] "Canary Islands:Lanzarote"

Llamada al mapa

map("world", fill = TRUE, col = "white", lty = 0, 
    bg = "lightblue", xlim = c(-15, 8), ylim = c(35, 45))
map("world", fill = FALSE, interior = TRUE, col = "grey", add = TRUE)

Otros mapas

head(map("world", fill = TRUE, col = "white", namesonly = TRUE, plot = FALSE), 30)

##  [1] "Aruba"                            
##  [2] "Afghanistan"                      
##  [3] "Angola"                           
##  [4] "Angola:Cabinda"                   
##  [5] "Anguilla"                         
##  [6] "Albania"                          
##  [7] "Finland:Aland Islands:Foglo"      
##  [8] "Finland:Aland Islands:Eckero"     
##  [9] "Finland:Aland Islands:Fasta Aland"
## [10] "Andorra"                          
## [11] "United Arab Emirates:Abu al Abyad"
## [12] "United Arab Emirates:Sir Baniyas" 
## [13] "United Arab Emirates:Marawwah"    
## [14] "United Arab Emirates:Al Aryan"    
## [15] "United Arab Emirates"             
## [16] "United Arab Emirates:6"           
## [17] "Argentina:Isla de los Estados"    
## [18] "Argentina:Tierra del Fuego"       
## [19] "Argentina:Isla Trinidad"          
## [20] "Argentina"                        
## [21] "Armenia:Artsvashen"               
## [22] "Armenia"                          
## [23] "Armenia:Enclave border"           
## [24] "American Samoa:Tutuila"           
## [25] "Antarctica:1"                     
## [26] "Antarctica:2"                     
## [27] "Antarctica:3"                     
## [28] "Antarctica:4"                     
## [29] "Antarctica:Henry Ice Rise"        
## [30] "Antarctica:6"

4. Paquete sp

Este es un paquete antiguo, aunque facilita algunas cosas, pore ejemplo a la hora de trabajar solo con ciertas comunidades autónomoas, provincias o municipios.

Enlace a página

if (!"sp" %in% installed.packages()) install.packages("sp")
library("sp")
library("rgdal")

Descarga de ficheros

Los ficheros de mapas en este caso se pueden encontrar en gadm.org

ESPAÑA <- read_rds("gadm36_ESP_4_sp.rds") 
# str(ESPAÑA)

En este caso vamos a centrarnos sólo en Andalucía, pafra probar, pero también porque el ficheros sp es grande.

andalucía <- ESPAÑA[ESPAÑA$NAME_1=="Andalucía",]

Los datos recibidos de GADM tienen el texto en una codificación que no es compatible con sp (muestra los acentos y las eñes de manera incorrecta), por lo que tenemos que transformarlos a un formato compatible con la función iconv. La variable NAME_1 contiene los nombres de las Comunidades Autónomas:

# andalucía@data

# andalucía@data$NAME_1 <- as.factor(iconv(as.character(andalucía@data$NAME_1), "latin1", "utf-8"))
# andalucía@data$NAME_2 <- as.factor(iconv(as.character(andalucía@data$NAME_2), "latin1", "utf-8"))
# andalucía@data$NAME_3 <- as.factor(iconv(as.character(andalucía@data$NAME_3), "latin1", "utf-8"))
# andalucía@data$NAME_4 <- as.factor(iconv(as.character(andalucía@data$NAME_4), "latin1", "utf-8"))

Mapa

spplot(andalucía, "NAME_4")

Añadir datos

prueba <-  tibble(NAME_4 = andalucía@data$NAME_4, capital = 1)

datos <- andalucía@data %>% 
  full_join(prueba)

Al fusionar los datos se mezclan las filas, por lo que hay que reordenarlas:

correct.ordering <- match(andalucía@data$NAME_4, datos$NAME_4)
andalucía@data <- datos[correct.ordering,]

Visualización

spplot(andalucía, "capital")

Elementos

Obtenemos las coordenadas

coordinates(andalucía)

Mapas de google

if (!"ggmap" %in% installed.packages()) install.packages("ggmap")
library(ggmap)

if (!"googleVis" %in% installed.packages()) install.packages("googleVis")
library(googleVis)

5. Paquete Siane

Aquí se explica como obtener los mapas: enlace a página

Descarga de datos

# install.packages("devtools")
library(devtools)
# install_github("rOpenSpain/Siane")
library(Siane)

Indicamos el path a la carpeta y lo seleccionamos.

obj <- register_siane("/Users/albtorval/Documents/Aprendizaje/Trabajos/Mapas/")
shp <- siane_map(obj = obj, level = "Municipios")

## Using default year as the latest year 
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/albtorval/Documents/Aprendizaje/Trabajos/Mapas/SIANE_CARTO_BASE_S_3M/anual/20190101", layer: "se89_3_admin_muni_a_x"
## with 8124 features
## It has 19 fields

Visualización

raster::plot(shp)

Incorporar datos del INE

df <- read_csv2("munCADIZ.csv")

df <- df %>% 
  separate(Municipios, into = c("Codigo","Mun"), sep = "\\s", extra = "merge") %>% 
  filter(Codigo != 11)
  # mutate(Codigo = replace(Codigo, Mun == "Cádiz", 11000) %>% as.numeric()) 

#Representamos datos totales para periodo 2019 
df <- df %>% 
  filter(Sexo == "Total" & Periodo == 2019)

by <- "Codigo"
value <- "Total"
level <- "Municipios"

shp_merged <- siane_merge(shp = shp, df = df, by = by, level = level, value = value)

library(RColorBrewer)
pallete_colour <- "OrRd" # Scale of oranges and reds
n <- 5
values_ine <- shp_merged@data[[value]] # Values we want to plot are stored in the shape@data data frame
colors <- brewer.pal(n, pallete_colour) # A pallete from RColorBrewer
style <- "quantile"
library(classInt)
brks <- classIntervals(values_ine, n = n, style = style)
my_pallete <- brks$brks # my_pallete is a vector of breaks
col <- colors[findInterval(values_ine, my_pallete,all.inside=TRUE)]


raster::plot(shp_merged,col=col)
title_plot <- "Población total por municipios en La Rioja"
title(main = title_plot)
library(maptools)
legend(legend = leglabs(round(my_pallete)), fill = colors, x = "bottomright")

6. Paquete plotly

• Enlace 1.
• Enlace 2
• Enlace 3
• También en el siguiente enlace se pueden encontrar distintos ejemplos muy útiles.

La sintáxis de este paquete será muy parecida a la del paquete leaftlet, por esto no se va a indicar, ver en los ejemplos.

Para mejorar el mapa haremos uso de los archivos shapefile ya vistos y otro paquete también para visualización de mapas sf. La instalación del paquete sf puede dar errores al cargar el paquete rgdal, para arreglarlo hacer brew install rgdalen un terminal.

# PAQUETE PLOTLY:
# devtools::install_github('ropensci/plotly')
# install.packages("plotly")
library(plotly)

# PAQUETE RGDAL:
# install.packages("rgdal")
# rgdal::getGDALVersionInfo(str = "--version")

# PAQUETE SF:
#devtools::install_github("r-spatial/sf")
#install.packages("sf")
library(sf)

Ejemplo covid-2019 por Comunidades Autónomas.

Tnedremos que modificar las bases de datos que tenemos.

• Para hacer web scrapping y conseguir los datos de long y lat de municipios.

library(readr)
df_covid <- read_csv("https://raw.githubusercontent.com/datadista/datasets/master/COVID%2019/ccaa_covid19_fallecidos.csv",
            col_types = cols(cod_ine = col_factor()))
df_covid <- df_covid %>% 
  select(cod_ine,CCAA,tail(names(.),1)) %>% 
  filter(CCAA != "Total")

library(readxl)
geo_municipios <- read_xls("listado-geo-municipios-espana.xls",skip=2)
geo_capitales <- read_xls("listado-capitales-CCAA.xls")
df_geo <- geo_capitales %>%
  left_join(geo_municipios, by = c("Ciudad"="Población")) %>% 
  mutate(cod_ine = as.factor(CA)) %>% 
  filter(Capital == "Sí") %>%
  select(cod_ine,Ciudad,Autonomía,Superficie,Habitantes,Latitud,Longitud)

df <- df_covid %>% 
  left_join(df_geo, by="cod_ine") %>% 
  rename(Muertes = 3)

fig <- plot_mapbox() %>% 
  add_markers(data = df, x = ~Longitud, y = ~Latitud, size = ~Muertes,
              color = "red", hoverinfo = "text",
              text = ~paste(df$CCAA, "\n", df$Muertes, " muertos")) 

fig <- fig %>% layout(title = 'Muertes Covid-2019 en España',
                      mapbox = list(style = "carto-darkmatter", zoom = 3,
                                    center = list(lon = -4.905, lat = 38.603))) 

fig 

# Modificaciones a tener en cuenta:  
#   corregir, hacer los habitantes totales de la CCAA
#   plot_geo(spainCCAA) | plot_ly(spainCCAA) | plot_mapbox(spainCCAA)
#   df$Muertes | round(df$Muertes/df$Habitantes,3)
#   style = "carto-darkmatter") # "carto-positron", "carto-darkmatter"

Con sintaxis ggplot2

spainCCAA <- sf::st_read("Comunidades_Autonomas_ETRS89_30N.shp", quiet = TRUE)

spainCCAA <- spainCCAA %>% 
  left_join(df, by = c("Codigo"="cod_ine")) %>% 
  mutate(Latitud  = round(Latitud,2),
         Longitud = round(Longitud,2))

mapa <- spainCCAA
# Puede que haya que modificar la proyección empleada
mapa <- st_transform(mapa, crs = "+proj=longlat +datum=WGS84")

p <- mapa %>% ggplot() +
  geom_sf(aes(fill = Muertes)) + 
  geom_text(aes(Longitud, Latitud, label = CCAA),size=3, colour = "black" ,position = "jitter") +
  scale_fill_gradient(low = "cornflowerblue", high = "blue", na.value = NA) +
  geom_point(aes(x = Longitud, y = Latitud,  size = Muertes), color= "red", alpha = 1/2) +
  scale_size_continuous(range = c(1,10)) +
  labs(title = "Muertes covid-19 en España", 
       subtitle = "Alberto Torrejón Valenzuela") + 
  theme_void() 

p

ggplotly(p) %>%  
  layout(title = 'Muertes Covid-2019 en España')

7. Paquete highmaps

Visto en clase de Inteligencia Artificial.

library(highcharter)

regiones <- read_excel("data/regiones.xlsx")
ccaacsv <- read_csv2("data/ccaa.csv", locale = locale(encoding="ISO-8859-1"),
                     col_types = cols(Codigo = col_character()))
mydata <- fromJSON("https://code.highcharts.com/mapdata/countries/es/es-all.geo.json", 
                   flatten = TRUE)
provincias <- tibble(clave = mydata$features$`properties.hc-key`)
provincias %<>% mutate(valor = sample.int(1000, 53))
(mes <- hcmap("countries/es/es-all", data = provincias,
      name = "clave", value = "valor", joinBy = c("hc-key", "clave"),
      borderColor = "transparent") %>%
  hc_colorAxis(dataClasses = color_classes(seq(0, 1000, by = 100))) %>%
  hc_legend(layout = "vertical", align = "right",
            floating = TRUE, valueDecimals = 0))

Myhcmap <- jsonlite::fromJSON("https://code.highcharts.com/mapdata/countries/es/es-all.geo.json",
                              simplifyVector = FALSE)

highchart() %>%
hc_add_series_map(Myhcmap, provincias, value = "valor", joinBy = c("hc-key", "clave"))

Myhcmap2 <- fromJSON("shapefiles_espana_ccaa_2.geojson", simplifyVector = FALSE)

ccaacsv <- read_csv2("data/ccaa.csv", locale = locale(encoding="ISO-8859-1"),
                     col_types = cols(Codigo = col_character()))
                     #col_factor(ordered = TRUE)))
ccaacsv %<>% mutate(valor = sample.int(1000, 19))

highchart() %>%
hc_add_series_map(Myhcmap2, ccaacsv, value = "valor", 
                  joinBy = c("codaut_sin", "Codigo"),
                  name="Fallecidos",
                    dataLabels = list(enabled = TRUE,
                                      format = '{point.properties.name}')) %>%
  hc_colorAxis(dataClasses = color_classes(seq(0, 1000, by = 100))) %>%
  hc_legend(layout = "vertical", align = "right",
            floating = TRUE, valueDecimals = 0)

8. Otros paquetes

• library(tmap)

Si conocéis cualquier otro paquete que sirva para trabajar con mapas dejad un comentario.