vec-prep.Rmd

---
title: "Vector preparation"
author: "Robert Schlegel"
date: "2019-06-06"
output: workflowr::wflow_html
editor_options:
  chunk_output_type: console
csl: FMars.csl
bibliography: MHWNWA.bib
---

```{r global_options, include = FALSE}
knitr::opts_chunk$set(fig.width = 8, fig.align = 'center',
                      echo = TRUE, warning = FALSE, message = FALSE, 
                      eval = TRUE, tidy = FALSE)
```


## Introduction

It would seem logical that the U, V, and W vectors for water movement would be prepared at the same time as all of the other variables, as seen in the [Variable preparation](https://robwschlegel.github.io/MHWNWA/var-prep.html) vignette. The issue we face with the vectors is that they are stored in the 3D version of the model, which has a 5 day temporal resolution. This means that they will need to be handled slightly differently from the rest of the variables. Largely though we should be able to adjust the functions previously created for the other variables.

```{r startup}
# Packages used in this vignette
library(jsonlite, lib.loc = "../R-packages/")
library(tidyverse) # Base suite of functions
library(ncdf4) # For opening and working with NetCDF files
library(lubridate) # For convenient date manipulation

# Load functions required below
source("code/functions.R")

# Set number of cores
doMC::registerDoMC(cores = 50)

# Disable scientific notation for numeric values
  # I just find it annoying
options(scipen = 999)

# Corners of the study area
NWA_corners <- readRDS("data/NWA_corners.Rda")

# The NAPA surface data location
NAPA_files <- dir("../../data/NAPA025/1d_grid_T_2D", full.names = T)

# The NAPA 5 day U vector data location
NAPA_U_files <- dir("../../data/NAPA025/5d_grid_U", full.names = T)

# The NAPA 5 day V vector data location
NAPA_V_files <- dir("../../data/NAPA025/5d_grid_V", full.names = T)

# The NAPA 5 day W vector data location
NAPA_W_files <- dir("../../data/NAPA025/5d_grid_W", full.names = T)

# The NAPA model lon/lat values
NAPA_coords <- readRDS("data/NAPA_coords.Rda")

# The NAPA model lon/lat values subsetted to the study area
NAPA_coords <- readRDS("data/NAPA_coords_sub.Rda")

# Load NAPA bathymetry/lon/lat
NAPA_bathy <- readRDS("data/NAPA_bathy.Rda")

# Load NAPA bathymetry/lon/lat subsetted to study area
NAPA_bathy <- readRDS("data/NAPA_bathy_sub.Rda")

# Load MHW results
NAPA_MHW_sub <- readRDS("data/NAPA_MHW_sub.Rda")

# MHW Events
NAPA_MHW_event <- NAPA_MHW_sub %>%
  select(-clims, -cats) %>%
  unnest(events) %>%
  filter(row_number() %% 2 == 0) %>%
  unnest(events)
```


## Vector variables

We obviously want the primary vectors from the model, but I noticed that there are other things stored in these 5 day NetCDf files, so let's have a closer look.

```{r ncdf-var-dump}
# The U variables
NAPA_U <- ncdump::NetCDF(NAPA_U_files[1])$variable[1:6]
NAPA_U

# The V variables
NAPA_V <- ncdump::NetCDF(NAPA_V_files[1])$variable[1:6]
NAPA_V

# The W variables
NAPA_W <- ncdump::NetCDF(NAPA_W_files[1])$variable[1:6]
NAPA_W
```

It appears that the U and V files also contain their corresponding downward stress (`utau` and `vtau`), though with the availability of the W vectors, I don't know that we need these, too. For now we will just stick with: `uoce`, `voce`, and `wo`.


## Synoptic states

As with the other variables, we will also need to create mean synoptic states of the vectors (U, V, W) during all of the recorded MHWs. From here out we will again be running on data stored on a local server, meaning the following code won't be of much use unless one has a copy of the NAPA model data. The resulting synoptic states we will create now are likely to be in excess of a couple of GB so will not be hosted on the public GitHub server either. They are of course available upon request.


### Date index for NAPA vector files

To create the index of dates to be found within each of the thousands of NAPA 5 day vector NetCDF files we will use the same simple for loop to crawl through the files and write down for us in one long spreadsheet which dates are to be found in which files. This worked well for the previous vignette, so I've kept this step the same. The vectors are stored in three different file directories, but their dates are the same, so we only need to create one date index file, not three individual ones.

```{r date-index, eval=FALSE}
# Pull out the dates
NAPA_vector_files_dates <- data.frame()
for(i in 1:length(NAPA_files)){
  file_name <- NAPA_U_files[i]
  date_start <- ymd(str_sub(basename(as.character(file_name)), start = 26, end = 33))
  date_end <- ymd(str_sub(basename(as.character(file_name)), start = 35, end = 43))
  date_seq <- seq(date_start, date_end, by = "day")
  # Take the middle date as this is the real data in the pentad
  date_info <- data.frame(file = file_name, date = date_seq[3])
  NAPA_vector_files_dates <- rbind(date_info, NAPA_vector_files_dates)
}

# Order by date, just for tidiness
NAPA_vector_files_dates <- dplyr::arrange(NAPA_vector_files_dates, date)

# Save
# saveRDS(NAPA_vector_files_dates, "data/NAPA_vector_files_dates.Rda")
```


### Variable climatologies

Part of the data packet we need to create for the SOMs is the anomaly values. In order to create anomalies however we need to first create climatologies for all of the variables. This may prove to be a somewhat daunting task, but it's what we are here to do! In order to create a climatology of values we will need to load all of the files and then pixel-wise go about getting the seasonal (daily) climatologies. This will be done with the same function (`ts2clm()`) that is used for the MHW climatologies. We will first create a function that extracts the desired variables from any NetCDF files fed to it. With that done it should be a routine matter to get the climatologies. Hold onto your hats, this is going to be RAM heavy...

```{r clim-var-all, eval=FALSE}
# Load functions required below
source("code/functions.R")

# U
clim_one_vec(NAPA_U_files)
gc()

# V
clim_one_vec(NAPA_V_files)
gc()

# W
clim_one_vec(NAPA_W_files)
gc()
```

And then combine them into one big clim file.

```{r clim-var-all-join, eval=FALSE}
# Load all variable climatologies
NAPA_clim_U <- readRDS("data/NAPA_clim_U.Rda")
NAPA_clim_V <- readRDS("data/NAPA_clim_V.Rda")
NAPA_clim_W <- readRDS("data/NAPA_clim_W.Rda")

# left join them
NAPA_clim_vecs <- left_join(NAPA_clim_U, NAPA_clim_V, by = c("lon", "lat", "doy")) %>% 
  left_join(NAPA_clim_W, by = c("lon", "lat", "doy"))

# Convert DOY to MM-DD for joining to daily data below
NAPA_clim_vecs$doy <- format(as.Date(NAPA_clim_vecs$doy, origin = "2015-12-31"), "%m-%d")

# Change column names to highlight that these are climatology values
colnames(NAPA_clim_vecs)[-c(1:3)] <- paste0(colnames(NAPA_clim_vecs)[-c(1:3)],"_clim")

# Save
saveRDS(NAPA_clim_vecs, "data/NAPA_clim_vecs.Rda")
rm(NAPA_clim_U, NAPA_clim_V, NAPA_clim_W); gc()
```


### Variable extractor

We needed a list of the dates present in each file so that we can easily load only the NetCDF files we need to extract our desired variables. The dates we want are the range of dates during each of the MHWs detected in the [SST preparation](https://robwschlegel.github.io/MHWNWA/sst-prep.html) vignette. In the chunk below we will create a function that decides which files should have their variables loaded and a function that binds everything up into tidy data packets that our SOM can ingest.

```{r extractor-funcs}
# Load NAPA file date index
NAPA_vector_files_dates <- readRDS("data/NAPA_vector_files_dates.Rda")

# Load full variable climatology file
NAPA_clim_vecs <- readRDS("data/NAPA_clim_vecs.Rda")

# Function for extracting the desired variables from a given NetCDF file
# testers...
# file_number <- 1376
extract_all_vec <- function(file_number){
  
  # Extract and join variables with a for loop
    # NB: This should be optimised...
  u_vecs <- extract_one_vec(NAPA_U_files[file_number])
  v_vecs <- extract_one_vec(NAPA_V_files[file_number])
  if(file_number == 276){
    w_vecs <- v_vecs[1,]
    colnames(w_vecs)[6] <- "wo"
  } else {
    w_vecs <- extract_one_vec(NAPA_W_files[file_number])
  }
  
  NAPA_vecs_extracted <- left_join(u_vecs, v_vecs, by = colnames(u_vecs)[1:5]) %>% 
    left_join(w_vecs, by = colnames(u_vecs)[1:5])
  
  # Exit
  return(NAPA_vecs_extracted)
}

# Function for extracting variables from as many files as a MHW event requires
# testers...
# event_sub <- NAPA_MHW_event[23,]
data_vec_packet <- function(event_sub){
  
  # Create date and file index for loading
  date_idx <- seq(event_sub$date_start, event_sub$date_end, by = "day")
  file_idx <- which(NAPA_vector_files_dates$date %in% date_idx)
  
  # Load required base data
  # system.time(
  packet_base <- plyr::ldply(file_idx, extract_all_vec) %>% 
    filter(t %in% date_idx) %>%
    mutate(doy = format(t, "%m-%d"))
  # ) # 16 seconds for seven files
  
  # Join to climatologies
  packet_join <- left_join(packet_base, NAPA_clim_vecs, by = c("lon", "lat", "doy"))
  
  # Create anomaly values and remove clim columns
  packet_anom <- packet_join %>% 
    mutate(uoce_anom = uoce - uoce_clim,
           voce_anom = voce - voce_clim,
           wo_anom = wo - wo_clim) %>% 
    dplyr::select(lon, lat, doy, uoce:wo_anom, 
                  -c(colnames(NAPA_clim_vecs)[-c(1:3)]))
  
  # Create mean synoptic values
  packet_mean <- packet_anom %>% 
    select(-doy) %>% 
    # NB: The lowest pixels are a forcing edge and shouldn't be included
      # We can catch these out by filtering pixels whose SST is exactly 0
    # filter(sst != 0) %>% 
    group_by(lon, lat) %>% 
    summarise_all(mean, na.rm = T) %>% 
    arrange(lon, lat) %>% 
    ungroup() %>% 
    nest(.key = "synoptic")
  
  # Combine results with MHW dataframe
  packet_res <- cbind(event_sub, packet_mean)
  
  # Test visuals
  # ggplot(packet_mean, aes(x = lon, y = lat)) +
  #   geom_point(aes(colour = uoce_anom)) +
  #   scale_colour_gradient2(low = "blue", high = "red") +
  #   coord_cartesian(xlim = NWA_corners[1:2],
  #                   ylim = NWA_corners[3:4])
  
  # Exit
  return(packet_res)
}
```

With our functions sorted, it is now time to create our data packets.

```{r synoptic-states, eval=FALSE}
# Set number of cores
  # NB: Was set to 25 as someone else was using the server at the time
doMC::registerDoMC(cores = 25)

# Create one big packet
system.time(
synoptic_vec_states <- plyr::ddply(NAPA_MHW_event, c("region", "sub_region", "event_no"), data_vec_packet, .parallel = T)
) # 28 seconds for first 2, 6125 seconds (102 minutes) for all

# Save
# saveRDS(synoptic_vec_states, "data/synoptic_vec_states.Rda")
```

With all of our synoptic snapshots for our chosen variables created it is now time to feed them to the [Self-organising map (SOM) analysis](https://robwschlegel.github.io/MHWNWA/som.html).