フルマラソンの完走タイム

あまりの運動不足から最近ランニングをはじめ、継続するための目標としてフルマラソンを走ってみたいと思っている。 フルマラソンをとても走れるとは思えなかったが、調べると意外と走れるようである。 走れるといってもそう簡単にできるものではないので、下調べをしていると 自分がどの程度で走れるのか、また、ほかの人がどの程度で走っているのかに興味が出てきた。

スポロクリザルトというサイトで一部マラソン大会の結果が公表されていたので、 一般市民が参加しそうで比較的有名そうな金沢マラソンの記録をもとに、マラソン大会出場者の完走タイムについて簡単に調べてみた。

# パッケージの読込
library(tidyverse)
library(rvest)
library(lubridate)
# データの取得
## 男子総合の結果を順次取得する
## 時間もかかるため、ここでは読み込んだ既に取得したデータを読み込むことにする
# # 男子
# adress = "https://sporoku.jp/result/kanazawa_20191027/list?race_id=696&div_id=584&page="
# pages = 1:209
# 
# 
# dt_male = list()
# for(i in pages){
#   dt_male[i] =
#   read_html(paste0(adress, pages[i])) %>%
#     html_table()
# }
# saveRDS(dt_male,"kanazawa_marathon_male.rds")
# 
# # 女子
# adress = "https://sporoku.jp/result/kanazawa_20191027/list?race_id=696&div_id=585&page="
# pages = 1:67
# 
# dt_female = list()
# for(i in pages){
#   dt_female[i] =
#   read_html(paste0(adress, pages[i])) %>%
#     html_table()
# }
# saveRDS(dt_female,"kanazawa_marathon_female.rds")


dt_male = 
  readRDS("kanazawa_marathon_male.rds")
dt_female =
  readRDS("kanazawa_marathon_female.rds")

data = 
  map(list(dt_male, dt_female), ~ do.call("rbind",.x)) %>% 
  map(~ .x[,1:4] %>% 
        as_tibble() %>% 
        mutate(time = lubridate::hms(タイム),
         ranking = str_replace(順位,"位","") %>% as.integer())) %>% 
  map2(c("male", "female"), ~ mutate(.x, gender = .y))

data =
  do.call("rbind",data)
data =
  group_by(data, gender) %>% 
  mutate(ranking_percentile = ranking/max(ranking)) %>% 
  ungroup()
# 順位パーセンタイル値と完走時間の関係をプロット
data %>% 
  ggplot(aes(ranking_percentile, time, color = gender))+
  geom_point()+
  scale_y_time()

f:id:jerrarrdan:20200101154447p:plain

4人に1人はサブ4を達成していることが分かる。 このままでは詳細が見ずらいので、順位のパーセンタイル値ごとに完走時刻をプロットする。

data %>% 
  mutate(percentile_factor = cut(ranking_percentile,breaks = seq(0,1,by=.1))) %>% 
  ggplot(aes(ranking_percentile, time, color=gender))+
  geom_point()+
  scale_y_time() +
  scale_x_continuous(breaks=seq(0,1,by=.1))+
  facet_wrap(~percentile_factor, scales = "free")

f:id:jerrarrdan:20200101154513p:plain

サブ3.5を達成すると上位10%に入ることができるようである。
また、サブ4, サブ5近辺でグラフの傾きが少し変わっていることが確認できる。 これは、それぞれの時間内で完走を目指すためだろう。

## 完走時間のパーセンタイル値
data %>% 
  select(time, gender) %>% 
  group_split(gender) %>%
  map(~ 
        pull(.x,time) %>% 
        as.numeric() %>% 
        quantile(seq(0,1,by=.05)) %>% 
        data_frame(percentage = names(.),value=.) %>% 
        mutate(h = floor(value/3600),
               m = floor((value-h*3600)/60),
               s = floor(value-h*3600-m*60)) %>% 
        mutate(time = paste(h,m,s, sep = ":") %>% hms()) %>% 
        select(percentage, time))  %>% 
  map2(.y=c("female","male"),~mutate(.x,gender=.y)) %>% 
  do.call(what="rbind",args=.) %>% 
  spread(gender, time, convert=T) %>% 
  mutate(num=str_extract(percentage,"\\d+") %>% as.integer()) %>% arrange(num) %>% select(-num) %>% 
  mutate_at(vars(ends_with("male")), hms) %>% 
  mutate(diff_time = as.duration(female - male)) %>% 
  knitr::kable(format="markdown")
