---
title: 'Tutorial on anomaly detection, Part R: hurricanes'
author: "Pavlo Mozharovskyi"
output:
  pdf_document: default
  html_document:
    df_print: paged
---

This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook on the tutorial 
on "Anomaly detection" given on Wednesday the 13th of April 2022.

# 1) Load hurricane data

## 1.1) Load libraries

```{r}
library(curveDepth)
library(fda)
library(abind)
library(stringi)
library(ggplot2)
library(maps)
library(rnaturalearth)
library(rnaturalearthdata)
source("depth_routines.R")
```

## 1.2) Load and pre-process the data (https://www.nhc.noaa.gov/data/)

Pre-read the data and set limits and counters
```{r}
# https://www.nhc.noaa.gov/data/
fName <- "../data/hurdat2-1851-2019-052520.txt"
nLines <- length(readLines(fName))
nLinesRead <- 0
all.hurricanes <- list()
nHurricanes <- 0
strengths <- c(" HU", " TS", " EX")
all.lim.lattitude <- c(+89, -89)
all.lim.longitude <- c(+179, -179)
theTimeOrigin <- as.Date(as.Date.character("1850-01-01", format = "%Y-%m-%d"))
```

Read the hurricanes one by one
```{r}
while(nLinesRead < nLines){
  all.hurricanes.tmp.header <- read.table(fName, 
                                          skip = nLinesRead, nrows = 1, sep = ",")
  nLinesRead <- nLinesRead + 1
  # all.hurricanes.tmp.header[3] is the length
  all.hurricanes.tmp <- read.table(fName, 
                                   skip = nLinesRead, 
                                   nrows = 
                                     as.numeric(all.hurricanes.tmp.header[3]), 
                                   sep = ",")
  nLinesRead <- nLinesRead + as.numeric(all.hurricanes.tmp.header[3])
  # Read hurricane data:
  # Read latitude
  latitude <- as.numeric(strtrim(all.hurricanes.tmp[[5]], 5))
  latitude.dir <- stri_detect(all.hurricanes.tmp[[5]], fixed = "S")
  latitude[latitude.dir] <- -latitude[latitude.dir]
  # Update latitude limits
  all.lim.lattitude[1] <- min(all.lim.lattitude[1], latitude)
  all.lim.lattitude[2] <- max(all.lim.lattitude[2], latitude)
  # Read longitude
  longitude <- as.numeric(strtrim(all.hurricanes.tmp[[6]], 6))
  longitude.dir <- stri_detect(all.hurricanes.tmp[[6]], fixed = "W")
  longitude[longitude.dir] <- -longitude[longitude.dir]
  # Update longitude limits
  all.lim.longitude[1] <- min(all.lim.longitude[1], longitude)
  all.lim.longitude[2] <- max(all.lim.longitude[2], longitude)
  # Read strength
  strengths.dir <- all.hurricanes.tmp[[4]]
  if (strengths[1] %in% strengths.dir){
    strength.cur <- strengths[1]
  }else{
    if (strengths[2] %in% strengths.dir){
      strength.cur <- strengths[2]
    }else{
      strength.cur <- strengths[3]
    }
  }
  # Read the time axis
  absTime <- ISOdatetime(as.integer(all.hurricanes.tmp[[1]] / 10000), 
                         as.integer(all.hurricanes.tmp[[1]] %% 10000 / 100), 
                         as.integer(all.hurricanes.tmp[[1]] %% 100),
                         as.integer(all.hurricanes.tmp[[2]] / 100), 
                         as.integer(all.hurricanes.tmp[[2]] %% 100),
                         0)
  relTime <- as.numeric((absTime - absTime[1]) / (60 * 60))
  # Save the hurricane
  nHurricanes <- nHurricanes + 1
  all.hurricanes[[nHurricanes]] <- list(
    name = as.character(all.hurricanes.tmp.header$V1),
    timeStart = absTime[1],
    strength = strength.cur,
    args = relTime,
    vals = cbind(latitude, longitude), 
    coords = cbind(latitude, longitude))
}
```

