Análisis de puntos de interés (POI) con OpenStreetMap y python
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?
# 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
nodeyare.Para aquellos elementos que cumplan la condición anterior, se almacena toda la información de sus tags.
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
# 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.
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)
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.
df_poi[df_poi.id == 21812709]
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.
estaciones_fuel = df_poi[df_poi['tag_key'].isin(['brand', 'operator'])]
estaciones_fuel = estaciones_fuel.groupby('id').first()
estaciones_fuel
# 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()
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.
replacements = {
'repsol moto stop': 'repsol',
'moto stop repsol': 'repsol',
'e. leclerc': 'e.leclerc'
}
estaciones_fuel['tag_value'] = estaciones_fuel['tag_value'].replace(replacements)
estaciones_fuel['tag_value'].value_counts()
# Creación del mapa con una ubicación inicial
# ==============================================================================
centro_madrid = [40.41678750299365, -3.7037881869889424]
#mapa = folium.Map(location=centro_madrid, zoom_start=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)
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.
# 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
# Añadir marcadores
# ==============================================================================
# for row in estaciones_fuel.itertuples():
# folium.Marker(
# location = [row.lat, row.lon],
# popup = row.tag_value,
# ).add_to(mapa)
# mapa
Identificar estaciones dentro de un radio¶
Para identificar elementos que están dentro de un determinado radio, hay que definir el centroide de referencia y calcular la distancia que hay desde él a cada uno de los elementos. Véase cómo se seleccionan, por ejemplo, todas las estaciones de servicio que están a menos de 8 km del centro de Madrid.
folium.Circle(
location = centro_madrid,
radius = 8_000, # metros
popup = '10 km radius',
fill = True,
color = '#B22222'
).add_to(mapa)
mapa
# Filtrado de elementos por distancia
# ==============================================================================
def calcular_distancia(ref_cord, target_cord):
'''
Función para calcular la distancia (km) que hay entre dos puntos de la superficie
terrestre.
'''
return distance(ref_cord, target_cord).km
# Cálculo de distancia al centro de Madrid
estaciones_fuel['distancias'] = estaciones_fuel.apply(
lambda x: calcular_distancia((x.lat, x.lon), centro_madrid),
axis = 1
)
# Filtrado de las que están a menos de 8 km del centro de Madrid
estaciones_fuel = estaciones_fuel[estaciones_fuel['distancias'] < 8]
estaciones_fuel
Mapa de calor¶
Los mapas de calor resultan muy útiles para identificar zonas con elevada densidad de elementos, en este caso, estaciones de servicio.
fig = Figure(width=600, height=400)
mapa = folium.Map(location=centro_madrid, zoom_start=11)
fig.add_child(mapa)
heatmap = HeatMap(
data = estaciones_fuel[['lat','lon']].to_numpy(),
radius = 20,
blur = 25,
)
heatmap.add_to(mapa)
mapa
Distancia mínima entre estaciones de servicio¶
A la hora de decidir la ubicación de un nuevo establecimiento, es importante conocer a qué distancia se encuentra el establecimiento más cercano que ofrece el mismo tipo de servicio. Con este objetivo, se procede a calcular, para cada estación de servicio, cuál es la distancia a la estación de servicio más próxima. Los pasos a seguir son:
Obtener todos los pares de estaciones de servicio.
Calcular la distancia entre todos los pares de estaciones de servicio.
Identificar la distancia mínima entre cada estación de servicio.
Para minimizar el número de cálculos, dado que la distancia de A a B es la misma que de B a A, es recomendable utilizar combinaciones en lugar de permutaciones. Esto significa que, una vez se ha calculado la distancia de la estación A a la estación B, no se repite el proceso de B a A.
# Cálculo de distancias entre todos los pares de estaciones
# ==============================================================================
combinaciones = itertools.combinations(estaciones_fuel.index, 2)
distancia = []
estacion_a = []
estacion_b = []
for combinacion in combinaciones:
coord_a = estaciones_fuel.loc[combinacion[0], ['lat', 'lon']]
coord_b = estaciones_fuel.loc[combinacion[1], ['lat', 'lon']]
estacion_a.append(combinacion[0])
estacion_b.append(combinacion[1])
distancia.append(distance(coord_a, coord_b).km)
Para facilitar el filtrado, se concatenan los resultados en ambas direcciones, es decir, si se ha calculado la distancia A-B, se añade también como B-A.
distancias = pd.DataFrame({
'estacion_a' : estacion_a + estacion_b,
'estacion_b' : estacion_b + estacion_a,
'distancia' : distancia + distancia
})
distancias.head(4)
distancia_minima = distancias.groupby('estacion_a')['distancia'].min()
distancia_minima.describe()
fig, ax = plt.subplots(figsize=(8, 4))
distancia_minima.hist(ax=ax)
ax.set_title('Distribución distancia mínima entre estaciones de servicio')
ax.set_xlabel('km');
Tras seleccionar únicamente las estaciones de servicio que se encuentran a un radio de menos de 8 km del centro de la ciudad, la distancia media es de aproximadamente 0.5 km. La estación más aislada está a 3.2 km de la estación más cercana y las dos estaciones más próximas están a 0.0039 km.
La distancia mínima de 3.9 metros es sospechosamente baja. Si se identifican cuáles son las estaciones, se observa que pertenecen a un nodo y a una área. Lo más probable es que se trate de un problema de duplicidad, en el que, el mismo elemento, se ha registrado como nodo y como área. En un estudio más detallado convendría eliminar duplicidades de este tipo, por ejemplo, detectando cuando la distancia no supera un mínimo razonable.
distancias.sort_values('distancia').head(1)
estaciones_fuel.loc[[2802306712, 551200196]]
Overpass-turbo¶
Extraer datos de OSM con overpass-turbo es muy útil cuando se quiere únicamente información sobre un subconjunto de elementos. Además de su asistente de consultas, se puede encontrar información detallada en su documentación.
De entre los formatos de descarga, se encuentra geojson, que puede leerse fácilmente con la función read_file de geopandas.
data = gpd.read_file('export_madrid_osm.geojson')
data.head(2)
Una vez cargados los datos con geopandas, con la función explore se crea un mapa sobre el que se representa la información disponible.
# data.explore(tooltip=['brand', 'operator'])
from branca.element import Figure
fig = Figure(width=600, height=400)
mapa = folium.Map(location=centro_madrid, zoom_start=11)
fig.add_child(mapa)
data.explore(
tooltip = ['brand', 'operator'],
m = mapa,
marker_kwds = {'radius':4, 'color': 'blue', 'fill': True}
)
Información de sesión¶
import session_info
session_info.show(html=False)
¿Cómo citar este documento?
Análisis de puntos de interes (POI) con OpenStreetMap y python por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/py40-puntos-interes-openstreetmap-python.html 
¿Te ha gustado el artículo? Tu ayuda es importante
Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊
Este contenido, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.
Se permite:
-
Compartir: copiar y redistribuir el material en cualquier medio o formato.
-
Adaptar: remezclar, transformar y crear a partir del material.
Bajo los siguientes términos:
-
Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.
-
NoComercial: No puedes utilizar el material para fines comerciales.
-
CompartirIgual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.