percentage female male diff_time
0% 2H 34M 52S 2H 11M 36S 1396s (~23.27 minutes)
5% 3H 44M 25S 3H 14M 22S 1803s (~30.05 minutes)
10% 4H 0M 18S 3H 30M 22S 1796s (~29.93 minutes)
15% 4H 15M 23S 3H 44M 18S 1865s (~31.08 minutes)
20% 4H 28M 6S 3H 54M 24S 2022s (~33.7 minutes)
25% 4H 37M 35S 4H 1M 3S 2192s (~36.53 minutes)
30% 4H 46M 56S 4H 12M 57S 2039s (~33.98 minutes)
35% 4H 55M 21S 4H 22M 52S 1949s (~32.48 minutes)
40% 5H 1M 54S 4H 29M 53S 1921s (~32.02 minutes)
45% 5H 13M 17S 4H 39M 50S 2007s (~33.45 minutes)
50% 5H 22M 27S 4H 48M 48S 2019s (~33.65 minutes)
55% 5H 30M 23S 4H 56M 20S 2043s (~34.05 minutes)
60% 5H 39M 5S 5H 4M 27S 2078s (~34.63 minutes)
65% 5H 47M 47S 5H 15M 21S 1946s (~32.43 minutes)
70% 5H 56M 13S 5H 26M 49S 1764s (~29.4 minutes)
75% 6H 6M 0S 5H 38M 33S 1647s (~27.45 minutes)
80% 6H 16M 25S 5H 49M 52S 1593s (~26.55 minutes)
85% 6H 28M 3S 6H 2M 59S 1504s (~25.07 minutes)
90% 6H 37M 51S 6H 19M 38S 1093s (~18.22 minutes)
95% 6H 49M 7S 6H 39M 38S 569s (~9.48 minutes)
100% 7H 0M 19S 7H 1M 55S 96s (~1.6 minutes)

上記の結果から、男子のタイムに30分を加えたタイムが、同じ走力を持つ女子と考えられる。

mapviewで躓いたので備忘録

mapviewの躓いた点の備忘録

久しぶりにmapviewパッケージをを触ったら躓いたので、メモ。

利用するライブラリ

library(tmap)
library(leaflet)
library(mapview)

利用するデータ

data("franconia")
class(franconia)
[1] "sf"         "data.frame"

描画対象物の色を白系にした場合に描画できない問題

描画対象物を白色を指定した場合に、地図が表示できなくなる場合は、以下のようにmap.types=にベースとなるタイルを指定してやるのが良い。 map.typesにはサイトから好みの名前を指定する。

mapview(franconia, map.types="OpenStreetMap.Mapnik", col.regions="snow")

tmapを利用した背景地図付きの地図の保存方法

下記だと描画モード(Rstudio上)では表示されるが、画像ファイルとして保存する際には背景画像が消えてしまう。

tmap_mode("view")
m =
  tm_shape(franconia)+
  tm_polygons()
m
# 下記では背景地図が保存されない
tmap_save(m, "test.png")

このように背景地図が描画されない f:id:jerrarrdan:20190404214145p:plain

解決方法としてはmapshot()関するを利用する。

mapshot(tmap_leaflet(m), file="test2.png")

f:id:jerrarrdan:20190404214321p:plain 表示できた!!!

nestしたデータフレームの扱い

nestしたデータフレームを扱い方(メモ用)

library(tidyverse)
iris %>% 
  group_by(Species) %>% 
  nest() %>%
  mutate(data = map2(data, Species, ~mutate(.x, Species_index = as.numeric(.y)))) %>%
  mutate(id = seq(1,3)) %>% 
  mutate(data = map2(data, id, ~slice(.x, seq.int(.y)))) %>% 
  mutate(MeanSepalLength = map_dbl(data, ~pull(.x, Sepal.Length) %>% mean))

上記の結果以下のデータを得ることができる。

  Species    data                id MeanSepalLength
  <fct>      <list>           <int>           <dbl>