Plot all hurricanes
```{r}
plot(NULL, xlim = all.lim.lattitude, ylim = all.lim.longitude)
for (i in 1:length(all.hurricanes)){
  lines(all.hurricanes[[i]]$vals)
}
```

# 2) Select hurricanes passing through the Gulf of Mexico

## 2.1) Prepare the structure and set the filter

```{r}
sel.hurricanes <- list()
year.min <- 1918
year.max <- 2021
coords.through <- c(24.51, -78.50)
coords.radius <- 65
desired.strength <- " HU"
min.points <- 0
```

## 2.2) Filter all the available hurricanes

```{r}
for (i in 1:length(all.hurricanes)){
  # Choose a circle for the trajectory
  if (as.numeric(format(all.hurricanes[[i]]$timeStart, "%Y")) >= year.min && 
      as.numeric(format(all.hurricanes[[i]]$timeStart, "%Y")) <= year.max){
    if (all.hurricanes[[i]]$strength == desired.strength){
      curHur.dists.through <- sqrt(rowSums((t((t(all.hurricanes[[i]]$vals) - coords.through) * c(60, 60)))^2))
      if (sum(curHur.dists.through < coords.radius) > 0.5) {
        cur.hurricane <- all.hurricanes[[i]]
        begind <- which(curHur.dists.through < coords.radius)[1]
        cur.hurricane$args <- cur.hurricane$args[begind:length(cur.hurricane$args)]
        cur.hurricane$vals <- cur.hurricane$vals[begind:nrow(cur.hurricane$vals)
                                                 ,,drop = FALSE]
        if (nrow(cur.hurricane$vals) >= min.points){
          sel.hurricanes[[length(sel.hurricanes) + 1]] <- cur.hurricane
        }
      }
    }
  }
}

```

## 2.3) Plot the selected hurricane tracks
```{r}
# Prepare a data frame with hurricanes
latitudeNs <- c()
longitudeEs <- c()
names <- c()
for (i in 1:length(sel.hurricanes)){
  latitudeNs <- c(latitudeNs, sel.hurricanes[[i]]$vals[,1])
  longitudeEs <- c(longitudeEs, sel.hurricanes[[i]]$vals[,2])
  names <- c(names, rep(sel.hurricanes[[i]]$name, 
                        nrow(sel.hurricanes[[i]]$vals)))
}
hurricanes.df <- data.frame(latitudeNs, longitudeEs, names)
colnames(hurricanes.df) <- c("latitudeN", "longitudeE", "name")
# Plot them on the map
world <- ne_countries(scale = 'medium', returnclass = 'sf')
ggulf <- ggplot(data = world) +
  geom_sf(aes(fill = NULL)) +
  annotate(geom = "text", x = -90, y = 26, label = "Gulf of Mexico", 
           fontface = "italic", color = "grey22", size = 6) +
  coord_sf(xlim = c(-150, -25), ylim = c(0, 75), expand = FALSE) +
  geom_path(aes(x = longitudeE, y = latitudeN, color = name), 
            data = hurricanes.df) + 
  scale_fill_viridis_d(option = "plasma") +
  theme(legend.position = "none", axis.title.x = element_blank(), 
        axis.title.y = element_blank(), panel.background = element_rect(fill = "azure"), 
        panel.border = element_rect(fill = NA))
plot(ggulf)

```

# 3) Calculate Tukey curve depth

## 3.1) Obtain parametrization
```{r}
paramsTime.data <- paramTime(sel.hurricanes, 100)
# Extract matrices for MFHD
y1 <- NULL
y2 <- NULL
# One observation is a column
for (i in 1:length(paramsTime.data)){
  rawData <- paramsTime.data[[i]]$vals
  y1 <- cbind(y1, rawData[,1])
  y2 <- cbind(y2, rawData[,2])
}
```

