analyse_saison_2020_en.Rmd

---
title: "Exploratory data analysis of spatial sailing data in R"
subtitle: "What can we learn from 230 dinghy miles?<br/>(Alternative title: That's how a sailing data scientist spends his time in autumn!)"
author: Sandro Raabe
date: '`r format(Sys.Date(), "%dth %B, %Y")`'
output:
  prettydoc::html_pretty:
    theme: architect
---

<!-- 
md_document:
variant: markdown_github
-->


```{r setup, include=FALSE}
checkpoint::checkpoint("2020-08-01", checkpointLocation = "/")
knitr::opts_chunk$set(echo = FALSE, warning=FALSE, message=FALSE)
Sys.setlocale("LC_TIME", "English")
```


```{r titelbild, out.width = "100%", fig.align = "default"}
knitr::include_graphics("for_import/2019-05-23 20.59.38.jpg")
```

```{r read_data}
library(prettydoc)
library(XML)
library(dplyr)
library(readr)
library(tidyr)
library(ggplot2)
library(purrr)
library(stringr)
library(lubridate)
library(readxl)
library(glue)
library(leaflet)
library(ggmap)
library(highcharter)
library(yarrr)

gpx_to_tibble <- function(filename){
  filename %>%
    xmlTreeParse(useInternalNodes = TRUE) %>%
    xmlRoot %>%
    xmlToList %>%
    .$trk %>%
    (function(x) unlist(x[names(x) == "trkseg"], recursive = FALSE)) %>%
    map_dfr(function(x) as_tibble(t(unlist(x)))) %>% 
    select(-ele) %>% 
    set_names(str_remove, ".attrs.") %>% 
    mutate(across(c("lon", "lat"), as.numeric)) %>% 
    mutate(across(time, ymd_hms))
}

waterspeed <- map_dfr(list.files("for_import", pattern = "csv", full.names=T), read_csv) %>% 
  select(Lon = Long, Lat, Time = Date, Course, Speed = `Speed (m/s)`) %>%
  rlang::set_names(str_to_lower) %>% 
  mutate(across(speed, ~.x*1.4384))
komoot <- list.files("for_import", pattern = ".gpx$", full.names=T) %>% map_dfr(gpx_to_tibble)

metadata1 <- read_csv("for_import/saison2020_meta.csv") %>% mutate(across(Date, date)) %>% 
  set_names(c("Datum", "Distance", "Duration", "mean_speed", "max_speed", "fort",
              "wind_speed", "wind_angle", "fort2")) %>%
  select(-contains("fort")) %>% group_by(Datum) %>% 
  summarise(across(c(Distance, Duration), sum), 
            across(c(mean_speed, max_speed, wind_speed), mean),
            across(wind_angle, first),
            .groups = "drop") %>% 
  mutate(across(c(mean_speed, max_speed, Distance), ~.x/1.852),  # km to nm
         across(Duration, ~.x/60/60)) %>%  # s to h
  mutate(across(wind_angle, factor, levels = c("North", "North East", "East", "South East", "South", "South West", "West", "North West")))

metadata <- read_excel("for_import/Metadaten.xlsx") %>%
  mutate(across(Datum, dmy)) %>%
  group_by(Datum) %>% 
  arrange(desc(Segelpartner)) %>% # damit nicht "allein" gewinnt
  mutate(across(Segelpartner, str_split, " & ")) %>% 
  slice(1) %>% ungroup %>% 
  left_join(metadata1, "Datum")

all_dat <- bind_rows(waterspeed, komoot) %>% filter(year(time) == 2020) %>% 
  arrange(time) %>% distinct() %>% 
  mutate(Datum = date(time), Wochentag = factor(weekdays(Datum), levels(wday(today(), label=T,abbr = F, week_start = 1)))) %>% 
  left_join(metadata, "Datum")

if (any(is.na(all_dat$Boot))) stop("Boot fehlt manchmal. Metadaten.xlsx nachtragen!")
```

# Management Summary

This is a case study on **exploratory data analysis** of manually collected GPS data. Through various data wrangling packages (`XML`, `tidyverse`) and visualization tools  (`ggplot2`, `highcharter`, `leaflet`) we explore relationships in the data, which, however, are not further investigated through hypothesis tests.

