R que R

Maps: Mapas regionales con el paquete {raster}

Thu, Nov 28, 2019
R Mapas
#Raster #maps #mapas


Introducción


En post anteriores hemos abordado diversos aspectos que pueden ser de utilidad a la hora de realizar mapas en R. Así, por ejemplo, hemos explicado cómo trabajar con las nomenclaturas de las unidades territoriales estadísticas (NUTS) de Eurostat y expusimos cómo realizar mapas 3D con el paquete {echarts4r} y e_geo_3d(). Previamente dedicamos un post para explicar cómo realizar mapas de España a diversas unidades territoriales utilizando archivos shapefile y otro para mostrar cómo realizar mapas en R utilizando coordenadas geográficas así como mapas animados utilizando {gganimate}. Por consiguiente, la presente entrada no tiene como objetivo específico representar en mapas alguna variable o indicador determinado, aspecto que ya se ha abarcado en los post mencionados, sino mostrar cómo podemos obtener mapas base de cualquier país a distintos niveles territoriales o geográficos. Para ello resulta de gran utilidad el paquete {raster}, a partir del cual podemos acceder a información geográfica de cualquier país del mundo y, gracias a ello, obtener cualquier mapa base que nos pueda interesar.

Paquetes



library(raster)
library(DT)
library(ggplot2)
library(tidyverse)


la función codes() nos ofrece la información geográfica relevante de 256 unidades territoriales como muestra la tabla siguiente.



datatable(ccodes())


Aunque los encontramos en la tabla anterior, también podemos identificar los código ISO3 de cada uno de los países con la función getData(“ISO3”):



datatable(getData("ISO3"))


Utilizando dichos códigos, con la información que provee gadm.org, podemos representar mapas a diferentes niveles de agregación regional. Veámos algunos ejemplos en los siguientes apartados.


España


Nivel: País


ESP_0 <- getData("GADM", country= "ESP", level=0)
ESP_0
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : -18.16153, 4.328195, 27.63736, 43.79153  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 2
## names       : GID_0, NAME_0 
## value       :   ESP,  Spain



ggplot(ESP_0, aes(long, lat, group=group))+
  labs(title = "País",
       subtitle = "") +
  geom_polygon(fill = "orange", 
               color= "white",
               size = 0.2,
               alpha = 0.8) +
  theme_bw()


Nivel: Comunidades Autónomas


ESP_1 <- getData("GADM", country= "ESP", level=1)
ESP_1
## class       : SpatialPolygonsDataFrame 
## features    : 18 
## extent      : -18.16153, 4.328195, 27.63736, 43.79153  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 10
## names       : GID_0, NAME_0,   GID_1,           NAME_1,                                                                                      VARNAME_1, NL_NAME_1,             TYPE_1,            ENGTYPE_1, CC_1, HASC_1 
## min values  :   ESP,  Spain, ESP.1_1,        Andalucía,                                                      Andalousie|Andaluc¡a|Andalusien|Andaluzia,        NA, Ciudades Autónomas,      Autonomous City,   01,  ES.AN 
## max values  :   ESP,  Spain, ESP.9_1, Región de Murcia, Valencia|Communauté de Valence|Comunidade Valenciana|Comunidad Valenciana|Comunitat Valenciana,        NA, Comunidad Autónoma, Autonomous Community,   19,  ES.VC



ggplot(ESP_1, aes(long, lat, group=group))+
  labs(title = "Comunidades Autónomas",
       subtitle = "") +
  geom_polygon(fill = "orange", 
               color= "white",
               size = 0.2,
               alpha = 0.8) +
  theme_bw()


Nivel: Provincias


