Extraer y analizar puntos de interés (POI) con OpenStreetMap OSM y python

Análisis de puntos de interés (POI) con OpenStreetMap y python

Joaquín Amat Rodrigo
Enero, 2022

Más sobre ciencia de datos: cienciadedatos.net

Introducción


OpenStreetMap (OSM) es un proyecto colaborativo para crear mapas libres. Su base de datos contiene información detallada sobre carreteras, caminos y muchos otros elementos de la cartografía mundial. Gracias a su principio de código abierto y libre, se ha convertido en una de las principales fuentes de datos para realizar proyectos en el campo de los Sistema de Información Geográfica (SIG o GIS).

Existen varias formas de acceder a estos datos, dos de las más fáciles son:

  • Descargando los ficheros .pbf con toda la información de una determinada región, por ejemplo desde la web Geofabrik.

  • Descargar únicamente los elementos de interés de una determinada región. Esto puede hacerse mediante el uso de overpass-turbo, una herramienta que, además de tener una interfaz web muy sencilla, dispone de un asistente para la creación de consultas. Las descargas desde overpass-turbo pueden ser en formato GeoJSON, GPX, KML y JSON.

En este documento, se muestra cómo descargar los datos disponibles en OSM y cómo extraer información sobre puntos de interés.

Caso de uso


Supóngase una empresa comercializadora de combustible que quiere abrir una nueva estación de servicio en la ciudad de Madrid. Para identificar potenciales ubicaciones, se quieren conocer características sobre la distribución geográfica de las estaciones de servicio que hay en la ciudad. Algunas de las preguntas que se quieren contestar son:

  • ¿Cuántas estaciones de servicio hay en la ciudad y dónde se encuentran?

  • ¿A qué compañías pertenecen las estaciones de servicio?

  • ¿Cuál es la distribución de distancias entre las estaciones de servicio en la ciudad?

  • ¿Qué estaciones están dentro de un determinado radio?



Librerías


Las librerías utilizadas en este documento son:

In [1]:
# Lectura datos osm
import osmium as osm

# Tratamiento de datos
import pandas as pd
import numpy as np
import itertools

# Tratamiento de datos geográficos y mapas
import matplotlib.pyplot as plt
import geopandas as gpd
from geopy.distance import distance
from shapely.geometry import Polygon
import folium
from folium.plugins import HeatMap
from branca.element import Figure

Extracción de puntos de interés de un fichero .pbf


Descarga de datos


Desde geofabrik pueden descargarse los archivos .pbf con la información relativa a múltiples regiones del mundo, entre ellas, las comunidades de España.

Elementos y tags en OSM


Para poder extraer la información almacenada en un archivo .pbf es importante entender cómo se registran los datos en OpenStreetMap. En OpenStreetMap existen 3 tipos de elementos principales:

  • nodes: representan puntos en el espacio definidos por una latitud y longitud.

  • ways: son listas ordenadas de nodos. Se utilizan para definir elementos lineales (carreteras, ríos...) o áreas (edificios, parques...). En este último caso, el primer y último nodo son el mismo.

  • relations: definen relaciones entre dos o más elementos.

Todo elemento (nodes, ways y relations) pueden tener asociado uno o más tags que describen sus características. Cada tag está formada por dos campos de texto llamados key y value. Por ejemplo, un nodo que tiene asociado el tag amenity=restaurant es un restaurante. OpenStreetMap permite registrar cualquier key y value, sin embargo, existen una serie de recomendaciones para tratar de estandarizar la información. Puede encontrarse un listado completo en Map_features.

Extracción de puntos de interés


El primer paso para extraer un determinado elemento cartográfico de un archivo .pbf es identificar cómo está almacenado dentro de OSM. En este ejemplo, se quiere extraer información sobre las estaciones de servicio. Consultando el listado Map_features se especifica que, las estaciones de servicio, están identificadas con el tag Tag:amenity=fuel y que están almacenadas como elementos de tipo node y way. Por lo tanto, de toda la información almacenada dentro del archivo .pbf, hay que extraer únicamente los nodes y ways que tengan el tag amenity con el valor fuel.

La librería Pyosmium contiene la clase SimpleHandler que permite leer y extraer información de ficheros OSM. Esta es una clase genérica, que se ha de completar con un método específico para cada elemento que se quiere extraer. En este caso, se necesita:

  • Para elementos de tipo node y area (este es un subtipo de way) comprobar si tienen un tag con key=amenity y value=fuel. De esto se encargan los métodos node y are.

  • Para aquellos elementos que cumplan la condición anterior, se almacena toda la información de sus tags.

In [2]:
from shapely.geometry import Polygon