# Introduction
In spring 2020 I joined a [local dinghy sailing club](http://www.osgev.de) on the Alster lake in Hamburg, Germany. This summer I used their boats extensively and built up my dinghy sailing skills. I recorded some of these trips with GPS and examined them with various visualization tools. These are some of my insights:

* Apparently I don't like sailing on **Wednesdays** and **Thursday** is Alster exploration day.
* With the boat type **Möwe** one apparently prefers to stay **close to the mooring**.
* **Kielzugvogel** should take part in **German Sailing League**
* **Corona** enhances **solo sailing skills**.
* GPS tracks of Regatta races look like **balls of wool**.
* The center of the Alster is (as expected) the **sailing hotspot**.

# Data overview

Since I did not record all sessions (of course), we start with an overview of the data used in this analysis:

```{r}
anz <- round(sum(metadata$Distance, na.rm=T) + 40)
```

* Period: `r str_flatten(format(range(date(all_dat$time)), '%A, %d.%m.%Y'), " to ")`
* Number of recorded **GPS points: `r format(nrow(all_dat), big.mark = ".")`**
* Number of recorded **days: `r n_distinct(all_dat$Datum)`**
* Total distance recorded: **`r anz` nautical miles**, which is approximately `r round(anz*1.8)` km
* Number of used **boat types: `r n_distinct(all_dat$Boot)`**
* Number of different **sailing partners: `r n_distinct(unlist(all_dat$Segelpartner))`**
* `r n_distinct(all_dat$Datum)` sessions in `r round(as.numeric(difftime(max(all_dat$Datum), min(all_dat$Datum), units="weeks")), 1)` weeks yields `r (n_distinct(all_dat$Datum)/as.numeric(difftime(max(all_dat$Datum), min(all_dat$Datum), units="weeks"))) %>% round(1)` sessions per average week, or **1 session each `r round(as.numeric(difftime(max(all_dat$Datum), min(all_dat$Datum), units="days"))/n_distinct(all_dat$Datum), 1)` days**

# Calendar View

The days in the data set can be displayed without aggregation:

```{r kalender, fig.height=2, fig.width=8, height=200}
dfr <- tibble(date = seq(floor_date(min(all_dat$Datum), "month"), 
                         ceiling_date(max(all_dat$Datum), 'month') - 1, "day")) %>% 
  mutate(comment = ifelse(date %in% all_dat$Datum, "ja", "nein"))
dfr$day <- factor(strftime(dfr$date,format="%a"), levels = rev(levels(wday(today(), label=T, week_start = 1))))
dfr$week <- factor(strftime(dfr$date,format="%V"))
dfr$month <- factor(strftime(dfr$date,format="%B"),levels = levels(month(today(), label=T, abbr = F)))
dfr$ddate <- factor(strftime(dfr$date,format="%d"))

ggplot(dfr,aes(x=week,y=day))+
  geom_tile(aes(fill=comment))+
  geom_text(aes(label=ddate))+
  scale_fill_manual(values=c("green","#ECEBEB"))+
  facet_grid(~month,scales="free",space="free")+
  labs(x="Week",y="")+
  theme_bw(base_size=10)+
  theme(legend.title=element_blank(),
        panel.grid=element_blank(),
        panel.border=element_blank(),
        axis.ticks=element_blank(),
        strip.background=element_blank(),
        legend.position="none",
        legend.justification="right",
        legend.direction="horizontal",
        legend.key.size=unit(0.3,"cm"),
        legend.spacing.x=unit(0.2,"cm"),
        plot.background = element_rect(fill = "transparent"),
        panel.background = element_rect(fill = "transparent")
  )
```

**Observation:** `r dfr %>% filter(comment=="ja") %>% count(month, sort=T) %>% slice_max(n, n=1) %>% mutate(text = glue("{month} was the most active sailing month ({n} training sessions).")) %>% pull(text)`

# Days of the week

Now we can split the dataset according to the contained variables, beginning with the weekday.

## Frequency distribution

```{r wochentage_histogramm}
add_colors <- function(df, colname){
  # make same color vector regardless of shape of levels
  set.seed(123)
  col <- tibble(sort(unique(as.character(df[[colname]])))) %>% 
    set_names(colname) %>%
    mutate(col = randomcoloR::distinctColorPalette(n(), runTsne = TRUE))
  if(is.factor(df[[colname]])) col[[colname]] <- factor(col[[colname]], levels = levels(df[[colname]]))
  df %>% left_join(col, colname)
}

make_hist <- function(fill_var, dat = all_dat){
  d <- all_dat %>% distinct_at(c("Datum", fill_var)) %>% set_names("Datum", "x") %>% 
    unnest(cols=c(x)) %>% 
    #mutate(across(x, ~reorder(.x, .x, function(x) -length(x)))) %>% 
    count(x)
  
  if(fill_var != "Wochentag") d <- d %>% arrange(desc(n))

  highchart() %>% hc_chart(type = "column") %>%
    hc_add_series(data = d$n, name = "number of sessions",
                  showInLegend = FALSE,
                  dataLabels = list(enabled = TRUE)) %>% 
    hc_title(text = glue("Number of sessions per {ifelse(fill_var == 'Wochentag', 'weekday', ifelse(fill_var == 'Boot', 'boat', fill_var))}")) %>% 
    hc_yAxis(text = NULL) %>% 
    hc_xAxis(text = NULL, categories = d$x) %>% 
    hc_colors("#9E9E9E") %>% 
    hc_plotOptions(column = list(colorByPoint = TRUE))
}
make_hist("Wochentag")
```

**Observation:** On Wednesdays, the big alster race *kangaroo regatta* seems to spoil my sailing, but I seem to like participating at the club's own Tuesday regatta. Alternatively: After our Tuesday's regatta, am I so tired that I rather stay home on Wednesdays?

## Spatial distribution

We use the `Leaflet` package for interactive visualisation of all GPS tracks (use mouse wheel or soft buttons for zooming, use legend for weekday switch):

```{r wochentage_karten, fig.height=6, fig.width=10}
library(ggmap)
alster_map <- get_stamenmap(c(left=9.99, right=10.02, bottom=53.555, top=53.583), zoom=15, maptype = "toner-lite")
draw_map <- function(col_var, nrow, colorize=FALSE){
  dat <- all_dat %>% unnest(cols=all_of(col_var))
  if(col_var != "Wochentag"){
    counts <- dat %>% distinct_at(c("Datum", all_of(col_var))) %>% count(get(col_var), sort=T) %>% set_names(c("x", "n"))
    dat <- dat %>% mutate(across(all_of(col_var), factor, levels = counts$x))
  }
  
  if(!colorize){
    dat$col <-   "#9E9E9E"
  }else{
    dat <- add_colors(dat, col_var)
  }
  
  # durch unnest werden Strecken doppelt bei 2 Segelpartnern gezeichnet -> ok!
  out <- ggmap(alster_map) + geom_point(data = dat,
                                        aes(lon,lat, color = I(col)), size=0.1, alpha=0.7) +
    theme(legend.position = "none", axis.ticks = element_blank(),
          axis.text = element_blank(), axis.title = element_blank()) +
    theme(
      plot.background = element_rect(fill = "transparent"),
      panel.background = element_rect(fill = "transparent")
    )
  if (col_var == "Wochentag"){
    out <- out + facet_wrap(col_var, nrow = nrow, labeller = as_labeller(setNames(unique(glue("Alle {t}e", t = all_dat[[col_var]])), unique(all_dat[[col_var]]))))
  }else{
    out <- out +
      facet_wrap(col_var, nrow = nrow)
  }
  out
}

draw_leaflet <- function(col_var){
  dat <- all_dat %>% unnest(cols=c(col_var))
  if(col_var != "Wochentag"){
    counts <- dat %>% distinct(across(c("Datum", all_of(col_var)))) %>% 
      count(get(col_var), sort=T) %>% set_names(c("x", "n"))
    dat <- dat %>% mutate(across(all_of(col_var), factor, levels = counts$x))
  }
  dat <- add_colors(dat, col_var) %>% arrange(across(c(all_of(col_var), "Datum")))
   
  all_polylines <- dat %>% group_by(across(all_of(col_var))) %>%
    summarize(across(c(lon, lat), list), .groups = "drop") %>% select(2:3)
  lons <- all_polylines$lon
  lats <- all_polylines$lat
  
  groups <- unique(dat[[col_var]])
  base <- leaflet() %>% addProviderTiles("CartoDB")  %>%  
    fitBounds(  lng1 =10.01, lng2 = 10.03,  lat1 = 53.555, lat2 = 53.58)
  
  for (i in seq_along(lons)){
    base <- base %>% addPolylines(lons[[i]], lats[[i]], color = "#9E9E9E", group = groups[i], weight = 1.5, opacity = 0.7)
  }
  base <- base %>% addLayersControl(baseGroups = groups, 
                                    position = "topleft",
                                    options = layersControlOptions(collapsed = FALSE))
  return(base)
}
draw_leaflet("Wochentag") #, nrow = 2)

```

**Observation:** On Saturdays I like to stay close to the mooring, whereas on Thursdays I sail all the way to the university. The regatta race tracks in the shape of triangles are clearly visible (using zoom).

# Headings and wind direction

For every GPS point we know the current heading (a.k.a. the "driving direction"). We can count these and visualise as histogram (frequency diagram) like a compass. Additionally, we draw the wind directions.

```{r kompassrose}
wind_grad <- c(North = 0, East = 90, South = 180, West = 270, "Nord East" = 45, 
               "South East" = 135, "South West" = 180+45, "North West" = 360-45) %>% 
  as.list() %>% as_tibble %>% t %>% as.data.frame %>% 
  tibble::rownames_to_column("wind_angle")

wind_data <- metadata %>% left_join(wind_grad, by = "wind_angle") %>% select(course = V1) %>% 
    mutate(group = "Wind Directions") %>% select(course, group) %>%
  map_df(., rep, 1200) %>% mutate(across(course, ~.x + rnorm(n(), 0, 12)))

all_dat %>% select(course) %>% mutate(group = "Headings") %>% 
  bind_rows(wind_data) %>% 
  mutate(bin = cut(course, breaks = seq(-1,361, 3))) %>% 
  filter(!is.na(bin)) %>% group_by(group) %>% count(bin) %>% 
  hchart("column", hcaes(x = bin, y = n, group = group)) %>% 
  hc_chart(polar = TRUE) %>% 
  hc_title(text = "Common headings and wind directions (in Degrees)") %>% 
  hc_yAxis(title = list(text = NULL), labels = list(enabled = F)) %>% 
  hc_xAxis(title = list(text = NULL), labels = list(enabled =F))

```

The heading around 15° (North-Northeast) and 165° (South-Southeast) seems to be rather popular with me. This is obvious considered the geographical shape of the Alster lake: With its slim North-South orientation one goes more "up and down" then "left and right". Furthermore, one can clearly see the main wind directions for Hamburgs geographical location: Southwest and Northwest.

**Observation:** The usual start from the OSG mooring is towards North-Northeast. Furthermore, we often have Southwesterly and Northwesterly wind, so exactly those headings should be rare.

# Boats

## Frequency

```{r histogramm_boote}
dat_boote <- all_dat %>% distinct(Datum, Boot) %>% count(Boot) %>% slice_max(n) 
make_hist("Boot")
```

With `r  dat_boote$n` sessions *`r dat_boote$Boot`* was my favourite boat. This has two reasons: On the one hand it's a very beginner-friendly boat (I only did my certification in autumn 2019 and have never touched a sail boat before) - the other reason you can find in the [analysis of sailing partners](#frequency-1) later in this article.

**Observation:** Conger boats are great for learning how to sail.

## Spatial distribution

```{r boote_karte, fig.height=6, fig.width=10}
draw_leaflet("Boot")#, nrow = 2)
```

**Observation:** *Conger and* *Kielzugvogel* can get you anywhere, using the *Möwe* you better stay close to the mooring.

## Which boat is the fastest?

An interesting question arises: Which boat performs best, and what is *perfomance*? Is it top speed, regardless of wind conditions? Is it speed compared to wind speed? To subtract the influence of different wind speeds from of the single sessions, we re-calculate the measured boat speed in comparison to the wind speed and thus make the speed dimensionless as **speed in % of wind speed** - 100% means that we went *as fast as the wind*, at 50% *half as fast* et cetera.


```{r pirateplot, fig.height=5.5, fig.width=12}
levels <- all_dat %>% group_by(Boot) %>% 
  summarize(meanspeed = mean(speed/wind_speed*100, na.rm=T), .groups = "drop") %>% 
  arrange(desc(meanspeed)) %>% filter(!is.na(meanspeed)) %>% pull(Boot)
dat <- all_dat %>% filter(Boot %in% levels) 
dat$Boat <- factor(dat$Boot, levels = levels) 

pirateplot(formula = speed/wind_speed*100 ~ Boat,# point.col = "#6B6B6B",
           data = dat, inf.f.o = 0, point.o = .005) 
title(main = "Boat speed relativ to Wind Speed")
```

Interesting to note: Small *Conger* was almost twice as fast as the wind at its top speed, but the racing dinghies *470er* and *Laser* were only half as fast as the wind on average. These only perform well at high wind speeds, other don't need as much wind, especially if they are lightweight and have large sails (which is the case with the two winners).

**Beobachtung:** Kielzugvogel and J70 (the boat used in Bundesliga/German Sailing League) get most power out of any given wind - Kielzugvogel should probably participate at Bundesliga races?

# Sailing Partners

## Frequency

```{r treemap}
lvl_opts <-  list(
  list(
    level = 1,
    color = "yellow",
    borderWidth = 5,
    borderColor = "white",
    dataLabels = list(
      enabled = TRUE,
      align = "left",
      verticalAlign = "top",
      style = list(fontSize = "14px", textOutline = FALSE, fontWeight = "regular")
    )
  )
)

dat <- all_dat %>% unnest(cols = c(Segelpartner)) %>% 
  distinct(Datum, Segelpartner) %>% 
  mutate(n=1) %>% 
  data_to_hierarchical(group_vars = Segelpartner, size_var = n) %>% 
  bind_rows() %>% 
  add_colors("name") %>%
  select(-color) %>%
  rename(color = col) %>%
  as.list %>% transpose()

hchart(dat, type = "treemap", layoutAlgorithm = "strip",
       levels = lvl_opts,
       tooltip = list(valueDecimals = 0, valueSuffix = " sessions")) %>% 
  hc_title(text = "Size corresponds to number of sessions") %>% 
  hc_size(height = 400)
```

**Observation:** With `r bind_rows(dat) %>% slice_max(value) %>% pull(value)` sessions I was most often alone. This was most importantly due to the COVID-19 restrictions in May and June 2020 - I had no chance but to learn one-handed sailing.

## Spatial distribution of sailing partners

We use a static (i.e. not interactive) visualisation of the GPS tracks:

```{r segelpartner_ggmap, fig.height=10, fig.width=10}
draw_map("Segelpartner", nrow = 3, colorize = TRUE)
```

**Observation:** Clearly visible are racing partners Christoph, Bernd and Jochem with their tracks shaped like balls of wool.

# Favourite Alster Regions

We consider the two-dimensional density function of all driven tracks. In simple words: We divide the Alster lake into small rectangles and count how often we cross each rectangle during our sailing sessions. Then we can colour the rectangles green to red according to their frequency - like a COVID-19 hotspot map!

```{r dichteplot}
ggmap(alster_map) + stat_density2d(aes(x=lon, y=lat,fill=..level.., alpha=..level..), data=all_dat, geom="polygon")+ 
  scale_fill_gradient(low = "green", high = "red") +
  theme(legend.position = "none", axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        axis.title = element_blank(),
        plot.background = element_rect(fill = "transparent"),
        panel.background = element_rect(fill = "transparent"))
```

**Observations:** The red hotspots could be the club's Tuesday race, the mooring on the lower left is clearly visible.

# Further Research

This is just a small, visualisation-driven exploratory data analysis of the 2020 season, mostly univariate and without test of the found hypotheses. Interesting questions arose:

* Can I see the correlation between boat length and boat speed in my data?
* Can I leverage the *yardstick* system for a further boat vs. speed research?
* Besides the `r anz` nautical miles on the Alster lake I also completed 280 nm on the Baltic sea in 2020. We didn't include these here, but that would be an interesting analysis, too.

These questions could easily be examined using hypothesis tests or other machine learning methods, which would be out of scope here. The computer chips in the basement are, however, ready and I still need something to do in 2021 :-)

