Estimando nodos OpenStreetMap por país

OpenStreetMap es un proyecto remarcable. Antes de su existencia, la producción de geoinformación era una actividad exclusiva de instituciones gubernamentales o agencias comerciales, de manera que los usuarios enfrentaban limitaciones en forma de licencias o costos. Desde 2004, en cambio, existe una plataforma que -con ayuda de las tecnologías que volvieron posible el compartir datos libremente a través de internet- recibe y entrega geoinformación de manera gratuita y ubicua.

No obstante, la utilización de OpenStreetMap implica ciertos inconvenientes derivados de la propia naturaleza de su funcionamiento. Un inconveniente es conocer su estado actual: como la base de datos cambia constantemente -con adiciones y modificaciones de usuarios en todo el mundo- es complicado conocer, por ejemplo, cuántas calles han sido ingresadas en el mapa hasta el momento. Con todo, la página OpenStreetMap Wiki Stats presenta diversos recursos (algunos de ellos incluso actualizados diariamente) que permiten tener una idea del estado actual del proyecto.

Mientras realizaba mi tesis de grado, en determinado momento necesité un recurso muy específico: un mapa coroplético de elementos OpenStreetMap por país; es decir, un mapa donde cada país está coloreado de acuerdo a cuántos elementos han sido contribuidos en su territorio. Para tener una idea de cómo luce tal mapa, comparto uno hallado en una de mis referencias (Jokar Arsanjani et al. 2015):

Nodos OpenStreetMap creados por país en octubre de 2014; © 2015 Springer International Publishing

Si bien para elaborar este mapa se contabilizó únicamente los nodos creados en un mes, de todas maneras se demuestra correctamente el panorama global de las contribuciones a OpenStreetMap: en los países europeos se encuentran las comunidades más activas; Canadá, Estados Unidos, Brasil, Rusia y Australia reciben gran cantidad de contribuciones, acorde a sus extensos territorios; y los niveles de contribución varían considerablemente en Latinoamérica, África y Asia.

Es importante tener presente que en OpenStreetMap existen tres elementos básicos: nodos (nodes), vías (ways) y relaciones (relations). Pero tanto las vías como las relaciones se componen, en última instancia, de nodos; por esto, un mapa que contabilice los nodos debería representar correctamente la distribución de las contribuciones. Una forma de generar tal mapa es descargar los datos de cada país -la página Geofabrik, por ejemplo, contiene extractos por país- y contar los nodos en cada uno, pero este proceso resultaría extenuante. Esta es la motivación detrás de la siguiente metodología para estimar el número de nodos por país aplicando web scraping, una técnica de minería de datos que consiste en utilizar software para extraer el contenido de una página web.

La fuente de datos es el sitio OSMstats desarrollado por Pascal Neis, cuya investigación en el ámbito de OpenStreetMap (Neis and Zipf 2012) también fue clave en mi tesis. Este sitio ofrece estadísticas sobre diferentes aspectos del proyecto, con una periodicidad diaria desde el 1 de noviembre de 2011; particularmente, la pestaña Countries presenta estimaciones del número de contribuyentes y de nodos creados / modificados / eliminados en cada país, cada día. El primer paso en esta metodología es escribir una función de extracción; en este caso, se utilizó los paquetes httr, xml2 y purrr:

library(purrr)

get_day = possibly(function(day, url){
   
   message("Day: ", day)
   
   df = paste0(url, day) %>%
      httr::GET() %>%
      httr::content(encoding = "UTF-8") %>%
      xml2::as_list() %>%
      plucking() %>%
      map(framing) %>%
      reduce(rbind)
   
   df$date = day
   
   return(df)}, NULL)

Las funciones plucking() y framing() no se encuentran en ningún paquete; en realidad, fueron declaradas previamente con el propósito de recortar y tabular los datos, respectivamente. Los pasos más importantes son aquellos entre httr::GET() y xml2::as_list() pues es allí donde se extrae y transforma el contenido de la página expuesta al web scraping. Adicionalmente, toda la función de extracción fue insertada en purrr::possibly(); se trata de una precaución que devuelve NULL en caso de que la ejecución falle.

Ahora, esta función extrae los datos de una fecha específica, pero su verdadera utilidad aparecerá al iterarla en varias fechas. Considerando que los datos empiezan el 1 de noviembre de 2011, se decidió iterar por 2923 días, hasta el 1 de noviembre de 2019; es decir, se contabilizará el crecimiento de OpenStreetMap en ocho años.

neis = "https://osmstats.neis-one.org/?item=countries&date="

days = as.Date("2011-11-1") + 0:2922

countries = map(days, get_day, url = neis) %>% reduce(rbind)

Los datos contienen 738001 observaciones, de 2917 días en 253 países o territorios. Seis de los días del año 2018 carecen de datos: 27 de julio, 3 y 24 de agosto, 10 y 29 de septiembre, y 22 de octubre; desconocemos la causa de este hecho. El siguiente paso será sumar los nodos creados -restando al mismo tiempo los eliminados- en cada país; esto es posible implementando dplyr:

library(dplyr)

countries_nodes = group_by(countries, country) %>%
   mutate(across(created:deleted, as.double)) %>%
   summarise(nodes = sum(created - deleted), dates = n_distinct(date))

summary(countries_nodes)

##    country              nodes               dates
##  Length:253         Min.   :       -2   Min.   :2917
##  Class :character   1st Qu.:   259873   1st Qu.:2917
##  Mode  :character   Median :  2756173   Median :2917
##                     Mean   : 16590188   Mean   :2917
##                     3rd Qu.: 13561870   3rd Qu.:2917
##                     Max.   :388753376   Max.   :2917