ESP_2 <- getData("GADM", country= "ESP", level=2)
ESP_2
## class       : SpatialPolygonsDataFrame 
## features    : 52 
## extent      : -18.16153, 4.328195, 27.63736, 43.79153  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 13
## names       : GID_0, NAME_0,   GID_1,           NAME_1, NL_NAME_1,     GID_2,   NAME_2, VARNAME_2, NL_NAME_2,          TYPE_2,       ENGTYPE_2, CC_2,   HASC_2 
## min values  :   ESP,  Spain, ESP.1_1,        Andalucía,        NA, ESP.1.1_1, A Coruña,   Alacant,        NA, Ciudad Autónoma, Autonomous City,   01, ES.AN.AM 
## max values  :   ESP,  Spain, ESP.9_1, Región de Murcia,        NA, ESP.9.1_1, Zaragoza,  València,        NA,       Provincia,        Province,   52, ES.VC.VN



ggplot(ESP_2, aes(long, lat, group=group))+
  labs(title = "Provincias",
       subtitle = "") +
  geom_polygon(fill = "orange", 
               color= "white",
               size = 0.2,
               alpha = 0.8) +
  theme_bw()


Nivel: Comarcas


ESP_3 <- getData("GADM", country= "ESP", level=3)
ESP_3
## class       : SpatialPolygonsDataFrame 
## features    : 369 
## extent      : -18.16153, 4.328195, 27.63736, 43.79153  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 16
## names       : GID_0, NAME_0,   GID_1,           NAME_1, NL_NAME_1,     GID_2,   NAME_2, NL_NAME_2,       GID_3,         NAME_3, VARNAME_3, NL_NAME_3,  TYPE_3, ENGTYPE_3, CC_3, ... 
## min values  :   ESP,  Spain, ESP.1_1,        Andalucía,        NA, ESP.1.1_1, A Coruña,        NA, ESP.1.1.1_1,       Alacanti,        NA,        NA, Comarca,   Comarca,   NA, ... 
## max values  :   ESP,  Spain, ESP.9_1, Región de Murcia,        NA, ESP.9.1_1, Zaragoza,        NA, ESP.9.1.9_1, Vinalopo Mitja,        NA,        NA, Comarca,   Comarca,   NA, ...



ggplot(ESP_3, aes(long, lat, group=group))+
  labs(title = "Comarcas",
       subtitle = "") +
  geom_polygon(fill = "orange", 
               color= "white",
               size = 0.2,
               alpha = 0.8) +
  theme_bw()


Nivel: Municipios


ESP_4 <- getData("GADM", country= "ESP", level=4)
ESP_4
## class       : SpatialPolygonsDataFrame 
## features    : 8302 
## extent      : -18.16153, 4.328195, 27.63736, 43.79153  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 14
## names       : GID_0, NAME_0,   GID_1,           NAME_1,     GID_2,   NAME_2,       GID_3,         NAME_3,         GID_4,   NAME_4, VARNAME_4,              TYPE_4,           ENGTYPE_4, CC_4 
## min values  :   ESP,  Spain, ESP.1_1,        Andalucía, ESP.1.1_1, A Coruña, ESP.1.1.1_1,       Alacanti, ESP.1.1.1.1_1, A Arnoia,  Abadiano,        Municipality,        Municipality,   NA 
## max values  :   ESP,  Spain, ESP.9_1, Región de Murcia, ESP.9.1_1, Zaragoza, ESP.9.1.9_1, Vinalopo Mitja, ESP.9.1.9.9_1,  Zurgena,      Zuya, Plazas de soberanía, Sovereign territory,   NA



ggplot(ESP_4, aes(long, lat, group=group))+
  labs(title = "Municipios",
       subtitle = "") +
  geom_polygon(fill = "orange", 
               color= "white",
               size = 0.2,
               alpha = 0.8) +
  theme_bw()


Perú


Nivel: País


PER_0 <- getData("GADM", country= "PER", level=0)
PER_0
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : -81.3307, -68.65311, -18.3518, -0.03747  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 2
## names       : GID_0, NAME_0 
## value       :   PER,   Peru



