Si alguna vez han sufrido como yo, con mapas que exceden los recursos de su ordenador (multiples capas que se sobre-empalman, tamaño de la geometría inoperable), entonces el paquete rmapshaper es para ustedes. A continuación comparto algunos pasos para editar un shapefile que tiene más de 3 mil capas, las cuales muchas de ellas se sobreponen. Aunado a ello, el archivo pesa más de 260 megas, por lo que mi limitado ordenador apenas si me deja generar mapas en un tiempo razonable. Aqui la solución:

Configuración general

Llamamos a toda la libreria que usaremos

library(sf) # lectura de shp files
## Warning: package 'sf' was built under R version 3.6.2
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(tidyverse) # Manipulacion de datos
## ── Attaching packages ──────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0           ✓ purrr   0.3.3      
## ✓ tibble  3.0.0           ✓ dplyr   0.8.99.9002
## ✓ tidyr   1.0.2           ✓ stringr 1.4.0      
## ✓ readr   1.3.1           ✓ forcats 0.5.0
## Warning: package 'tibble' was built under R version 3.6.2
## ── Conflicts ─────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(dplyr) # En caso de querer usar dplyr 1.0
library(pacman) # Instala y carga librerias no existentes

Revisión general del archivo SHP

Llamamos y revisamos el shp

shp_original<-st_read("territorios_for_whg_NEW.shp")
## Reading layer `territorios_for_whg_NEW' from data source `/Users/dante/MEGA/JMP/Historia/Bases/Mapas/Colonia/territorios_for_whg_NEW.shp' using driver `ESRI Shapefile'
## Simple feature collection with 3054 features and 12 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -126.9816 ymin: -52.34173 xmax: -48.32736 ymax: 49.88676
## CRS:            4326
head(shp_original) # Variables contenidas
## Simple feature collection with 6 features and 12 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -101.13 ymin: -13.39983 xmax: -71.99898 ymax: 23.09539
## CRS:            4326
##         Tipo       Label                Nombre
## 1 Marquesado     Oropesa               Oropesa
## 2    Condado Casa Bayona           Casa Bayona
## 3     Ducado     Atrisco               Atrisco
## 4 Marquesado     Bejucal San Felipe y Santiago
## 5 Marquesado    El Valle              El Valle
## 6 Marquesado       Guiza                 Guiza
##                                       Variantes Entidad_ID   Nivel  AbrA  AbrB
## 1                        Marquesado de Oropesa;   SEPEMO00 Senorio Marq. Marq.
## 2  Condado de Casa Bayona; Casa-Bayona| Rosario   SESDRO00 Senorio Cond. Cond.
## 3                            Ducado de Atlixco;   SENEDA00 Senorio  Duc.  Duc.
## 4  Marquesado de San Felipe y Santiago; Bejucal   SESDFS00 Senorio Marq. Marq.
## 5                         Marquesado del Valle;   SENEMV00 Senorio Marq. Marq.
## 6 Marquesado de Guiza; Guisa| San Jose de Guiza   SESDGU00 Senorio Marq. Marq.
##   Niv_Ent_id                                          Wikilink START END_
## 1       <NA>                   per_cuzco_marquesado de oropesa  1701 1740
## 2       <NA>              sdo_la habana_condado de casa bayona  1732 1808
## 3       <NA>                     nes_senorio_ducado de atlixco  1710 1808
## 4       <NA> sdo_la habana_marquesado de san felipe y santiago  1713 1808
## 5       <NA>                  nes_senorio_marquesado del valle  1701 1808
## 6       <NA>                      sdo_cuba_marquesado de guiza  1774 1808
##                         geometry
## 1 MULTIPOLYGON (((-72.1507 -1...
## 2 MULTIPOLYGON (((-82.25132 2...
## 3 MULTIPOLYGON (((-96.52511 1...
## 4 MULTIPOLYGON (((-82.36249 2...
## 5 MULTIPOLYGON (((-95.32976 1...
## 6 MULTIPOLYGON (((-76.53549 2...
class(shp_original) # Vemos que tipo de objeto es (sf y df)
## [1] "sf"         "data.frame"
str(shp_original) # Características de las variables del df
## Classes 'sf' and 'data.frame':   3054 obs. of  13 variables:
##  $ Tipo      : Factor w/ 59 levels "Alcaldia mayor",..: 37 16 23 37 37 37 37 31 31 27 ...
##  $ Label     : Factor w/ 803 levels "Abancay","Acaponeta",..: 476 119 52 77 232 270 431 550 785 344 ...
##  $ Nombre    : Factor w/ 873 levels "Abancay","Acaponeta",..: 508 119 55 625 249 287 463 586 855 369 ...
##  $ Variantes : Factor w/ 823 levels "Abancay; Anta",..: 395 189 228 396 397 393 394 78 807 342 ...
##  $ Entidad_ID: Factor w/ 892 levels "AUBA0000","AUCA0000",..: 882 885 880 883 881 884 886 806 776 816 ...
##  $ Nivel     : Factor w/ 6 levels "Audiencia","Jurisdiccion",..: 6 6 6 6 6 6 6 5 5 5 ...
##  $ AbrA      : Factor w/ 25 levels "A.","AM","Ao.",..: 16 10 12 16 16 16 16 14 14 14 ...
##  $ AbrB      : Factor w/ 24 levels "AM","Arzob.",..: 15 9 10 15 15 15 15 5 5 19 ...
##  $ Niv_Ent_id: Factor w/ 14 levels "AUCC0000","AUPA0000",..: NA NA NA NA NA NA NA NA NA 9 ...
##  $ Wikilink  : Factor w/ 886 levels "cha_audiencia_charcas",..: 613 780 489 784 490 772 783 752 452 793 ...
##  $ START     : num  1701 1732 1710 1713 1701 ...
##  $ END_      : num  1740 1808 1808 1808 1808 ...
##  $ geometry  :sfc_MULTIPOLYGON of length 3054; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:149, 1:2] -72.2 -72.1 -72.1 -72.1 -72.1 ...
##   ..- attr(*, "class")= chr  "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr  "Tipo" "Label" "Nombre" "Variantes" ...
nrow(shp_original) # Número de capas
## [1] 3054
st_crs(shp_original) # CRS de la geometría
## Coordinate Reference System:
##   User input: 4326 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]

