En este artículo, combino trabajos anteriores sobre accesibilidad urbana o caminabilidad con datos de fuente abierta sobre la ubicación de dispositivos desfibriladores públicos. Además, incorporo datos de población global y el sistema de red H3 de Uber para estimar la proporción de la población al alcance razonable de cualquier dispositivo dentro de Budapest y Viena.
La raíz de la accesibilidad urbana, o transitabilidad, radica en un cálculo basado en gráficos que mide la distancia euclidiana (transformándola en minutos de caminata, suponiendo una velocidad constante y sin atascos ni obstáculos). Los resultados de dichos análisis pueden decirnos qué tan fácil es llegar a tipos específicos de servicios desde cada lugar de la ciudad. Para ser más precisos, de cada nodo dentro de la red de carreteras de la ciudad, pero debido a la gran cantidad de cruces de carreteras, esta aproximación es prácticamente insignificante.
En este estudio de caso actual, me centro en un tipo particular de punto de interés (POI): la ubicación de los dispositivos desfibriladores. Si bien el Portal de Datos Abiertos del gobierno austriaco comparte registros oficiales al respecto, en Hungría solo pude obtener un conjunto de datos de fuentes colectivas de cobertura inferior a la mitad, que, con suerte, crecerá más adelante tanto en tamaño absoluto como en cobertura de datos.
En la primera sección de mi artículo crearé el mapa de accesibilidad de cada ciudad, visualizando el tiempo necesario para llegar a las unidades desfibriladoras más cercanas en un radio de 2,5 km a una velocidad de carrera de 15 km/h. Luego, dividiré las ciudades en cuadrículas hexagonales utilizando la biblioteca H3 de Uber para calcular el tiempo promedio de accesibilidad al desfibrilador para cada celda de la cuadrícula. También estimo el nivel de población en cada celda hexagonal siguiendo mi anterior artículo. Finalmente, los combino y calculo la fracción de la población accesible en función del tiempo de accesibilidad (en ejecución).
Como descargo de responsabilidad, quiero enfatizar que de ninguna manera soy un experto médico capacitado y no tengo la intención de adoptar una postura sobre la importancia de los dispositivos desfibriladores en comparación con otros medios de desfibrilador. vida apoyo. Sin embargo, partiendo del sentido común y de los principios de planificación urbana, supongo que cuanto más fácil sea llegar a dichos dispositivos, mejor.
Como siempre, me gusta empezar explorando los tipos de datos que uso. Primero, recopilaré los límites administrativos de las ciudades en las que estudio: Budapest, Hungría y Viena, Austria.
Luego, basándose en un artículo anterior mío sobre cómo procesar datos de población rasterizados, agrego información de población a nivel de ciudad desde el centro WorldPop. Finalmente, incorporo datos gubernamentales oficiales sobre dispositivos desfibriladores en Viena y mi propia versión extraída de la web de los mismos, aunque fuentes abarrotadas e intrínsecamente incompletas, para Budapest.
1.1. Límites administrativos
Primero, consulto los límites administrativos de Budapest y Viena desde Abrir mapa de calles utilizando el OSMNx biblioteca:
import osmnx as ox # version: 1.0.1
import matplotlib.pyplot as plt # version: 3.7.1admin = {}
cities = ['Budapest', 'Vienna']
f, ax = plt.subplots(1,2, figsize = (15,5))
# visualize the admin boundaries
for idx, city in enumerate(cities):
admin[city] = ox.geocode_to_gdf(city)
admin[city].plot(ax=ax[idx],color='none',edgecolor= 'k', linewidth = 2) ax[idx].set_title(city, fontsize = 16)
El resultado de este bloque de código:
1.2. Datos de población
En segundo lugar, siguiendo los pasos de este artículo, Creé la cuadrícula de población en formato de datos vectoriales para ambas ciudades, basándose en la base de datos demográfica en línea WorldPop. Sin repetir los pasos, simplemente leo los archivos de salida de ese proceso que contienen información de población de estas ciudades.
Además, para que las cosas se vean bien, creé un mapa de colores del color de 2022, Very Peri, usando Matplotlib y un script rápido de ChatGPT.
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormapvery_peri = '#8C6BF3'
second_color = '#6BAB55'
colors = [second_color, very_peri ]
n_bins = 100
cmap_name = "VeryPeri"
colormap = LinearSegmentedColormap.from_list(cmap_name, colors, N=n_bins)
import geopandas as gpd # version: 0.9.0demographics = {}
f, ax = plt.subplots(1,2, figsize = (15,5))
for idx, city in enumerate(cities):
demographics[city] = gpd.read_file(city.lower() + \
'_population_grid.geojson')[['population', 'geometry']]
admin[city].plot(ax=ax[idx], color = 'none', edgecolor = 'k', \
linewidth = 3)
demographics[city].plot(column = 'population', cmap = colormap, \
ax=ax[idx], alpha = 0.9, markersize = 0.25)
ax[idx].set_title(city)
ax[idx].set_title('Population density\n in ' + city, fontsize = 16)
ax[idx].axis('off')
El resultado de este bloque de código:
1.3. Ubicaciones de desfibriladores
En tercer lugar, recopilé datos de ubicación sobre los desfibriladores disponibles en ambas ciudades.
Para Viena, descargué este conjunto de datos del Portal oficial de datos abiertos del gobierno austriaco. que contiene la ubicación del punto de 1044 unidades:
Si bien en Budapest/Hungría no existe un portal oficial de datos abiertos, la Fundación Nacional Húngara del Corazón gestiona un sitio web colaborativo donde los operadores pueden actualizar la ubicación de sus unidades desfibriladoras. Su base de datos a nivel nacional consta de 677 unidades; sin embargo, su descargo de responsabilidad dice que conocen al menos mil unidades funcionando en el país y están esperando que sus propietarios las carguen. Con un simple rastreador web, descargué la ubicación de cada una de las 677 unidades registradas y filtré los datos establecidos hasta aquellas en Budapest, lo que resultó en un conjunto de 148 unidades.
# parse the data for each city
gdf_units= {}gdf_units['Vienna'] = gpd.read_file('DEFIBRILLATOROGD')
gdf_units['Budapest'] = gpd.read_file('budapest_defibrillator.geojson')
for city in cities:
gdf_units[city] = gpd.overlay(gdf_units[city], admin[city])
# visualize the units
f, ax = plt.subplots(1,2, figsize = (15,5))
for idx, city in enumerate(cities):
admin[city].plot(ax=ax[idx],color='none',edgecolor= 'k', linewidth = 3)
gdf_units[city].plot( ax=ax[idx], alpha = 0.9, color = very_peri, \
markersize = 6.0)
ax[idx].set_title('Locations of defibrillator\ndevices in ' + city, \
fontsize = 16)
ax[idx].axis('off')
El resultado de este bloque de código:
A continuación, concluí este gran artículo escrito por Nick Jones en 2018 sobre cómo calcular la accesibilidad de los peatones:
import os
import pandana # version: 0.6
import pandas as pd # version: 1.4.2
import numpy as np # version: 1.22.4
from shapely.geometry import Point # version: 1.7.1
from pandana.loaders import osmdef get_city_accessibility(admin, POIs):
# walkability parameters
walkingspeed_kmh = 15
walkingspeed_mm = walkingspeed_kmh * 1000 / 60
distance = 2500
# bounding box as a list of llcrnrlat, llcrnrlng, urcrnrlat, urcrnrlng
minx, miny, maxx, maxy = admin.bounds.T[0].to_list()
bbox = [miny, minx, maxy, maxx]
# setting the input params, going for the nearest POI
num_pois = 1
num_categories = 1
bbox_string = '_'.join([str(x) for x in bbox])
net_filename = 'data/network_{}.h5'.format(bbox_string)
if not os.path.exists('data'): os.makedirs('data')
# precomputing nework distances
if os.path.isfile(net_filename):
# if a street network file already exists, just load the dataset from that
network = pandana.network.Network.from_hdf5(net_filename)
method = 'loaded from HDF5'
else:
# otherwise, query the OSM API for the street network within the specified bounding box
network = osm.pdna_network_from_bbox(bbox[0], bbox[1], bbox[2], bbox[3])
method = 'downloaded from OSM'
# identify nodes that are connected to fewer than some threshold of other nodes within a given distance
lcn = network.low_connectivity_nodes(impedance=1000, count=10, imp_name='distance')
network.save_hdf5(net_filename, rm_nodes=lcn) #remove low-connectivity nodes and save to h5
# precomputes the range queries (the reachable nodes within this maximum distance)
# so, as long as you use a smaller distance, cached results will be used
network.precompute(distance + 1)
# compute accessibilities on POIs
pois = POIs.copy()
pois['lon'] = pois.geometry.apply(lambda g: g.x)
pois['lat'] = pois.geometry.apply(lambda g: g.y)
pois = pois.drop(columns = ['geometry'])
network.init_pois(num_categories=num_categories, max_dist=distance, max_pois=num_pois)
network.set_pois(category='all', x_col=pois['lon'], y_col=pois['lat'])
# searches for the n nearest amenities (of all types) to each node in the network
all_access = network.nearest_pois(distance=distance, category='all', num_pois=num_pois)
# transform the results into a geodataframe
nodes = network.nodes_df
nodes_acc = nodes.merge(all_access[[1]], left_index = True, right_index = True).rename(columns = {1 : 'distance'})
nodes_acc['time'] = nodes_acc.distance / walkingspeed_mm
xs = list(nodes_acc.x)
ys = list(nodes_acc.y)
nodes_acc['geometry'] = [Point(xs[i], ys[i]) for i in range(len(xs))]
nodes_acc = gpd.GeoDataFrame(nodes_acc)
nodes_acc = gpd.overlay(nodes_acc, admin)
nodes_acc[['time', 'geometry']].to_file(city + '_accessibility.geojson', driver = 'GeoJSON')
return nodes_acc[['time', 'geometry']]
accessibilities = {}
for city in cities:
accessibilities[city] = get_city_accessibility(admin[city], gdf_units[city])
for city in cities:
print('Number of road network nodes in ' + \
city + ': ' + str(len(accessibilities[city])))
Este bloque de código genera el número de nodos de la red de carreteras en Budapest (116.056) y en Viena (148.212).
Ahora visualiza los mapas de accesibilidad:
for city in cities:
f, ax = plt.subplots(1,1,figsize=(15,8))
admin[city].plot(ax=ax, color = 'k', edgecolor = 'k', linewidth = 3)
accessibilities[city].plot(column = 'time', cmap = 'RdYlGn_r', \
legend = True, ax = ax, markersize = 2, alpha = 0.5)
ax.set_title('Defibrillator accessibility in minutes\n' + city, \
pad = 40, fontsize = 24)
ax.axis('off')
Este bloque de código genera las siguientes figuras:
En este punto, tengo datos tanto de población como de accesibilidad; Sólo tengo que juntarlos. El único truco es que sus unidades espaciales difieren:
- La accesibilidad se mide y se adjunta a cada nodo dentro de la red vial de cada ciudad.
- Los datos de población se derivan de una cuadrícula ráster, ahora descrita por el punto de interés del centroide de cada cuadrícula ráster.
Si bien rehabilitar la cuadrícula ráster original puede ser una opción, con la esperanza de lograr una universalidad más pronunciada (y agregar un poco de mi gusto personal), ahora mapeo estos dos tipos de conjuntos de datos de puntos en el Sistema de red H3 de Uber Para aquellos que no lo han usado antes, por ahora basta con saber que se trata de un sistema de indexación espacial elegante y eficiente que utiliza mosaicos hexagonales. Y para leer más, ¡haga clic en este enlace!
3.1. Creando celdas H3
Primero, cree una función que divida una ciudad en hexágonos en cualquier resolución dada:
import geopandas as gpd
import h3 # version: 3.7.3
from shapely.geometry import Polygon # version: 1.7.1
import numpy as npdef split_admin_boundary_to_hexagons(admin_gdf, resolution):
coords = list(admin_gdf.geometry.to_list()[0].exterior.coords)
admin_geojson = {"type": "Polygon", "coordinates": [coords]}
hexagons = h3.polyfill(admin_geojson, resolution, \
geo_json_conformant=True)
hexagon_geometries = {hex_id : Polygon(h3.h3_to_geo_boundary(hex_id, \
geo_json=True)) for hex_id in hexagons}
return gpd.GeoDataFrame(hexagon_geometries.items(), columns = ['hex_id', 'geometry'])
resolution = 8
hexagons_gdf = split_admin_boundary_to_hexagons(admin[city], resolution)
hexagons_gdf.plot()
El resultado de este bloque de código:
Ahora, vea algunas resoluciones diferentes:
for resolution in [7,8,9]:admin_h3 = {}
for city in cities:
admin_h3[city] = split_admin_boundary_to_hexagons(admin[city], resolution)
f, ax = plt.subplots(1,2, figsize = (15,5))
for idx, city in enumerate(cities):
admin[city].plot(ax=ax[idx], color = 'none', edgecolor = 'k', \
linewidth = 3)
admin_h3[city].plot( ax=ax[idx], alpha = 0.8, edgecolor = 'k', \
color = 'none')
ax[idx].set_title(city + ' (resolution = '+str(resolution)+')', \
fontsize = 14)
ax[idx].axis('off')
El resultado de este bloque de código:
¡Mantengamos la resolución 9!
3.2. Asignar valores a celdas h3
Ahora tengo nuestras dos ciudades en un formato de cuadrícula hexagonal. A continuación, asignaré los datos de población y accesibilidad a las celdas hexagonales en función de en qué celdas de la cuadrícula se encuentre cada geometría de punto. Para esto, la función de unión de GeoPandasa, que realiza una buena unión espacial, es una buena opción.
Además, como tenemos más de 100.000 nodos de red de carreteras en cada ciudad y miles de centroides de cuadrícula de población, lo más probable es que haya múltiples puntos de interés asignados en cada celda de la cuadrícula hexagonal. Por lo tanto, será necesaria la agregación. Como la población es una cantidad aditiva, agregaré niveles de población dentro del mismo hexágono resumiéndolos. Sin embargo, la accesibilidad no es extensa, por lo que calcularía el tiempo promedio de accesibilidad al desfibrilador para cada mosaico.
demographics_h3 = {}
accessibility_h3 = {}for city in cities:
# do the spatial join, aggregate on the population level of each \
# hexagon, and then map these population values to the grid ids
demographics_dict = gpd.sjoin(admin_h3[city], demographics[city]).groupby(by = 'hex_id').sum('population').to_dict()['population']
demographics_h3[city] = admin_h3[city].copy()
demographics_h3[city]['population'] = demographics_h3[city].hex_id.map(demographics_dict)
# do the spatial join, aggregate on the population level by averaging
# accessiblity times within each hexagon, and then map these time score # to the grid ids
accessibility_dict = gpd.sjoin(admin_h3[city], accessibilities[city]).groupby(by = 'hex_id').mean('time').to_dict()['time']
accessibility_h3[city] = admin_h3[city].copy()
accessibility_h3[city]['time'] = \
accessibility_h3[city].hex_id.map(accessibility_dict)
# now show the results
f, ax = plt.subplots(2,1,figsize = (15,15))
demographics_h3[city].plot(column = 'population', legend = True, \
cmap = colormap, ax=ax[0], alpha = 0.9, markersize = 0.25)
accessibility_h3[city].plot(column = 'time', cmap = 'RdYlGn_r', \
legend = True, ax = ax[1])
ax[0].set_title('Population level\n in ' + city, fontsize = 16)
ax[1].set_title('Defibrillator reachability time\n in ' + city, \
fontsize = 16)
for ax_i in ax: ax_i.axis('off')
Los resultados de este bloque de código son las siguientes figuras:
En este último paso, estimaré la fracción de la población accesible desde la unidad desfibriladora más cercana dentro de un cierto período de tiempo. Aquí sigo aprovechando el ritmo de carrera relativamente rápido de 15 km/h y el límite de distancia de 2,5 km.
Desde la perspectiva técnica, fusiono los marcos de datos de tiempo de accesibilidad y población de nivel H3 y luego hago un umbral simple en la dimensión de tiempo y una suma en la dimensión de población.
f, ax = plt.subplots(1,2, figsize = (15,5))for idx, city in enumerate(cities):
total_pop = demographics_h3[city].population.sum()
merged = demographics_h3[city].merge(accessibility_h3[city].drop(columns =\
['geometry']), left_on = 'hex_id', right_on = 'hex_id')
time_thresholds = range(10)
population_reached = [100*merged[merged.time<limit].population.sum()/total_pop for limit in time_thresholds]
ax[idx].plot(time_thresholds, population_reached, linewidth = 3, \
color = very_peri)
ax[idx].set_xlabel('Reachability time (min)', fontsize = 14, \
labelpad = 12)
ax[idx].set_ylabel('Fraction of population reached (%)', fontsize = 14, labelpad = 12)
ax[idx].set_xlim([0,10])
ax[idx].set_ylim([0,100])
ax[idx].set_title('Fraction of population vs defibrillator\naccessibility in ' + city, pad = 20, fontsize = 16)
El resultado de este bloque de código son las siguientes figuras:
Al interpretar estos resultados, me gustaría enfatizar que, por un lado, la accesibilidad al desfibrilador puede no estar directamente relacionada con la tasa de supervivencia al ataque cardíaco; juzgar ese efecto está más allá de mi experiencia y del alcance de este proyecto. Además, los datos utilizados para Budapest son fuentes conscientemente incompletas y abarrotadas, a diferencia de la fuente de datos oficial de Austria.
Después de los descargos de responsabilidad, ¿qué vemos? Por un lado, vemos que en Budapest alrededor del 75-80% de la población puede acceder a un dispositivo en 10 minutos, mientras que en Viena alcanzamos una cobertura casi completa en unos 6-7 minutos. Además, debemos leer atentamente estos valores de tiempo: si nos encontramos en un incidente desafortunado, debemos llegar al dispositivo, recogerlo, regresar (haciendo que el tiempo de viaje duplique el tiempo de accesibilidad), instalarlo, etc. en una situación en la que cada minuto puede ser una cuestión de vida o muerte.
Entonces, desde una perspectiva de desarrollo, las conclusiones son garantizar que tengamos datos completos y luego usar los mapas de accesibilidad y población, combinarlos, analizarlos y aprovecharlos al implementar nuevos dispositivos y nuevas ubicaciones para maximizar la población efectiva alcanzada. .