1 setosa     <tibble [1 x 5]>     1            5.1 
2 versicolor <tibble [2 x 5]>     2            6.7 
3 virginica  <tibble [3 x 5]>     3            6.40

ルート探索

# 必要なライブラリのロード
library(tidyverse)
library(sf)
library(mapview)
library(stplanr)
library(OpenStreetMap)

# 2地点間のポイントデータ作成
sfd = 
  tribble(~lat,~lon,~location,
          36.075411, 136.212213,"福井大学",
          36.062244, 136.223214, "福井駅") %>% 
    st_as_sf(coords = c("lon", "lat"), crs = 6668)
# linestringへ変換
sf_line = 
  sfd %>% 
  st_combine() %>% 
  st_cast("LINESTRING") %>% 
  st_sf(a = 1)

# ラインからルートデータへの変換
line2route(sf_line, "route_osrm") %>% 
  mapview()

f:id:jerrarrdan:20181113225205p:plain

library(ggspatial)
sf_route = line2route(sf_line, "route_osrm")
ggplot(sf_route) +
  annotation_map_tile(zoomin = 0) +
  geom_sf()

f:id:jerrarrdan:20181113224409p:plain

upperLeft = st_bbox(sf_route) %>% as.vector() %>% .[c(2, 1)]
lowerRight = st_bbox(sf_route) %>% as.vector() %>% .[c(4, 3)]
osmap = OpenStreetMap::openmap(upperLeft = upperLeft,
                               lowerRight = lowerRight,
                               type = "stamen-toner",
                               mergeTiles = T) 
osmaplocation = OpenStreetMap::openproj(osmap, st_crs(sfd)$proj4string)
autoplot(osmaplocation) + 
  geom_path(data = 
              st_coordinates(sf_route)[,1:2] %>% 
              as.data.frame() %>%
              `colnames<-`(c("x", "y")),
            color="blue", size = 3)

f:id:jerrarrdan:20181113224431p:plain

GeoPackageが便利

GISで利用するデータのファイル形式について、 これまでいろいろな不満がありながらも、 そのほかの代替手段を知らずやむを得ずシェープファイルを利用していた。

たまたま、sfパッケージの入出力 について調べていたら、GeoPackageというファイル形式があること知った。

少し調べたところ、最近普及し始めているようなのでメモ。

参考としたWebSite

GeoPackageでポータブルOSMは約2年前(2016年)に書かれた記事である。 シェープファイルを利用していて常に悩まされていることが書かれており、 大体のファイル形式としてGeoPackageが挙げられている。

以下のような悩みが解決されるらしい。

  • 文字コードによる文字化け
  • ファイルが1つ(.gpkg)だけ
  • 列名の文字数

既にQGISArcGIS、GDAL/OGRで利用できるらしい。 また、QGIS3.0からは保存の際のデフォルトのファイル形式になるらしい。 ということで、全然動向を知らなかったが、GeoPackageに移行する時期に 既になっているらしいので、どんどん使っていこうと思う。

sf vignette 6

Rのsfパッケージの vignette6について、 これまでのvignetteで書けなかったことなどを雑多に書いてあるので、 あまりまとまりがない感じである。

EPSGとは

EPSGは空間座標参照システムとしてよく知られており、 IOGP (International Association of Oil & Gass Producers)によってメンテナンスされている。 下記の関数を利用することでEPSG一覧を見ることができる。

rgdal::make_EPSG()

sfで第二のgeometry columnを扱うには

sfオブジェクトは一つ以上のgeometry list-columnを持つことができるが、有効となるのは一つだけであり、 st_geometry()関数によって取得できる。 printメソッドを利用することでgeometry columnが有効か確認することができる。

library(sf)
demo(nc, ask =F, echo = F)
nc$geom2 = st_centroid(st_geometry(nc))
print(nc, n = 2)
Simple feature collection with 100 features and 14 fields
Attribute-geometry relationship: 0 constant, 8 aggregate, 6 identity, 1 NA's
Active geometry column: geom
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
epsg (SRID):    4267
proj4string:    +proj=longlat +datum=NAD27 +no_defs
First 2 features:
   AREA PERIMETER CNTY_ CNTY_ID      NAME  FIPS FIPSNO CRESS_ID BIR74