class POIHandler(osm.SimpleHandler):
    '''
    Clase para extraer información de un archivo osm.pbf. Únicamente se extraen
    elementos identificados como 'node' o 'area'. Además, se puede aplicar un filtrado
    para seleccionar únicamente aquellos que tengan tags con un determinado key y
    value.
    
    La posición de las áreas se obtiene calculando el centroide del polígono que
    forman sus nodos.
    
    TODO: 'all' value to include all available values of a tag.
    
    Arguments
    ---------
    
    custom_filter: dict
        Diccionario con los tags y valores que han de tener los elementos para
        ser extraídos. Por ejemplo:
        
        `{'amenity': ['restaurant', 'bar']}` selecciona únicamente aquellos
        elementos que tengan el tag 'amenity' con valor 'restaurant' o 'bar'.
        
        `{'amenity': ['restaurant', 'bar'], 'building': ['car']}` selecciona
        únicamente aquellos elementos que tengan el tag 'amenity' con valor
        'restaurant' o 'bar', o los que tengan el tag 'building' con valor 'hotel'.
    '''
    
    def __init__(self, custom_filter=None):
        osm.SimpleHandler.__init__(self)
        self.osm_data = []
        self.custom_filter = custom_filter
        
        if self.custom_filter:
            for key, value in self.custom_filter.items():
                if isinstance(value, str):
                    self.custom_filter[key] = [value]
             
    def node(self, node):
        if self.custom_filter is None:
            name = node.tags.get('name', '')
            self.tag_inventory(node, 'node', name)
        else:
            if any([node.tags.get(key) in self.custom_filter[key] for key in self.custom_filter.keys()]):
                name = node.tags.get('name', '')
                self.tag_inventory(node, 'node', name)
                
    def area(self, area):
        if self.custom_filter is None:
            name = area.tags.get('name', '')
            self.tag_inventory(area, 'area', name)
        else:
            if any([area.tags.get(key) in self.custom_filter[key] for key in self.custom_filter.keys()]):
                name = area.tags.get('name', '')
                self.tag_inventory(area, 'area', name)

    def tag_inventory(self, elem, elem_type, name):
        if elem_type == 'node':
            for tag in elem.tags:
                self.osm_data.append([elem_type, 
                                       elem.id,
                                       name,
                                       elem.location.lon,
                                       elem.location.lat,
                                       pd.Timestamp(elem.timestamp),
                                       len(elem.tags),
                                       tag.k, 
                                       tag.v])
        if elem_type == 'area':
            try:
                # Se crea un Polygon con los nodos que forman el area para calcular
                # su centroide.
                nodes = list(elem.outer_rings())[0]
                polygon = Polygon([(node.lon, node.lat) for node in nodes])
                for tag in elem.tags:
                    self.osm_data.append([elem_type, 
                                           elem.id,
                                           name,
                                           polygon.centroid.x,
                                           polygon.centroid.y,
                                           pd.Timestamp(elem.timestamp),
                                           len(elem.tags),
                                           tag.k, 
                                           tag.v])
            except:
                pass
In [3]:
# Extracción de elementos identificados con 'amenity'='fuel'
# ==============================================================================
poi_handler = POIHandler(custom_filter={'amenity':['fuel']})
poi_handler.apply_file('madrid-latest.osm.pbf')

Como resultado de la extracción se obtiene una lista en la que, cada elemento, es a su vez una lista con la información asociada a un nodo o área. Para facilitar el manejo de los datos extraídos, se almacenan en formato de pandas dataframe.

In [4]:
colnames = ['type', 'id', 'name', 'lon', 'lat', 'timestamp','n_tags', 'tag_key',
            'tag_value']
df_poi = pd.DataFrame(poi_handler.osm_data, columns=colnames)
df_poi.head(4)
Out[4]:
type id name lon lat timestamp n_tags tag_key tag_value
0 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 name Repsol AutoGas
1 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 note Only LPG / Sólo butano
2 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 brand Repsol
3 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 source yahoo_maps;survey

Cada fila contiene un atributo (tag_key y tag_value) asociado a un nodo. OpenStreetMap tiene un sistema de tags libre por lo que, un mismo nodo, puede tener un número ilimitado de atributos asociados. Véase por ejemplo todos los 14 tags asociados al nodo con id 21812709.

In [5]:
df_poi[df_poi.id == 21812709]
Out[5]:
type id name lon lat timestamp n_tags tag_key tag_value
0 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 name Repsol AutoGas
1 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 note Only LPG / Sólo butano
2 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 brand Repsol
3 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 source yahoo_maps;survey
4 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 amenity fuel
5 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 fuel:e10 no
6 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 fuel:lpg yes
7 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 operator Repsol
8 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 fuel:diesel no
9 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 brand:wikidata Q174747
10 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 fuel:octane_91 no
11 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 fuel:octane_95 no
12 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 fuel:octane_98 no
13 node 21812709 Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 brand:wikipedia es:Repsol

Número de estaciones de servicio por compañía