## 3.2) Calculate the Tukey curve depth
```{r}
library(curveDepth)
for (i in 1:length(paramsTime.data)){
  paramsTime.data[[i]]$coords = paramsTime.data[[i]]$vals
}
dssTCD <- depthc.Tukey(paramsTime.data, paramsTime.data, m = 250)
print(dssTCD)
```

# 4) Plot the two depths for comparison

## 4.1) Prepare a data frame with hurricanes
```{r}
latitudeNs <- c()
longitudeEs <- c()
color <- c()
clalp <- c()
lltype <- c()
track <- c()
for (i in 1:length(paramsTime.data)){
  # Write the latitude and longitude
  latitudeNs <- c(latitudeNs, paramsTime.data[[i]]$vals[,1])
  longitudeEs <- c(longitudeEs, paramsTime.data[[i]]$vals[,2])
  # Put placeholder for color/linetype
  color <- c(color, rep("a", 
                        nrow(paramsTime.data[[i]]$vals)))
  clalp <- c(clalp, rep(0, 
                        nrow(paramsTime.data[[i]]$vals)))
  lltype <- c(lltype, rep(0, 
                        nrow(paramsTime.data[[i]]$vals)))
  # Define grouping to sidtinguish single trajectories
  track <- c(track, rep(i, nrow(paramsTime.data[[i]]$vals)))
}
# Rescale the data to fit on the map
hurricanes.df <- data.frame(latitudeNs, longitudeEs, color, clalp, 
                            lltype, track)
colnames(hurricanes.df) <- c("latitudeN", "longitudeE", "color", "transp", 
                             "ltype", "track")
```

## 4.2) Plot the depth-colored curve-boxplot hurricanes on the map for TCD

```{r}
hurricanes.df.curve <- hurricanes.df
hurricanes.df.curve$color <- 
  as.vector(t(matrix(rep(getColorsBoxplot(dssTCD)$colors, 
                         nrow(paramsTime.data[[i]]$vals)), 
                     ncol = nrow(paramsTime.data[[i]]$vals), byrow = FALSE)))
hurricanes.df.curve$transp <- 
  as.vector(t(matrix(rep(getColorsBoxplot(dssTCD)$alphas, 
                         nrow(paramsTime.data[[i]]$vals)), 
                     ncol = nrow(paramsTime.data[[i]]$vals), byrow = FALSE)))
hurricanes.df.curve <- hurricanes.df.curve[order(hurricanes.df.curve$color),]
world <- ne_countries(scale = 'medium', returnclass = 'sf')
ggulf2 <- ggplot() +
  geom_sf(data = world, aes(fill = NULL, color = "border1")) +
  coord_sf(xlim = c(-110, -65), ylim = c(15, 50), expand = FALSE) +
  geom_path(aes(x = longitudeE, y = latitudeN, color = color, alpha = transp), 
            group = track, size = 1, linetype = 1, data = hurricanes.df.curve) +
  geom_path(aes(x = longitudeE, y = latitudeN, color = "order3"), 
            size = 1, data = hurricanes.df.curve[
              hurricanes.df.curve$color == "order3",]) + 
  scale_color_manual(values = c("order4" = "red", "order1" = "royalblue",
                                "order2" = "royalblue",
                                "order3" = "royalblue4",
                                "border1" = "black")) +
  theme(legend.position = "none", axis.title.x = element_blank(), 
        panel.grid = element_line(colour = "gray", linetype = 4), 
        axis.title.y = element_blank(), 
        panel.background = element_rect(fill = "azure"), 
        panel.border = element_rect(fill = NA))
plot(ggulf2)
```

# ?) Try further depths

## For comparison, implement further depths, e.g. from R-package ddalpha

Add a new chunk by clicking the *Insert Chunk* button on the toolbar or by pressing *Cmd+Option+I*.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the *Preview* button or press *Cmd+Shift+K* to preview the HTML file). 

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike *Knit*, *Preview* does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.
