Manipulando el mapa de los cantones y distritos de Costa Rica con R

Trabajar con datos geoespaciales anteriores a 2017 a menudo resulta en mapas con espacios vacíos o zonas grises, debido a cambios en fronteras, actualización de códigos o datos ausentes. Este tutorial muestra cómo resolver estos problemas y unir regiones usando el paquete ‘sf’.

Autor/a

Afiliación

Carlos Aguero

Aprende Tidyverse

Hace un par de semanas, en nuestras sesiones de consultas los domingos a las 6 p.m., nos plantearon el siguiente problema: ¿qué pasa si, al descargar el mapa actualizado, como vimos en el post anterior, mis datos no están actualizados? Dicho de otra forma, en mi mapa existen los nuevos cantones como Río Cuarto, Puerto Jiménez y Monteverde, pero no en mis datos, y esas zonas aparecerán en gris. Bueno, en esta guía veremos paso a paso cómo regresar el mapa a la forma que tenía en 2017, uniendo estos nuevos cantones al cantón al que pertenecían, y así poder visualizar nuestro mapa sin regiones sin datos. Esto nos servirá también para crear nuestras propias agrupaciones de áreas, por ejemplo, un mapa de las regiones socioeconómicas de Costa Rica.

Puedes descargar los datos de los cantones y distritos en este enlace

Unir cantones

Paso 1: Cargar nuestros datos

En este caso, asumiremos que ya leíste y ejecutaste lo visto en Actualizando el mapa de cantones de Costa Rica con R, así que importaremos directamente el archivo cantones.geojson.

library(tidyverse)
library(readxl)
library(sf)
library(janitor)

cantones <- read_sf("cantones.geojson") |> 
  clean_names() |> 
  mutate(codigo_canton = as.character(codigo_canton))

glimpse(cantones)

