Это продолжение серии статей об анализе данных персональных тренировок из набора 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. Наполнение и способы изображения можно комбинировать, но в то же время не стоит перегружать картографическое изображение одновременным отображением разных показателей.