Removamos la geometría para quedarnos solo con el data.frame:

shp_sin_geometria <- shp_original  %>%  st_set_geometry(NULL)

Transformación y selección de datos y capas

Identificación de objetos de interés dentro de nuestro shapefile

Lo primero que debemos hacer es identificar las capas que son de nuestro interés en nuestro análisis. (los lugares que pueden ser parte del actual México en 1804).

lista_lugares1<-list("nes","gdj","sdo")
lista_lugares2<-list("gua_ch")

Filtramos por año para reducir el numero de capas

shp_modificado<-shp_sin_geometria%>%filter(START==1804)

Ahora incluimos solo las capas con los lugares definidos previamente

shp_modificado<-shp_modificado%>%
  select(Wikilink,START)%>%
  filter(substr(Wikilink,1,3) %in% lista_lugares1 | substr(Wikilink,1,6)%in% lista_lugares2)%>%
  unique()

Filtramos y unimos el shp con el df sin geometria para reducir el numero de capas. En este caso nos interesan las Provincias, Obispados, y Jurisdicciones definidas en 1804.

shp_modificado<-shp_original%>%inner_join(shp_modificado) %>% filter(Nivel=="Provincia" |Nivel=="Jurisdiccion"|Nivel=="Obispado" )
## Joining, by = c("Wikilink", "START")

Re-escalamiento del objeto sf

Primero vamos a instalar y cargar rmapshaper que va a disminuir el tamaño de la geometria

p_load(rmapshaper)

Veamos ahora cuanto se tarda R en mapear nuestro objeto original

system.time(ggplot (data=shp_modificado) + geom_sf(aes(fill=NULL)))
##    user  system elapsed 
##   0.008   0.000   0.009
ggplot (data=shp_modificado) + geom_sf(aes(fill=NULL))

Ahora vamos a reducir el tamaño del objeto con rmapshaper con sus opciones definidas por default. Si tienes curiosidad checa el detalle en:https://github.com/ateucher/rmapshaper

shp_GIS_reducido <- ms_simplify(shp_modificado)
system.time(ggplot (data=shp_GIS_reducido) + geom_sf(aes(fill=NULL)))
##    user  system elapsed 
##   0.003   0.000   0.005
ggplot (data=shp_GIS_reducido) + geom_sf(aes(fill=NULL))

Como puedes ver disminuyo el tiempo de la operacion ya que la geometria esta mas “simplificada”. A continuación seamos más agresivos en el cambio:

shp_GIS_reducido <- ms_simplify(shp_modificado, keep = 0.01)
system.time(ggplot (data=shp_GIS_reducido) + geom_sf(aes(fill=NULL)))
##    user  system elapsed 
##   0.001   0.000   0.001
ggplot (data=shp_GIS_reducido) + geom_sf(aes(fill=NULL))

Como ultimo ejercicio, creemos una variable y peguemosla en nuestro objeto:

IDMexico<-data.frame(ID=c(1,0,1,0,0,1,0,1,1,1,1),
Entidad_ID=c("PRGDZA00","PRSDFL00","PRGDGD00","PRSDLH00","PRNECL00","OBNEYU00","JUNECLCV","JUNECLCN","PMGDCOEN",
"PTGDCOMO","JUNEMICL")) %>% mutate(ID=as.numeric(ID))
shp_GIS_reducido<-inner_join(shp_GIS_reducido,IDMexico)
## Joining, by = "Entidad_ID"
ggplot(data = shp_GIS_reducido) +
    geom_sf(aes(fill = ID))

Finalmente exportemos nuestro objeto reducido

st_write(shp_GIS_reducido,"Ejemplo.shp",driver = "ESRI shapefile", overwrite=T)
## Writing layer `Ejemplo' to data source `Ejemplo.shp' using driver `ESRI shapefile'
## Writing 11 features with 13 fields and geometry type Multi Polygon.

Resultado el archivo paso de 260 megas a menos de 1 mega!!! Tamaño{width=75%, center}