Rows: 84
Columns: 5
$ codigo_de_provincia <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ canton              <chr> "San José", "Escazú", "Desamparados", "Puriscal", …
$ provincia           <chr> "San José", "San José", "San José", "San José", "S…
$ codigo_canton       <chr> "101", "102", "103", "104", "105", "106", "107", "…
$ geometry            <MULTIPOLYGON [°]> MULTIPOLYGON (((-84.15569 9..., MULTI…

Por otro lado, utilizaremos los datos del IDS de Mideplan de 2017, intencionalmente, para no tener datos para nuestros nuevos cantones.

ids <- read_excel(
  path = "IDS 2017.xlsx",
  sheet = "IDS cantonal",
  skip = 2
) |> 
clean_names() |> 
select(codigo_canton = codigo, ids_2017)

glimpse(ids)

Rows: 83
Columns: 2
$ codigo_canton <chr> "102", "409", "407", "115", "408", "406", "118", "109", …
$ ids_2017      <dbl> 100.00000, 99.93229, 98.10882, 96.76766, 96.70938, 91.37…

Paso 2: Unir nuestros datos

En este caso siempre debemos recordar como se menciona en Actualizando el mapa de cantones de Costa Rica con R que debemos tener presente los posibles cambios de nombre de los cantones en el tiempo.

cantones <- left_join(cantones, ids)
glimpse(cantones)

Rows: 84
Columns: 6
$ codigo_de_provincia <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ canton              <chr> "San José", "Escazú", "Desamparados", "Puriscal", …
$ provincia           <chr> "San José", "San José", "San José", "San José", "S…
$ codigo_canton       <chr> "101", "102", "103", "104", "105", "106", "107", "…
$ geometry            <MULTIPOLYGON [°]> MULTIPOLYGON (((-84.15569 9..., MULTI…
$ ids_2017            <dbl> 75.39888, 100.00000, 71.88500, 53.31167, 23.69962,…

Paso 3: Visualizar nuestro mapa

Ahora que tenemos nuestros datos, vamos a crear un mapa donde representaremos en colores el IDS de cada cantón. Además, para identificarlos mejor, pondremos en rojo aquellos cantones que no tengan datos, es decir, su IDS es NA.

ggplot(
  data = cantones,
  mapping = aes(
    fill = ids_2017
  )
) +
  geom_sf(color = "white", linewidth = 0.1) +
  scale_fill_continuous(na.value = "red") +
  theme_minimal() +
  theme(
    legend.position = "bottom"
  )

Podemos ver a simple vista que tenemos 2 zonas en rojo, pero apliquemos un filtro para verlas con más claridad.

cantones |> 
  as.data.frame() |> 
  filter(if_any(everything(), is.na)) |> 
  select(provincia, canton,codigo_canton, ids_2017) 

   provincia         canton codigo_canton ids_2017
1 Puntarenas     Monteverde           612       NA
2 Puntarenas Puerto Jiménez           613       NA

Podemos ver que efectivamente son los cantones de Puerto Jiménez¹ y Monteverde² los que no presentan datos.

Para términos de nuestro ejemplo, vamos a “regresar” Monteverde al cantón de Puntarenas y Puerto Jiménez al cantón de Golfito.

cantones |> 
  as.data.frame() |> 
  filter(canton %in% c("Golfito","Puntarenas")) |> 
  select(canton, codigo_canton)

      canton codigo_canton
1 Puntarenas           601
2    Golfito           607

Paso 3: Unificar nombre y códigos de los cantones a unir

Para esto, debemos asegurarnos de que tanto el nombre del cantón como el código del cantón sean iguales para los cantones que queremos unir.

cantones <- read_sf("cantones.geojson") |> 
  clean_names() |> 
  mutate(codigo_canton = as.character(codigo_canton))


cantones <- cantones |> 
  mutate(
    canton = case_match(
      canton,
      "Puerto Jiménez" ~ "Golfito",
      "Monteverde" ~ "Puntarenas",
      .default = canton
    ),
    codigo_canton = case_match(
      codigo_canton,
      "613" ~ "607",
      "612" ~ "601",
      .default = codigo_canton
    )
  )

Ahora que tenemos los nombres correctos, podemos utilizar la función st_union para unificar esas regiones en un solo polígono.

Paso 4: Unimos los polígonos de las regiones que queremos ‘reunificar’

cantones <- cantones |> 
  group_by(codigo_de_provincia, provincia, codigo_canton, canton) |> 
  summarise(
    geometry = st_union(geometry)
  )

`summarise()` has grouped output by 'codigo_de_provincia', 'provincia',
'codigo_canton'. You can override using the `.groups` argument.

Ahora realizamos el proceso anterior, unimos los datos del IDS y realizamos el gráfico.

cantones <- left_join(cantones, ids)
glimpse(cantones)

Rows: 82
Columns: 6
Groups: codigo_de_provincia, provincia, codigo_canton [82]
$ codigo_de_provincia <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ provincia           <chr> "San José", "San José", "San José", "San José", "S…
$ codigo_canton       <chr> "101", "102", "103", "104", "105", "106", "107", "…
$ canton              <chr> "San José", "Escazú", "Desamparados", "Puriscal", …
$ geometry            <POLYGON [°]> POLYGON ((-84.15569 9.96885..., POLYGON ((…
$ ids_2017            <dbl> 75.39888, 100.00000, 71.88500, 53.31167, 23.69962,…

ggplot(
  data = cantones,
  mapping = aes(
    fill = ids_2017
  )
) +
  geom_sf(color = "white", linewidth = 0.1) +
  scale_fill_continuous(na.value = "red") +
  theme_minimal() +
  theme(
    legend.position = "bottom"
  )

Para verificar que realmente hemos unificado los polígonos, veamos solamente el cantón de Golfito, en el que podemos ver a Puerto Jiménez (la de la izquierda para los que no conocen bien la geografía de la zona sur).

golfo <- cantones |> 
  filter(canton == "Golfito")

ggplot(
  data = golfo,
  mapping = aes(
    fill = ids_2017
  )
) +
  geom_sf(color = "white", linewidth = 0.1) +
  scale_fill_continuous(na.value = "red") +
  theme_minimal() +
  theme(
    legend.position = "bottom"
  )

Unir distritos en regiones de planificación

En este ejemplo, utilizaremos los datos de los distritos y los agruparemos según su región administrativa en Costa Rica, que incluye Brunca, Central, Chorotega, Huetar Caribe, Huetar Norte y Pacífico Central.

La Reforma de la División Regional del Territorio de Costa Rica, Decreto N° 9501, establece las regiones oficiales para la investigación y planificación del desarrollo socioeconómico ³

Puedes descargar los datos de los cantones y distritos en este enlace

distritos <- read_sf("distritos_cr.geojson") |> 
clean_names()

glimpse(distritos)

Rows: 491
Columns: 16
$ objectid         <int> 1131, 1132, 1133, 1134, 1135, 1136, 1137, 1138, 1139,…
$ provincia        <chr> "Limón", "Heredia", "Heredia", "Alajuela", "Heredia",…
$ canton           <chr> "Pococí", "Sarapiquí", "Sarapiquí", "San Carlos", "Sa…
$ region           <chr> "Huetar Caribe", "Huetar Norte", "Huetar Norte", "Hue…
$ oficializacion   <dttm> 1971-07-02, 1970-11-18, 1999-09-16, 1948-11-05, 1999…
$ codigo_canton    <int> 702, 410, 410, 210, 410, 213, 210, 213, 214, 213, 214…
$ global_id        <chr> "{59F87CBA-7CA7-4E12-918D-822CB26FDB38}", "{96EC74A5-…
$ distrito         <chr> "Colorado", "Puerto Viejo", "Llanuras del Gaspar", "P…
$ gml_id           <chr> "limitedistrital_5k.1", "limitedistrital_5k.2", "limi…
$ codigo           <int> 160105, 160105, 160105, 160105, 160105, 160105, 16010…
$ codigo_provincia <int> 7, 4, 4, 2, 4, 2, 2, 2, 2, 2, 2, 5, 5, 5, 2, 6, 5, 6,…
$ codigo_dta       <int> 70206, 41001, 41004, 21006, 41005, 21307, 21013, 2130…
$ legislacion      <chr> "Decreto  1825-G", "Ley 4671", "Decreto  28137-G", "D…
$ area             <dbl> 948.43468, 428.51741, 267.21675, 379.26541, 369.70036…
$ version          <dbl> 20240502001, 20240502001, 20240502001, 20240502001, 2…
$ geometry         <MULTIPOLYGON [°]> MULTIPOLYGON (((-83.41636 1..., MULTIPOL…

Si realizamos un gráfico con los datos completos, veremos que se incluyen los 491 distritos de Costa Rica (a junio de 2024).

ggplot(data = distritos) +
  geom_sf(color = "white", linewidth = 0.1, fill = "steelblue") +
  theme_minimal()

Ahora, podemos tomar este mapa y utilizar la función st_union junto con la variable región para crear un mapa de las regiones administrativas. Esto nos brinda mucha flexibilidad para adaptar los datos disponibles, permitiéndonos realizar las agrupaciones que mejor nos convengan.

regiones <- distritos |> 
  group_by(region) |> 
  summarise(
    geometry = st_union(geometry)
  )

ggplot(
  data = regiones,
  mapping = aes(
    fill = region
  )) +
  geom_sf(color = "white", linewidth = 0.1) +
  theme_minimal()

Espero que esta pequeña guía te sea de utilidad. Si tienes alguna consulta, no dudes en contactarme a carlos.aguero@aprendetidyverse.com o unirte a nuestros Domingos de Preguntas y Respuestas, un espacio público y gratuito donde conversamos y resolvemos juntos problemas con código R. Este espacio está diseñado para brindar acceso a toda la comunidad de usuarios de R, ofreciendo consultoría y apoyo en su proceso de aprendizaje sin costo alguno.

Notas

1. Delfino.cr. (21 de octubre de 2021). Expediente 22749: Creación del Cantón de Puerto Jiménez, Cantón XIII de la Provincia de Puntarenas. Delfino.cr. Leer más

2. Delfino.cr. (25 de septiembre de 2019). Expediente 21618: Creación del Cantón de Monteverde, Cantón XII de la Provincia de Puntarenas. Delfino.cr. Leer más

3. Costa Rica. Decreto Ejecutivo N° 9501. Reforma División Regional del Territorio de Costa Rica, para los Efectos de Investigación y Planificación del Desarrollo Socioeconómico. (11 de febrero de 2015). Procuraduría General de la República. Leer más