Это продолжение серии статей об анализе данных персональных тренировок из набора 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).
Последние используются в основном для фиксирования границы карты при визуализации тренировки. Для понимания приведу рисунок:
![Производные координаты для одной активности, которые сохраняются в сообщении Session: координаты старта активности и экстента трека Производные координаты для одной активности, которые сохраняются в сообщении Session: координаты старта активности и экстента трека](https://habrastorage.org/getpro/habr/upload_files/b6b/7fa/9ab/b6b7fa9ab5bea6cca3644a923b96b893.png)
Для хранения координат Garmin и другие производители устройств используют 32-bit integer и оригинальные координаты выглядят подобным образом:
![Пример хранения записей пары координат в таблице Session FIT-файла Пример хранения записей пары координат в таблице Session FIT-файла](https://habrastorage.org/getpro/habr/upload_files/22f/63b/4e7/22f63b4e78ce5162219c326ab25d6a63.png)
Для приведения координат к более привычному виду – географическим координатам (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), приведенных к привычному виду Пример координат (latitude/longitude), приведенных к привычному виду](https://habrastorage.org/getpro/habr/upload_files/29a/aed/0f3/29aaed0f3e19faef296b025a470bb373.png)
Как визуализировать трек тренировки на интерактивной карте
Для начала создадим датафрейм с данными одной активности, включая пары координат 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)
![Структура данных из таблицы Record, включая пары координат latitude/longitude Структура данных из таблицы Record, включая пары координат latitude/longitude](https://habrastorage.org/getpro/habr/upload_files/dd2/773/f7c/dd2773f7cffd392fa0ea8b3147c13b83.png)
Далее мы будем использовать библиотеку 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 позволяет работать с координатными системами и их трансформацией.
![Пример структуры геодатафрейма, где содержатся только координаты линейного объекта Пример структуры геодатафрейма, где содержатся только координаты линейного объекта](https://habrastorage.org/getpro/habr/upload_files/9b9/298/666/9b9298666c3521188dc9a7822baed902.png)
Для отображения трека активности на карте необходимо подготовить координаты центроида и экстента линии. При использовании этих параметров, мы получим аккуратную вид с автоматически подобранным масштабом. Используем стандартные функции модуля 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. Другие доступные картографические основы можно посмотреть здесь.
В результате получаем следующий вид:
![Визуализация трека велотренировки (синия линия) с помощью Folium на картографической основе Stamen Terrain, сфокусированной относительно центроида и по границам экстента трека Визуализация трека велотренировки (синия линия) с помощью Folium на картографической основе Stamen Terrain, сфокусированной относительно центроида и по границам экстента трека](https://habrastorage.org/getpro/habr/upload_files/8ab/f87/253/8abf87253ae766769525fd1db6922176.png)
Как отобразить динамику показателя на карте
В следующих примерах я покажу как можно отобразить значение скорости на картографической основе. Первый способ – через точечные объекты, второй через линейный.
Точечный способ
Этот способ отображения я нашел здесь. Для начала подготовим исходные данные, аналогично предыдущему способу:
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
![Визуализация показателя скорости точечным способом с помощью Folium на картографической основе cartodbpositron, сфокусированной относительно центроида и по границам экстента Визуализация показателя скорости точечным способом с помощью Folium на картографической основе cartodbpositron, сфокусированной относительно центроида и по границам экстента](https://habrastorage.org/getpro/habr/upload_files/2e2/6b1/6d6/2e26b16d65ab20c0563507ed2dd87af4.png)
Линейный способ
Этот способ отличается от предыдущего тем, что значение скорости отображается не для точки в момент измерения, а для сегмента (линии) между двумя соседними точками. Для значения скорости я буду использовать среднее значение скорости для двух соседних точек.
Подготавливаем исходные данные и объект 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)
Получаем таблицу со следующей структурой:
![Таблица Record с подставленными координатами предыдущей точки. Для первой строчки координат быть не должно Таблица Record с подставленными координатами предыдущей точки. Для первой строчки координат быть не должно](https://habrastorage.org/getpro/habr/upload_files/082/717/2ab/0827172ab6199b4c0665407365cf2962.png)
Далее рассчитаем значение средней скорости как округленное и усредненное текущего и предыдущего значений:
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:,:]
В итоге получаем таблицу следующего вида:
![Рассчитанные сегменты Рассчитанные сегменты](https://habrastorage.org/getpro/habr/upload_files/6ea/631/b3d/6ea631b3da1b65559f6bc87d10be6da4.png)
Создадим геометрию для точки от (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
![Визуализация показателя скорости линейным способом с помощью Folium на картографической основе cartodbpositron, сфокусированной относительно центроида и по границам экстента Визуализация показателя скорости линейным способом с помощью Folium на картографической основе cartodbpositron, сфокусированной относительно центроида и по границам экстента](https://habrastorage.org/getpro/habr/upload_files/7cd/49b/181/7cd49b1817352e7af0fe548d60e7ecfc.png)
Как найти самое частое место тренировок
H3 Uber's Hexagonal Hierarchical Spatial Index
Я вдохновился этой статьей перед тем как приступил к поиску решения на вопрос о самом частом месте тренировок.
В основе предложенного способа лежит метод разбиения поверхности земного шара на участки – гексагоны. Метод H3 Hexagonal Hierarchical Spatial Index разработан Uber для проведения детального анализа в установке динамических цен на разных уровнях исследования.
![H3 позволяет пользователям разбивать участки земного шара на гексагоны разного уровня для детального анализа. Источник:https://eng.uber.com/h3/ H3 позволяет пользователям разбивать участки земного шара на гексагоны разного уровня для детального анализа. Источник:https://eng.uber.com/h3/](https://habrastorage.org/getpro/habr/upload_files/a2d/16c/d6a/a2d16cd6a94c4dedec23688ce1ff7fa9.png)
Существует готовая библиотека 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
![В зависимости от порядкового номера результата поиска можно получить разные полигоны для Москвы. Слева: текущие административные границы города (which_result=2), справа: границы старой Москвы (which_result=1) В зависимости от порядкового номера результата поиска можно получить разные полигоны для Москвы. Слева: текущие административные границы города (which_result=2), справа: границы старой Москвы (which_result=1)](https://habrastorage.org/getpro/habr/upload_files/9d2/79d/871/9d279d87176dd88f3d8feec50903fc8e.png)
Отрисовка гексагонов 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
![В таком виде получается заполненный гексагонами h3 уровня 9 полигон для Москвы В таком виде получается заполненный гексагонами h3 уровня 9 полигон для Москвы](https://habrastorage.org/getpro/habr/upload_files/815/5d2/8a0/8155d28a00c17d7837d33dd62da3b6ac.png)
Для каждого из отрисованных гексагонов я хочу рассчитать количество раз, когда трек моей тренировки проходил через него. Для этого необходимо подготовить слой с треками всех тренировок.
Слой со всеми треками тренировок
В данном примере я анализировал треки моих велотренировок на шоссейном велосипеде, а также велопрогулок на МТБ и велосипедах городского велопроката, поэтому выбрал все уникальные номера активностей из таблицы 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')
![Подготовленный слой с треками всех активностей, связанных с велосипедом в пределах Москвы Подготовленный слой с треками всех активностей, связанных с велосипедом в пределах Москвы](https://habrastorage.org/getpro/habr/upload_files/510/c7b/820/510c7b8201c9185abfad2a432a01cb5f.png)
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'])
Получаем геодатафрейм со следующей структурой:
![Геодатафрейм с гексагонами и рассчитанным количеством раз прохождения треков тренировок в поле activity_id Геодатафрейм с гексагонами и рассчитанным количеством раз прохождения треков тренировок в поле activity_id](https://habrastorage.org/getpro/habr/upload_files/2b1/58c/55c/2b158c55cf7e849644bf1ce82edabc05.png)
Исходя из распределения значений, добавим функцию для шкалы цветовых значений:
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 и добавляем подготовленный слой с гексагонами с цветовой шкалой значений.
![Итоговая карта с наиболее часто посещенными территориями во время активностей связанных с велосипедом в пределах Москвы Итоговая карта с наиболее часто посещенными территориями во время активностей связанных с велосипедом в пределах Москвы](https://habrastorage.org/getpro/habr/upload_files/8d1/2d5/1ce/8d12d51ce56a4d4bb5218874243ad6c3.png)
![](https://habrastorage.org/getpro/habr/upload_files/07a/17d/b89/07a17db89cd31f559277113abf3504e9.png)
Визуально четко выделяются три наиболее посещаемых района: 1 – Крылатское, где проходят регулярные велотренировки на шоссейном велосипеде, 2 – район Краснопресненской набережной, где я часто пользовался услугами городского велопроката для того, чтобы добраться на работу, 3 – район Зюзино, где проживал продолжительное время.
Полный код доступен здесь.
Итоги
Я рассмотрел некоторые варинты визуализации данных тренировок на интерактивных картах исключительно с помощью библиотек Python. Наполнение и способы изображения можно комбинировать, но в то же время не стоит перегружать картографическое изображение одновременным отображением разных показателей.