ggplot() +
   geom_polygon(data= PER_0,
                aes( x= long, y=lat, group= group),
                color= "white",
                fill = "palegreen3", 
                size = 0.2,
               alpha = 0.8) +
  labs(title = "País",
       subtitle = "") +
  theme_bw()

Nivel: Regiones / Departamentos


En esta ocasión vamos a realizar un mapa a nivel regional/departamental (level=1), al que añadiremos los nombres de cada una de las regiones.



PER_1 <- getData("GADM", country= "PER", level=1)
PER_1
## class       : SpatialPolygonsDataFrame 
## features    : 26 
## extent      : -81.3307, -68.65311, -18.3518, -0.03747  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 10
## names       : GID_0, NAME_0,   GID_1,   NAME_1, VARNAME_1, NL_NAME_1,    TYPE_1, ENGTYPE_1, CC_1, HASC_1 
## min values  :   PER,   Peru, PER.1_1, Amazonas,   Ancachs,        NA, Provincia,  Province,   01,  PE.AM 
## max values  :   PER,   Peru, PER.9_1,  Ucayali,    Tumbez,        NA,    Región,    Region,   25,  PE.UC


Los nombres de las regiones se encuentran en una columna con nombre “NAME_1”. Realizamos una serie de modificaciones tras las cuales obtenemos un dataset con las regiones y unas coordenadas geográficas que utilizaremos en el mapa posterior.



PERU_df <- broom::tidy(PER_1, region = "NAME_1")
lapply(PERU_df, class)
## $long
## [1] "numeric"
## 
## $lat
## [1] "numeric"
## 
## $order
## [1] "integer"
## 
## $hole
## [1] "logical"
## 
## $piece
## [1] "factor"
## 
## $group
## [1] "factor"
## 
## $id
## [1] "character"

head(PERU_df)
## # A tibble: 6 x 7
##    long   lat order hole  piece group      id      
##   <dbl> <dbl> <int> <lgl> <fct> <fct>      <chr>   
## 1 -78.0 -6.97     1 FALSE 1     Amazonas.1 Amazonas
## 2 -78.0 -6.97     2 FALSE 1     Amazonas.1 Amazonas
## 3 -78.0 -6.96     3 FALSE 1     Amazonas.1 Amazonas
## 4 -78.0 -6.96     4 FALSE 1     Amazonas.1 Amazonas
## 5 -78.0 -6.96     5 FALSE 1     Amazonas.1 Amazonas
## 6 -78.0 -6.96     6 FALSE 1     Amazonas.1 Amazonas

region_names <- aggregate(cbind(long, lat) ~id, data = PERU_df, FUN = mean)

head(region_names)
##          id      long        lat
## 1  Amazonas -78.18920  -4.577463
## 2    Ancash -78.14134  -9.579794
## 3  Apurímac -73.07987 -13.990960
## 4  Arequipa -73.06728 -16.170506
## 5  Ayacucho -74.02473 -13.923152
## 6 Cajamarca -78.80095  -6.229225



ggplot() +
  geom_polygon(data = PER_1,
               aes(x = long, y = lat, group=group),
               colour = "white",
               fill = "palegreen3")+
  labs(title = "Regiones / Departamentos",
       subtitle = "")+ 
  geom_text(data=region_names, 
            aes(x= long, y=lat, label =id), 
            size = 2, color = "gray20") + 
  theme_bw()


Provincias

 


PER_2 <- getData("GADM", country= "PER", level=2)
PER_2
## class       : SpatialPolygonsDataFrame 
## features    : 195 
## extent      : -81.3307, -68.65311, -18.3518, -0.03747  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 13
## names       : GID_0, NAME_0,   GID_1,   NAME_1, NL_NAME_1,     GID_2,    NAME_2,            VARNAME_2, NL_NAME_2,     TYPE_2,  ENGTYPE_2, CC_2,   HASC_2 
## min values  :   PER,   Peru, PER.1_1, Amazonas,        NA, PER.1.1_1,   Abancay,             Asuncitn,        NA,  Provincia,   Province,   NA, PE.AM.BG 
## max values  :   PER,   Peru, PER.9_1,  Ucayali,        NA, PER.9.7_1, Zarumilla, Rodreguez de Mendoza,        NA, Water body, Water body,   NA, PE.UC.PU