1 0.114     1.442  1825    1825      Ashe 37009  37009        5  1091
2 0.061     1.231  1827    1827 Alleghany 37005  37005        3   487
  SID74 NWBIR74 BIR79 SID79 NWBIR79                           geom
1     1      10  1364     0      19 MULTIPOLYGON (((-81.47276 3...
2     0      10   542     3      12 MULTIPOLYGON (((-81.23989 3...
                       geom2
1  POINT (-81.49826 36.4314)
2 POINT (-81.12515 36.49101)

geometry columnを変更するにはst_geometry<-またはst_set_geometry関数を利用する。

plot(st_geometry(nc))

f:id:jerrarrdan:20180224104751p:plain

st_geometry(nc) <- "geom2"
plot(st_geometry(nc))

f:id:jerrarrdan:20180224104912p:plain

st_simplifyでのトポロジー関係の維持について

st_simpligyトポロジー保持する関数であるが、これは個々のfeature geometrieで行われる。 つまり、ポリゴンに対して実行したものはポリゴンとなる。 しかし、2つのfeaturesが境界を共有しているときにst_simplifyを実行した場合には2つのポリゴンが同じ境界を共有しているとは限らない。これは個々のポリゴンごとに単純化が行われるためである。

sfオブジェクトに対してdplyr verbsが使えない

多くの開発者はパッケージをロードせずに、sf::のように関数の属するパッケージを明示して下記のように記述する。

i = sf::st_intersects(sf1, sf2)

しかし、sfオブジェクトに対して、dplyr::selectは利用できません。 このケースではsfパッケージをロードする必要がある。

表示されるerror/warning/messageについて

座標系は緯度経度であるが xxxは平面を仮定しています(although coordinates are longitude/latitude, xxx assumes that they are planar)

sfで利用されるgeometryに関する計算においてはGEOSライブラリを利用している。このライブラリは2次元のフラットなユクリッド距離の空間座標を考慮する。しかし、例えば極を含む緯度経度データを扱い場合にはこれだけでは不十分です。

下記に例を示すが、これは間違った結果を返すこととなる。

polygon = st_sfc(
  st_polygon(
    list(rbind(c(0, 80), c(120, 80),
               c(240, 80), c(0, 80)))), 
  crs =4326)
pole = st_sfc(st_point(c(0, 90)), crs = 4326)
st_intersects(polygon, pole)
Sparse geometry binary predicate list of length 1, where the predicate was `intersects'
 1: (empty)

st_centroidは緯度経度データに対しては正しい値を返しません(

st_centroid does not give correct centroids for longitude/latitude data)

先ほどの例と同じですが、中心は2次元空間を前提としています。

st_centroid(polygon)[[1]]

この場合中心は極である必要があります。

距離は10進数を前提としています(

dist is assumed to be in decimal degrees (arc_degrees).)

このメッセージはsfでは距離の値が度で与えられていることを示しています。このメッセージを避けるには正しいユニット(単位)を与えてやる必要があります。

pt = st_sfc(st_point(c(0,0)), crs = 4326)
buf = st_buffer(polygon, 1)
buf = st_buffer(polygon, units::set_units(1, degree))

両方の式のオペランドがunitオブジェクトもしくはアトミックベクトルである必要があります(

both operands of the expression should be "units" objects or _Error: $ operator is invalid for atomic vectors")

このメッセージはsfパッケージをロードした後にspatstatパッケージをロードした時に現れる。パッケージ単位で作成された単位オブジェクトの場合、unitクラスも実装するspatstatのメソッドが選択される。この問題はspatstatの後にunitパッケージをロードすることで解決できる。なお、バージョン1.53-2.012のspatstatパッケージではこれが解決されています。

Rではじめる機械学習

パラパラと読んだ。 表紙に「人気セミナー講師によるクリアな解説」とあるとおり、 様々な機械学習手法に関してわかりやすくまとめられていた。 ベイズ推定やディープラーニングにも関して述べられており、 とりあえず一通り網羅しているという感じであった。

所々詳しく学びたい方にお勧めの参考書が記載されており、 つかみどころだけを紹介(記載)しているスタイルが徹底されていてよいと思う。 また、ほとんど数式等は現れず、 これからRによる機械学習を始めたいという人にとってわかりやすい書籍だと思う。