Это продолжение серии статей об анализе данных персональных тренировок из набора FIT-файлов, которые создаются при использовании носимых устройств (фитнес-браслеты, часы, смартфоны, велокомпьютеры). В предыдущих двух я рассказывал о том, как получить доступ к данным совершенных тренировок и как можно визуализировать показатели на графике.
В этой статье я расскажу как визуализировать пространственные данные одной или нескольких тренировок одновременно на интерактивной карте, как отобразить значение показателя скорости на треке тренировки и как определить самое частое место тренировок за период времени.
Для решения задач я использовал несколько библиотек Python, предназначенных для работы с пространственными данными и их визуализации – GeoPandas, PyProj, Shapely, Folium, Branca, H3-Pandas, OSMnx.
Где взять координаты
Для таких уличных тренировок как бег, езда на велосипеде, плавание в открытой воде, катание на лыжах и др. носимые устройства при наличии датчика GPS и/или ГЛОНАСС и устойчивом сигнале с обнаруженных спутников сохраняют текущие координаты местоположения в момент времени наряду с другими показателями (пульс, скорость и т.п.)
В 95% случаев приемники, работающие с глобальными навигационными системами, обеспечивают точность определения местоположения до 15 метров.
Координаты сохраняются в сообщении Record FIT-файла для каждой точки записи в виде пары Latitude/Longitude. Производные от всех координат одной уличной активности сохраняются в сообщении Session – такие как местоположение старта активности (первая записанная точка из набора Record, start_position_lat/start_position_long) и координаты экстента всего трека тренировки (крайняя северо-восточная и крайняя юго-западная вершины минимального прямоугольника, описывающего весь трек активности nec_lat/nec_long и swc_lat/swc_long).
Последние используются в основном для фиксирования границы карты при визуализации тренировки. Для понимания приведу рисунок:
Для хранения координат Garmin и другие производители устройств используют 32-bit integer и оригинальные координаты выглядят подобным образом:
Для приведения координат к более привычному виду – географическим координатам (55.7558, 37.6173) – можно использовать следующую формулу:
# Coordinate conversion example
import psycopg2
import pandas as pd
activity_id = 100007015961
conn = psycopg2.connect(host="localhost", database="garmin_data", user="postgres", password="******")
df = pd.read_sql_query("""select * from session
where activity_id ={}""".format(activity_id), conn)
c_value = 1/((2**32) / 360)
df['start_latitude'] = df.start_position_lat*c_value
df['start_longitude'] = df.start_position_long*c_value
Я сделал аналогичные расчеты для координат из таблицы Record.
Как визуализировать трек тренировки на интерактивной карте
Для начала создадим датафрейм с данными одной активности, включая пары координат latitude/longitude:
import psycopg2
import pandas as pd
activity_id = 100007015961
conn = psycopg2.connect(host="localhost", database="garmin_data", user="postgres", password="*****")
df = pd.read_sql_query("""select * from record
where activity_id ={} order by timestamp asc""".format(activity_id), conn)
Далее мы будем использовать библиотеку GeoPandas, предназначенную для работы с пространственными данными в Python. По названию можно отметить, что библиотека сочетает в себе функции популярной Pandas и другой Shapely для обработки геометрии объектов.
Основной объект модуля GeoPandas – геодатафрейм, он полностью совпадает с определением датафрейма Pandas, но также включает в себя поле Geometry, в котором содержится определение пространственного объекта.
В векторном представлении все реальные объекты можно представить с помощью трех типов данных – точка, линия и полигон - в библиотеке Shapely Point, LineString и Polygon соответственно.
Трек активности (линия) представляет собой упорядоченный набор пар координат. Создадим новый геодатафрейм, где будет храниться только геометрия линейного объекта:
line = gpd.GeoDataFrame(crs = CRS('EPSG:4326'), geometry=[LineString(zip(df.longitude, df.latitude))])Параметр crs содержит сведения о координатной системе. При работе с данными, полученными от глобальных навигационных систем, будем использовать стандартную WGS-84 (ее код EPSG:4326). Используемый здесь модуль PyProj позволяет работать с координатными системами и их трансформацией.
Для отображения трека активности на карте необходимо подготовить координаты центроида и экстента линии. При использовании этих параметров, мы получим аккуратную вид с автоматически подобранным масштабом. Используем стандартные функции модуля GeoPandas – centroid и bounds для расчета координат центроида и его экстента:
line_centroid = line.centroid
line_bounds = line.boundsСледующим шагом настраиваем картографическую основу интерактивной карты и добавляем слой с треком активности. Импортируем библиотеку Folium, позволяющую создавать несколько типов карт Leaflet:
import folium
m = folium.Map([line_centroid.y, line_centroid.x], tiles='Stamen Terrain')
folium.GeoJson(line).add_to(m)
folium.FitBounds([[line_bounds.miny[0], line_bounds.minx[0]],[line_bounds.maxy[0], t_bounds.maxx[0]]]).add_to(m)
mДля команды Map первый параметр (location) указывает центроид отображаемого объекта, в качестве основы (tiles) в этом примере я использовал вариант Stamen Terrain. Другие доступные картографические основы можно посмотреть здесь.
В результате получаем следующий вид:
Как отобразить динамику показателя на карте
В следующих примерах я покажу как можно отобразить значение скорости на картографической основе. Первый способ – через точечные объекты, второй через линейный.
Точечный способ
Этот способ отображения я нашел здесь. Для начала подготовим исходные данные, аналогично предыдущему способу:
import psycopg2
import pandas as pd
import geopandas as gpd
from shapely.geometry import LineString
activity_id = 97199718403
conn = psycopg2.connect(host="localhost", database="garmin_data", user="postgres", password="*****")
df = pd.read_sql_query("""select * from record
where activity_id ={} order by timestamp asc""".format(activity_id), conn)Для расчета центроида и границ экстента создадим линейный объект, соответствующий треку тренировки:
line = gpd.GeoDataFrame(index=[0], crs = 'EPSG:4326', geometry=[LineString(zip(df.longitude, df.latitude))]) Для отображения показателя скорости точечным способом создадим еще один геодатафрейм, в который помимо координат добавим атрибут скорости (speed):
gdf = gpd.GeoDataFrame(df[['speed', 'latitude', 'longitude']], crs='EPSG:4326',geometry=gpd.points_from_xy(df.longitude, df.latitude))Создадим шкалу цветовых значений для параметра speed, учитывая, что минимальная возможная скорость 0 км/ч, максимальная – 40 км/ч (после анализа всей таблицы Record).
Импортируем библиотеку Branca и создаем шкалу цветовых значений на основе линейного градиента, в который войдут цвета ['cyan', 'green','yellow', 'orange', 'red'], соответствующие значениям скорости [0, 10, 20, 30, 40]:
import branca
colormap = branca.colormap.LinearColormap(colors=['cyan', 'green','yellow', 'orange', 'red'],
index=[0, 10, 20, 30, 40], vmin=0, vmax=40,
caption='Average Speed, km/h')Импортируем библиотеку Folium и добавляем строки с отображением точек через CircleMarker и подготовленную заранее шкалу цветовых значений и выводим результат.
import folium
m = folium.Map([line.centroid.y, line.centroid.x], tiles='cartodbpositron')
for point, alt in gdf[['geometry', 'speed']].values:
folium.CircleMarker([point.y, point.x], radius=6, fill=True, color='r', fill_color=colormap(alt), fill_opacity=1.0).add_to(m)
m.add_child(colormap)
folium.FitBounds([[line.bounds.miny[0], line.bounds.minx[0]],[line.bounds.maxy[0], line.bounds.maxx[0]]]).add_to(m)
m
Линейный способ
Этот способ отличается от предыдущего тем, что значение скорости отображается не для точки в момент измерения, а для сегмента (линии) между двумя соседними точками. Для значения скорости я буду использовать среднее значение скорости для двух соседних точек.
Подготавливаем исходные данные и объект line для расчета центроида и экстента карты:
import psycopg2
import pandas as pd
from pyproj import Proj, CRS, transform
import geopandas as gpd
from shapely.geometry import Point, LineString, Polygon
activity_id = 97199718403
conn = psycopg2.connect(host="localhost", database="garmin_data", user="postgres", password="*****")
df = pd.read_sql_query("""select * from record
where activity_id ={} order by timestamp asc""".format(activity_id), conn)
line = gpd.GeoDataFrame(index=[0], crs = 'EPSG:4326', geometry=[LineString(zip(df.longitude, df.latitude))]) Выберем нужные колонки из исходного датафрейма и создадим несколько новых колонок, в которых мы рассчитаем среднюю скорость сегмента avg_speed и его порядковый номер segment:
df = df[['record_id','speed', 'latitude', 'longitude']]
df['avg_speed'], df['segment'] = 0, 0Скопируем координаты предыдущей точки для каждой записи в упорядоченном датафрейме, используя функцию shift в новые колонки prev_lat/prev_long:
df['prev_lat'] = df.latitude.shift(1)
df['prev_long'] = df.longitude.shift(1)Получаем таблицу со следующей структурой:
Далее рассчитаем значение средней скорости как округленное и усредненное текущего и предыдущего значений:
for row in range(1, len(df)):
df.loc[row, 'avg_speed'] = round((df.loc[row-1, 'speed']+df.loc[row, 'speed'])/2)Рассчитаем порядковый номер сегмента. Здесь под сегментом я буду понимать отрезок или несколько отрезков, где среднее значение скорости постоянно. Удалим первую строку, так как вся информация о записи есть во второй:
for i in range(2, len(df)):
df.loc[i, 'segment'] = [df.loc[i-1, 'segment'] if df.loc[i, 'avg_speed'] == df.loc[i-1, 'avg_speed'] else df.loc[i-1, 'segment']+1]
df = df.iloc[1:,:]В итоге получаем таблицу следующего вида:
Создадим геометрию для точки от (from) и точки до (to) для каждой записи и соберем линейный геодатафрейм, состоящий из набора линий:
df['from'] = [Point(xy) for xy in zip(df.prev_long, df.prev_lat)]
df['to'] = [Point(xy) for xy in zip(df.longitude, df.latitude)]
df['polyline'] = df.apply(lambda row: LineString([row['from'], row['to']]), axis=1)
gdf = gpd.GeoDataFrame(df[['segment', 'avg_speed']], geometry=df.polyline, crs='EPSG:4326')Затем укрупним линии до отдельных сегментов, упрощая задачу для отображения отрезков с последовательными одинаковыми значениями. Используем для этого функцию dissolve:
segments = gdf.dissolve(by=['segment', 'avg_speed'])
segments = segments.reset_index()Импортируем Folium и Branca, добавляем аналогичную шкалу цветовых значений на базе линейного градиента:
import folium
import branca
colormap = branca.colormap.LinearColormap(colors=['cyan', 'green','yellow', 'orange', 'red'],
index=[0, 10, 20, 30, 40], vmin=0, vmax=40,
caption='Average Speed, km/h')
m = folium.Map([line.centroid.y, line.centroid.x], tiles='cartodbpositron')
folium.GeoJson(segments, style_function=lambda x:{'color': colormap(x['properties']['avg_speed'])}).add_to(m)
folium.FitBounds([[line.bounds.miny[0], line.bounds.minx[0]],[line.bounds.maxy[0], line.bounds.maxx[0]]]).add_to(m)
m.add_child(colormap)
m
Как найти самое частое место тренировок
H3 Uber's Hexagonal Hierarchical Spatial Index
Я вдохновился этой статьей перед тем как приступил к поиску решения на вопрос о самом частом месте тренировок.
В основе предложенного способа лежит метод разбиения поверхности земного шара на участки – гексагоны. Метод H3 Hexagonal Hierarchical Spatial Index разработан Uber для проведения детального анализа в установке динамических цен на разных уровнях исследования.
Существует готовая библиотека h3 для расчета гексагонов в Python, я буду использовать аналогичную H3-Pandas, так как она позволяет работать в том числе с мультиполигональными объектами (объект состоит из двух и более полигонов).
Данные OSM
В данном примере я буду искать частые места велотренировок в границах Москвы. Для начала подготовим полигон города, в пределах которого я буду проводить пространственный анализ.
Используем базу готовых слоев пространственных данных OpenStreetMap (OSM). Для доступа к данным импортируем библиотеку OSMnx:
import geopandas as gpd
import osmnx as ox
# gdf = ox.geocode_to_gdf('Москва', which_result=1) # Old Moscow polygon
gdf = ox.geocode_to_gdf('Москва', which_result=2) # Old Moscow + New Moscow polygon
import folium
m = folium.Map([gdf.centroid.y, gdf.centroid.x], tiles='cartodbpositron')
folium.GeoJson(gdf).add_to(m)
folium.FitBounds([[gdf.bounds.miny[0], gdf.bounds.minx[0]],[gdf.bounds.maxy[0], gdf.bounds.maxx[0]]]).add_to(m)
m
Отрисовка гексагонов h3
Заполним подготовленный полигон с текущими границами Москвы гексагонами h3. Для своего исследования я использовал уровень 9 (средняя площадь гексагона ~0,1 кв км). Он позволяет анализировать на уровне участков района города. Для более детального анализа (улиц, дорог) можно взять уровень 10 или 11:
import geopandas as gpd
import osmnx as ox
import h3pandas
gdf = ox.geocode_to_gdf('Москва', which_result=2) # Old Moscow + New Moscow polygon
moscow_hex = moscow.h3.polyfill_resample(9)
import folium
m = folium.Map([gdf.centroid.y, gdf.centroid.x], tiles='cartodbpositron')
folium.GeoJson(gdf).add_to(m)
folium.GeoJson(moscow_hex).add_to(m)
folium.FitBounds([[gdf.bounds.miny[0], gdf.bounds.minx[0]],[gdf.bounds.maxy[0], gdf.bounds.maxx[0]]]).add_to(m)
m
Для каждого из отрисованных гексагонов я хочу рассчитать количество раз, когда трек моей тренировки проходил через него. Для этого необходимо подготовить слой с треками всех тренировок.
Слой со всеми треками тренировок
В данном примере я анализировал треки моих велотренировок на шоссейном велосипеде, а также велопрогулок на МТБ и велосипедах городского велопроката, поэтому выбрал все уникальные номера активностей из таблицы Session с помощью следующего запроса:
select activity_id from session
where sport = 'cycling' and (sub_sport = 'generic' or sub_sport = 'road')
order by timestamp ascДалее выбираем все строки для указанных активностей из таблицы Record, и создаем трек для каждой из них, добавляя все в один геодатафрейм. Существует опция для сохранения слоя в формате geojson для того, чтобы не пересоздавать его каждый раз:
activities_list = activities['activity_id'].to_list()
gdf = gpd.GeoDataFrame()
for activity_id in activities_list:
df = pd.read_sql_query("""select activity_id, longitude, latitude from record
where activity_id ={} order by timestamp asc""".format(activity_id), conn)
line = gpd.GeoDataFrame(df[['activity_id']], index=[0], crs = CRS('EPSG:4326'), geometry=[LineString(zip(df.longitude, df.latitude))])
gdf = gdf.append(line)
gdf.to_file(filename='geodata/cycling_routes.geojson', driver='GeoJSON')
Spatial Join
Для расчета количества пересечений линий треков и полигонов гексагонов будем использовать функцию spatial join, позволяющую связывать данные на основе их положения в пространстве. Логика аналогична обыкновенному join. Рассчитаем количество раз через параметр агрегирования aggfunc='count':
sjoin = moscow_hex.sjoin(routes, how='left')
sjoin = sjoin[['h3_polyfill', 'activity_id']]
sjoin = pd.pivot_table(sjoin, values ='activity_id', index='h3_polyfill', aggfunc='count')
sjoin = sjoin.reset_index()
moscow_hex = moscow_hex.merge(sjoin, how='left', on=['h3_polyfill'])Получаем геодатафрейм со следующей структурой:
Исходя из распределения значений, добавим функцию для шкалы цветовых значений:
def colorcode(x):
if x == 0:
return 'whitesmoke'
elif x == 1:
return 'lightgray'
elif x in [2,3]:
return 'sandybrown'
elif x in range(4,21):
return 'chocolate'
else:
return 'darkred'Импортируем Folium и Branca и добавляем подготовленный слой с гексагонами с цветовой шкалой значений.
Визуально четко выделяются три наиболее посещаемых района: 1 – Крылатское, где проходят регулярные велотренировки на шоссейном велосипеде, 2 – район Краснопресненской набережной, где я часто пользовался услугами городского велопроката для того, чтобы добраться на работу, 3 – район Зюзино, где проживал продолжительное время.
Полный код доступен здесь.
Итоги
Я рассмотрел некоторые варинты визуализации данных тренировок на интерактивных картах исключительно с помощью библиотек Python. Наполнение и способы изображения можно комбинировать, но в то же время не стоит перегружать картографическое изображение одновременным отображением разных показателей.