ggplot(PER_2, aes(long, lat, group=group))+
  labs(title = "Provincias",
       subtitle = "") +
  geom_polygon(fill = "palegreen3", 
               color= "white",
               size = 0.1,
               alpha = 0.8) +
  theme_bw()


En esta ocasión vamos a identificar grupos de regiones. Si observamos el tipo de provincias peruanas comprobamos que una de las unidades geográficas que incluye el dataset en lugar de ser provincia es una “Water body”. Sin duda dicha unidad geográfica debe corresponder al Lago Titicaca. Consecuentemente vamos a realizar un mapa donde distingamos el Lago Titicaca del resto del país. Para ello debemos realizar una serie de modificaciones como veremos a continuación.



PER_2@data$TYPE_2
##   [1] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##   [6] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [11] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [16] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [21] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [26] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [31] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [36] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [41] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [46] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [51] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [56] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [61] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [66] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [71] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [76] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [81] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [86] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [91] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
##  [96] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [101] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [106] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [111] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [116] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [121] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [126] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [131] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [136] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [141] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [146] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [151] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [156] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [161] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [166] "Water body" "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [171] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [176] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [181] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [186] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia" 
## [191] "Provincia"  "Provincia"  "Provincia"  "Provincia"  "Provincia"



PERU_2_df <- broom::tidy(PER_2, region = "NAME_2")
lapply(PERU_2_df, class)
## $long
## [1] "numeric"
## 
## $lat
## [1] "numeric"
## 
## $order
## [1] "integer"
## 
## $hole
## [1] "logical"
## 
## $piece
## [1] "factor"
## 
## $group
## [1] "factor"
## 
## $id
## [1] "character"

PER_2@data <- PER_2@data %>% mutate(isProvince = as.factor(ifelse(TYPE_2 == "Provincia", yes = "Provincia", no = "Water body")))

PERU_2_df <- broom::tidy(PER_2, region = "isProvince")



ggplot() + 
  geom_polygon(data = PERU_2_df, 
               aes(x = long, y = lat, 
                   group = group, 
                   fill = id), 
               colour = "white") +
  labs(title = "Provincias",
       subtitle = "") +
  theme_bw() +
  scale_fill_manual(values = c("palegreen3",  "deepskyblue"))


Distritos



PER_3 <- getData("GADM", country= "PER", level=3)
PER_3
## class       : SpatialPolygonsDataFrame 
## features    : 1815 
## extent      : -81.3307, -68.65311, -18.3518, -0.03747  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 16
## names       : GID_0, NAME_0,   GID_1,   NAME_1, NL_NAME_1,     GID_2,    NAME_2, NL_NAME_2,       GID_3,         NAME_3, VARNAME_3, NL_NAME_3,     TYPE_3,  ENGTYPE_3, CC_3, ... 
## min values  :   PER,   Peru, PER.1_1, Amazonas,        NA, PER.1.1_1,   Abancay,        NA, PER.1.1.1_1, 3 de Diciembre,        NA,        NA,   Distrito,   District,   NA, ... 
## max values  :   PER,   Peru, PER.9_1,  Ucayali,        NA, PER.9.7_1, Zarumilla,        NA, PER.9.7.9_1,         Zurite,        NA,        NA, Water body, Water body,   NA, ...



ggplot(PER_3, aes(long, lat, group=group))+
  labs(title = "Distritos",
       subtitle = "") +
  geom_polygon(fill = "palegreen3", 
               color= "white",
               size = 0.1,
               alpha = 0.8) +
  theme_bw()