Mapeando

Siempre me ha gustado el tema de los Sistemas de información geográfica y derivados. Ya cuando trabajaba en el IESA fui a un curso en Vigo sobre gvSIG y luego aprendí cosas con QGIS , el cual me sigue pareciendo un software más que excelente.

Hoy en día se pueden hacer muchísimas cosas con Python y con R, incluso tienen conectores con librerías de javascript para hacer cosas resultonas y con interactividad, veáse esto por ejemplo. En realidad tanto python como R tiran de las mismas librerías de bajo nivel, a saber, gdal, libproj, libgeos y similares, las cuales, para ser honestos, pueden meterte en un infierno de dependencias en según qué sistemas unix. Eso sí, una vez sales de ahí es una gozada.

Total, el caso es que el pasado jueves me llega una duda del gran Leonardo Hansa, que reproduzco a continuación.

Vicepresi, tú sabes cómo se puede saber si dos códigos postales son adyancentes? Con datos, claro, no a ojo

Y bueno, como no tenía mucho más que hacer esa tarde y soy un picado de la vida, pues me dije.

Creo que tengo un shapefile de códigos postales de 2019 que no recuerdo como llegó a mis manos, quizá pueda hacer algo.

Y me puse manos a la obra a investigar a ver cómo se podría hacer eso de encontrar polígonos adyacentes a uno dado usando Rstats. El caso es que llegué a varios hilos dónde se comentaban diferentes formas, una era usando la función poly2nb de la librería spdep y también cómo hacerlo usando la función st_intersects de la librería sf . En este issue incluso comentaban los grandes de estas cosas, Edzer Pebesma , Roger Bivand y nuestro conocico Virgilio Gómez Rubio.

Bueno, vamos al código, que de eso se trata.

library(tidyverse)
library(sf)

mapa <- st_read(here::here("data/shapefiles/cod_postales/cp_19_dis.shp"))
#> Reading layer `cp_19_dis' from data source 
#>   `/media/hd1/canadasreche@gmail.com/blog_quarto/data/shapefiles/cod_postales/cp_19_dis.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 10809 features and 4 fields (with 1 geometry empty)
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -2021660 ymin: 3203377 xmax: 481778.5 ymax: 5433002
#> Projected CRS: WGS 84 / Pseudo-Mercator

head(mapa)
#> Simple feature collection with 6 features and 4 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -1536953 ymin: 3373964 xmax: -41802.13 ymax: 5247186
#> Projected CRS: WGS 84 / Pseudo-Mercator
#>      cp cp_num cp_2 cp_2_num                       geometry
#> 1 35560  35560   35       35 MULTIPOLYGON (((-1518970 33...
#> 2 27330  27330   27       27 MULTIPOLYGON (((-821864.3 5...
#> 3 46680  46680   46       46 MULTIPOLYGON (((-51610.46 4...
#> 4 49706  49706   49       49 MULTIPOLYGON (((-641488.4 5...
#> 5 21120  21120   21       21 MULTIPOLYGON (((-776955.2 4...
#> 6 16623  16623   16       16 MULTIPOLYGON (((-256256.7 4...

Pinto códigos postales de la provincia de Madrid por ejemplo, coloreando por área de cada polígono.

mapa |>
    
    # calculo area usando st_area
    mutate(area_m2 = st_area(mapa) |> 
               as.numeric()) |>
    filter(cp_2 == 28) |>
    ggplot() +
    geom_sf(aes(fill = area_m2)) +
    scale_fill_viridis_c()

Para encontrar los polígonos adyacentes uso st_intersects que permite por ejemplo saber qué puntos están dentro de un polígono y cosas así. Al aplicarlo sobre una geometría de tipo polígono lo que encuentra son los polígonos adyacentes.

# mido a ver cuánto tarda. 
tictoc::tic()
lista_adyacentes <- st_intersects(mapa)
tictoc::toc()
#> 2.531 sec elapsed

Ahora si quisiera saber qué polígonos son adyacentes a uno dado, es simplemente seleccionar en lista adyacentes, por ejemplo

lista_adyacentes[1]
#> [[1]]
#> [1]    1 3684 4621 6345 7017 7996

Para ver los adyacentes a mi código postal en Madrid.

(fila_mi_cp <-  mapa |> 
    rownames_to_column() |> 
    filter(cp == "28043") |> 
    pull(rowname))
#> [1] "9138"


(mis_vecinos_fila <- lista_adyacentes[as.numeric(fila_mi_cp)])
#> [[1]]
#> [1]  493 4831 4998 5127 6941 7662 9138 9427

Pero me devuelve el número de fila, para ver el cp sería

mi_cp <-  "28043"
(mis_vecinos <- mapa$cp[mis_vecinos_fila[[1]]])
#> [1] "28028" "28042" "28002" "28022" "28033" "28016" "28043" "28027"

(adyacentes <-  setdiff(mis_vecinos, mi_cp))
#> [1] "28028" "28042" "28002" "28022" "28033" "28016" "28027"

Pintamos

  mapa |> 
    filter(cp_2 == 28) |>
        mutate(
            tipo = case_when(
                cp == mi_cp ~ "mi_cp", 
                cp %in% adyacentes ~ "mis_copostales_vecinos", 
                TRUE ~ "resto de codpostales"
            )
        ) |> 
        ggplot() +
        geom_sf(aes(fill = tipo)) +
        scale_fill_viridis_d()

y listo.. Lo ponemos en unas funcioncitas todo

get_adyacentes <- function(mapa = mapa, id_col = "cp_num") {
  
  # quiza sacar el st_intersects de la función sea mejor
  nb <-  st_intersects(mapa)
  
  get_nb <-  function(x){
      res <- mapa[[id_col]][x]
      res
  }
  
  adjacency_list <-  lapply(nb, get_nb)
  adjacency_list_names <-  mapa[[id_col]][1:length(adjacency_list)]
  names(adjacency_list) <- adjacency_list_names
  adjacency_list
}




get_mis_vecinos <-  function(mi_cp, cps_adyacentes){
    cp_simbol <-  as.symbol(mi_cp) # a simbolo para poder llamara cps_adyacentes[[`18814`]]
    mis_vecinos <-  cps_adyacentes[[cp_simbol]]
}



# mapa seleccionando solo tu provincia 


plot_cp_vecinos <-  function(mi_cp, cps_adyacentes, mapa){
    cp_simbol <-  as.symbol(mi_cp)
    mis_vecinos <-  cps_adyacentes[[cp_simbol]]
    mi_prop <-  stringr::str_sub(mi_cp, 1, 2) |> as.numeric()
    adyacentes <-  setdiff(mis_vecinos, mi_cp)
    mapa |> 
        filter(cp_2_num == mi_prop) |> 
        mutate(
            tipo = case_when(
                cp_num == mi_cp ~ "mi_cp", 
                cp_num %in% adyacentes ~ "mis_copostales_vecinos", 
                TRUE ~ "resto de codpostales"
            )
        ) |> 
        ggplot() +
        geom_sf(aes(fill = tipo)) +
        scale_fill_viridis_d()
}

Y ya podemos usarlo

cps_adyacentes <-  get_adyacentes(mapa, id_col = "cp_num")

(mis_cps_vecinos_pueblo <- get_mis_vecinos(18814, cps_adyacentes))
#> [1] 18815 18800 18817 18811 18814 18818 18816

Pintamos

plot_cp_vecinos(18814, cps_adyacentes, mapa) +
    labs(title = "CP núcleo principal Cortes de Baza")

# plaza cascorro

plot_cp_vecinos(28005, cps_adyacentes, mapa) +
    labs(title = "cp Plaza Cascorro")

# Pongo Carmona por ser el de mayor área 

(get_mis_vecinos(41410, cps_adyacentes))
#>  [1] 41620 41300 41449 41339 41420 41359 41610 41600 41320 41309 41520 41429
#> [13] 41410 41510 41310 41016 41500 41440

plot_cp_vecinos(41410, cps_adyacentes, mapa) +
    labs(title= "Carmona")

Notas

• Los códigos postales que uso están desactualizados. La capa de shapefile es un producto que vende Correos a precio no barato. Antes se podían descargar de CartoCiudad, pero ya no. Algún enlace interesante:
  • post de Raúl Vaquerizo
  • post de nosolosig
• Libros online sobre GIS en R:
  • r-spatial
  • geocomputation with R
• Applied Spatial Data Analysis with R . Un poco antiguo ya, no sé si Virgilio, Roger y Edzer están trabajando sobre nueva edición o no