# Credits

Thanks [OSG](https://osgev.de/) for the great community, the crazy boats and the fun things we do together!

# The technical stuff

Data collection was carried out using [Komoot](https://www.komoot.com) and [Waterspeed](https://www.waterspeedapp.com/) mobile apps, this analysis was conducted using [R 4.0](https://www.r-project.org/) and the following useful helpers:

* Data Input:
  * `readr`
  * `readxl`

* Data Wrangling:
  * `dplyr`
  * `purrr`
  * `tidyr`
  * `lubridate`
  * `glue`

* Graphics:
  * `highcharter`
  * `ggplot2`
  * `randomcoloR`
  * `yarrr`

* Spatial Analysis:
  * [`Leaflet`](https://CRAN.R-project.org/package=leaflet)
  * [`ggmap`](http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf)
  * `XML`

* Output:
  * `rmarkdown`
  * `knitr`
  * `prettydoc`

The Code for calculations and visualisations can be downloaded in my Github Repository: [https://github.com/shosaco/sailing_analyses](https://github.com/shosaco/sailing_analyses), this page is available at [https://shosaco.github.io/sailing_analyses](https://shosaco.github.io/sailing_analyses).


```{r img.show = "hold", out.width = "50%", fig.align = "default"}
knitr::include_graphics("for_import/2020-05-19 18.09.50.jpg")
knitr::include_graphics("for_import/WhatsApp Image 2020-09-22 at 21.48.55.jpeg")
```