Finalmente, para generar el mapa, los datos calculados en el paso anterior deben ser relacionados con la representación geográfica de los países. El dataset World contiene polígonos simplificados de territorios en todo el mundo; al ser un dataset de clase sf, el paquete homónimo deberá ser incluido para garantizar un procesamiento correcto. Luego es posible trazar el mapa vía ggplot2:

library(sf)
library(ggplot2)

data("World", package = "tmap")

left_join(World, countries_nodes, by = c("name" = "country")) %>%
   ggplot(aes(fill = nodes)) + geom_sf(color = NA) + theme_minimal()

Figure 1: Primer intento de generar nuestro mapa

En este primer intento se detectó que algunos países aparecen sin datos, como las dos repúblicas de Korea y las dos del Congo. El problema se desprende de los diferentes nombres con que estos países figuran en World y en OSMstats; por ende, para solucionar este problema se debe averiguar cuáles son esos nombres, ejecutando setdiff(World$name, countries_nodes$country):

##  [1] "Fr. S. Antarctic Lands" "Bahamas"                "Bosnia and Herz."
##  [4] "Central African Rep."   "Cote d'Ivoire"          "Dem. Rep. Congo"
##  [7] "Congo"                  "N. Cyprus"              "Czech Rep."
## [10] "Dominican Rep."         "Falkland Is."           "Gambia"
## [13] "Eq. Guinea"             "Korea"                  "Kosovo"
## [16] "Lao PDR"                "Myanmar"                "Dem. Rep. Korea"
## [19] "Palestine"              "W. Sahara"              "S. Sudan"
## [22] "Solomon Is."            "Somaliland"             "Timor-Leste"
## [25] "Tanzania"

Entonces los nombres obtenidos de OSMstats deben ser modificados para coincidir con los de arriba.

countries_nodes = mutate(countries_nodes, across(country, ~case_when(
   . == "Bosnia and Herzegovina"              ~ "Bosnia and Herz.",
   . == "Central African Republic"            ~ "Central African Rep.",
   . == "Congo-Brazzaville"                   ~ "Congo",
   . == "Congo-Kinshasa"                      ~ "Dem. Rep. Congo",
   . == "Czech Republic"                      ~ "Czech Rep.",
   . == "Dominican Republic"                  ~ "Dominican Rep.",
   . == "Equatorial Guinea"                   ~ "Eq. Guinea",
   . == "Falkland Islands (Islas Malvinas)"   ~ "Falkland Is.",
   . == "French Southern and Antarctic Lands" ~ "Fr. S. Antarctic Lands",
   . == "Ivory Coast"                         ~ "Cote d'Ivoire",
   . == "Laos"                                ~ "Lao PDR",
   . == "Myanmar (Burma)"                     ~ "Myanmar",
   . == "North Korea"                         ~ "Dem. Rep. Korea",
   . == "Republic of Kosovo"                  ~ "Kosovo",
   . == "Solomon Islands"                     ~ "Solomon Is.",
   . == "South Korea"                         ~ "Korea",
   . == "South Sudan"                         ~ "S. Sudan",
   . == "The Bahamas"                         ~ "Bahamas",
   . == "The Gambia"                          ~ "Gambia",
   . == "United Republic of Tanzania"         ~ "Tanzania",
   . == "Western Sahara"                      ~ "W. Sahara",
   TRUE ~ .)))

setdiff(World$name, countries_nodes$country)

## [1] "N. Cyprus"   "Palestine"   "Somaliland"  "Timor-Leste"

Chipre del Norte, Palestina, Somalilandia y Timor Oriental son países que, tal vez debido al limitado reconocimiento, no figuran en OSMstats y evidentemente persistirán vacíos en el mapa. En cambio, un total de 80 territorios que sí existen en OSMstats, no aparecen en World por tener una pequeña extensión: Andorra, Singapur y las Antillas Menores son algunos ejemplos de esta situación. Como consecuencia, del total de 4197.32 millones de nodos extraídos, 4180.61 millones sí se encuentran representados en el mapa, lo cual es el 99.6 %. Habiendo aclarado esto, podemos unir los datos nuevamente y trazar el mapa con algunas mejoras.

Figure 2: Segundo intento de generar nuestro mapa 🎉

Para concluir, es necesario recordar que OpenStreetMap inició en 2004; por ende, esta metodología es incapaz de contabilizar todos los nodos existentes. ¿Cuántos nodos faltan? Bueno, el propio sitio OSMstats reportó, el 1 de noviembre de 2019, la existencia de 5559.59 millones de nodos en la base de datos; quiere decir que esta metodología contabilizó el 75.49 % de todos los nodos. El mismo sitio reportó 1248.29 millones de nodos el 1 de noviembre de 2011 (lo cual es el 22.45 % del valor de 2019). Sumando ambos porcentajes alcanzamos el 97.94 %; este es el porcentaje de nodos OpenStreetMap cuya existencia podemos justificar.

A través del web scraping en el sitio de Pascal Neis, no solo tracé el mapa coroplético que necesitaba; también generé un dataset considerablemente grande e interante. Por supuesto, en futuros artículos lo seguiré utilizando para aprender más de este maravilloso proyecto. Por lo pronto, averigüemos cuáles son los diez países con más nodos:

slice_max(countries_nodes, nodes, n = 10)

## # A tibble: 10 x 3
##    country           nodes dates
##    <chr>             <dbl> <int>
##  1 United States 388753376  2917
##  2 Russia        294001095  2917
##  3 Canada        231464569  2917
##  4 France        230694245  2917
##  5 Germany       220406582  2917
##  6 Indonesia     154173689  2917
##  7 Italy         143362158  2917
##  8 Poland        120794989  2917
##  9 Brazil        112646013  2917
## 10 Japan          98296294  2917