Acorde a la documentación de OSM, la marca de la compañía a la que pertenece la estación de servicio está identificada en los tags brand o operator. Algunos nodos tienen informados los dos tags, y otros solo uno. Para evitar duplicidades en el caso de que se disponga de ambos, se selecciona únicamente el primero.

In [6]:
estaciones_fuel = df_poi[df_poi['tag_key'].isin(['brand', 'operator'])]
estaciones_fuel = estaciones_fuel.groupby('id').first()
estaciones_fuel
Out[6]:
type name lon lat timestamp n_tags tag_key tag_value
id
3039509 area -3.725623 40.372760 2020-11-17 13:26:15+00:00 2 operator CEPSA
21812709 node Repsol AutoGas -3.732854 40.399588 2020-04-13 00:08:46+00:00 14 brand Repsol
21812714 node Cepsa -3.761717 40.385122 2021-05-05 13:00:45+00:00 11 brand Cepsa
22554715 node -3.743247 40.373537 2017-04-28 15:54:13+00:00 3 operator ESSO
23484045 area Cepsa -3.627248 40.450873 2020-10-11 00:55:36+00:00 9 brand Cepsa
... ... ... ... ... ... ... ... ...
7617537917 node Plenoil -4.004376 40.637964 2021-01-03 10:27:19+00:00 6 brand Plenoil
8340046013 node Shell -3.482732 40.310970 2021-01-19 21:26:41+00:00 5 brand Shell
8669317316 node AS 24 -3.665193 40.164882 2021-04-26 17:32:00+00:00 6 brand AS 24
8858593415 node -3.760212 40.654533 2021-06-22 23:44:33+00:00 2 brand Galp
9010043453 node BP -3.633494 40.370196 2021-08-16 11:22:22+00:00 5 brand BP

600 rows × 8 columns

In [7]:
# Para evitar problemas de mayúsculas y minúsculas, se normalizan los nombres
# convirtiéndolos en minúsculas.
estaciones_fuel['tag_value'] = estaciones_fuel['tag_value'].str.lower()
estaciones_fuel['tag_value'].value_counts()
Out[7]:
repsol                246
cepsa                 106
bp                     65
galp                   53
shell                  51
carrefour              14
ballenoil              12
gas natural fenosa      8
campsa                  6
alcampo                 6
q8                      3
plenoil                 3
petroprix               3
esso                    3
top-oil                 2
petronor                2
meroil                  2
e.leclerc               2
as 24                   1
moto stop repsol        1
e. leclerc              1
ghc                     1
sp                      1
fast fuel               1
simon grup              1
repsol moto stop        1
supercor                1
simply                  1
la ballena              1
agip                    1
aliara energía          1
Name: tag_value, dtype: int64

Algunas comercializadoras están registradas con distintos nombres, por ejemplo: repsol, repsol moto stop y moto stop repsol. Se renombran los valores para unificar las marcas.

In [8]:
replacements = {
    'repsol moto stop': 'repsol',
    'moto stop repsol': 'repsol',
    'e. leclerc': 'e.leclerc'
}
estaciones_fuel['tag_value'] = estaciones_fuel['tag_value'].replace(replacements)
In [9]:
estaciones_fuel['tag_value'].value_counts()
Out[9]:
repsol                248
cepsa                 106
bp                     65
galp                   53
shell                  51
carrefour              14
ballenoil              12
gas natural fenosa      8
alcampo                 6
campsa                  6
q8                      3
plenoil                 3
petroprix               3
e.leclerc               3
esso                    3
petronor                2
meroil                  2
top-oil                 2
la ballena              1
simply                  1
simon grup              1
aliara energía          1
fast fuel               1
agip                    1
ghc                     1
as 24                   1
supercor                1
sp                      1
Name: tag_value, dtype: int64

Representación en mapas


La librería folium permite crear mapas interactivos basados en leaflet. Además de la cartografía del mapa, permite superponer elementos como puntos y polígonos.

In [10]:
# Creación del mapa con una ubicación inicial
# ==============================================================================
centro_madrid = [40.41678750299365, -3.7037881869889424]

#mapa = folium.Map(location=centro_madrid, zoom_start=11)
In [11]:
# Para evitar los espacios en blanco alrededor del mapa se puede utilizar la
# librería branca
fig = Figure(width=600, height=400)
mapa = folium.Map(location=centro_madrid, zoom_start=11)
fig.add_child(mapa)
Out[11]:

Una vez creado el mapa, se añaden los puntos de interés. Folium dispone de varias opciones a la hora de identificar elementos en un mapa, dos de ellas son Marker y CircleMarker.

In [12]:
# Añadir círculos de tamaño fijo
# ==============================================================================
for row in estaciones_fuel.itertuples():
    folium.CircleMarker(
      location = [row.lat, row.lon],
      radius = 4,
      color = 'blue',
      fill = True,
      popup = row.tag_value,
   ).add_to(mapa)
    
mapa
Out[12]: