data-prep.Rmd

---
title: "Preparing the data"
author: "Robert Schlegel"
date: "2020-02-25"
output: workflowr::wflow_html
editor_options:
  chunk_output_type: console
csl: FMars.csl
bibliography: MHWflux.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

The purpose of this vignette is to walk through the steps taken to prepare the data from the NOAA OISST, ERA5, and GLORYS products for the following analyses. The statistical analyses may be found in the [MHWs vs. heat flux]() vignette, and the SOM analysis may be found in the [SOM]() vignette. 

All of the libraries and functions used in this vignette, and the project more broadly, may be found [here](https://github.com/robwschlegel/MHWflux/blob/master/code/functions.R).

```{r satrtup}
# get everything up and running in one go
source("code/functions.R")
library(ggpubr)
library(gridExtra)
```

## Study area

A reminder of what the study area looks like. It has been cut into 6 regions, adapted from work by @Richaud2016.

```{r region-fig}
frame_base +
  geom_polygon(data = NWA_coords, alpha = 0.7, size = 2,
               aes(fill = region, colour = region)) +
  geom_polygon(data = map_base, aes(group = group))
```

## Pixels per region

In this study it was decided to use the higher resolution 1/12th degree GLORYS data. This means we will need to re-calculate which pixels fall within which region so we can later determine how to create our average SST time series per region as well as the other averaged heat flux term time series.

```{r grid-points, eval=FALSE}
# Load one OISST file to extract the lon/lat coords
OISST_grid <- readRDS("../data/OISST/daily/1982/daily.1982-01-01.Rda") %>% 
  dplyr::select(lon, lat) %>% 
  unique()

# Load one high-res GLORYS file to extract the lon/lat coords
GLORYS_grid <- tidync(GLORYS_files[1]) %>% 
  hyper_tibble() %>% 
  dplyr::rename(lon = longitude, lat = latitude) %>% 
  dplyr::select(lon, lat) %>% 
  unique()

# Load one ERA5 file to get the lon/lat coords
ERA5_grid <- tidync(ERA5_files[1]) %>%
  hyper_filter(latitude = dplyr::between(latitude, min(NWA_coords$lat), max(NWA_coords$lat)),
               longitude = dplyr::between(longitude, min(NWA_coords$lon)+360, max(NWA_coords$lon)+360),
               time = index == 1) %>%
  hyper_tibble() %>%
  dplyr::rename(lon = longitude, lat = latitude) %>%
  dplyr::select(lon, lat) %>%
  unique() %>%
  mutate(lon = lon-360) %>%  # Change back to +-180 scale
  mutate(lon = lon+0.125, lat = lat-0.125) # Regrid to match OISST coords
  
# Function for finding and cleaning up points within a given region polygon
points_in_region <- function(region_in, product_grid){
  region_sub <- NWA_coords %>% 
    filter(region == region_in)
  coords_in <- product_grid %>% 
    mutate(in_grid = sp::point.in.polygon(point.x = product_grid[["lon"]], point.y = product_grid[["lat"]], 
                                          pol.x = region_sub[["lon"]], pol.y = region_sub[["lat"]])) %>% 
    filter(in_grid >= 1) %>% 
    dplyr::select(-in_grid) %>% 
    mutate(region = region_in)
  return(coords_in)
}

# Run the function
OISST_regions <- plyr::ldply(unique(NWA_coords$region), points_in_region, 
                              .parallel = T, product_grid = OISST_grid)
saveRDS(OISST_regions, "data/OISST_regions.Rda")
GLORYS_regions <- plyr::ldply(unique(NWA_coords$region), points_in_region, 
                              .parallel = T, product_grid = GLORYS_grid)
saveRDS(GLORYS_regions, "data/GLORYS_regions.Rda")
ERA5_regions <- plyr::ldply(unique(NWA_coords$region), points_in_region, 
                            .parallel = T, product_grid = ERA5_grid)
saveRDS(ERA5_regions, "data/ERA5_regions.Rda")
```

```{r grid-points-visual}
OISST_regions <- readRDS("data/OISST_regions.Rda")
GLORYS_regions <- readRDS("data/GLORYS_regions.Rda")
# ERA5_regions <- readRDS("data/ERA5_regions.Rda")

# Combine for visual
ALL_regions <- rbind(OISST_regions, GLORYS_regions) %>% #, ERA5_regions) %>% 
  mutate(product = c(rep("OISST", nrow(OISST_regions)),
                     rep("GLORYS", nrow(GLORYS_regions))))#,
                     # rep("ERA5", nrow(ERA5_regions))))

# Visualise to ensure success
ggplot(NWA_coords, aes(x = lon, y = lat)) +
  geom_point(data = ALL_regions, aes(colour = region)) +
  geom_polygon(data = map_base, aes(group = group), show.legend = F) +
    coord_cartesian(xlim = NWA_corners[1:2],
                  ylim = NWA_corners[3:4]) +
  labs(x = NULL, y = NULL) +
  facet_wrap(~product)
```

## Average time series per region

With our pixels per region sorted we may now go about creating the average time series for each region from the OISST, GLORYS, and ERA5 data. First we will load a brick of the data constrained roughly to the study area into memory before assigning the correct pixels to their regions. Once the pixels are assigned we will summarise them into one mean time series per variable per region. These mean time series are what the rest of the analyses will depend on.

The code for prepping the OISST data.

```{r OISST_prep, eval=FALSE}
# The files with data in the study area
OISST_files_sub <- data.frame(files = OISST_files,
                              lon = c(seq(0.125, 179.875, by = 0.25), 
                                      seq(-179.875, -0.125, by = 0.25))) %>% 
  filter(lon >= min(NWA_coords$lon), lon <= max(NWA_coords$lon)) %>% 
  mutate(files = as.character(files))

# Load the full grid for the study area
system.time(
OISST_all <- plyr::ldply(OISST_files_sub$files, .fun = load_OISST, .parallel = TRUE)
) # 16 seconds

# Create the region time series and save
system.time(
OISST_all_ts <- OISST_all %>%
  right_join(OISST_regions, by = c("lon", "lat")) %>% 
  group_by(region, t) %>%
  summarise(temp = mean(temp, na.rm = T), .groups = "drop")
) # 14 seconds
saveRDS(OISST_all_ts, "data/OISST_all_ts.Rda")
```

The code for loading and processing the GLORYS data.

```{r GLORYS-prep, eval=FALSE}
# Set number of cores
  # NB: This is very RAM heavy, be careful with core use
doParallel::registerDoParallel(cores = 26)

# The GLORYS file location
system.time(
GLORYS_all_ts <- plyr::ldply(GLORYS_files, load_GLORYS_region, .parallel = T) %>% 
  dplyr::arrange(region, t) %>% 
  mutate(cur_spd = round(sqrt(u^2 + v^2), 2),
         cur_dir = round((270-(atan2(v, u)*(180/pi)))%%360))
) # 103 seconds on 26 cores
saveRDS(GLORYS_all_ts, "data/GLORYS_all_ts.Rda")
```

The code for the ERA5 data. Note that the ERA5 data are on an hourly 0.25 degree spatiotemporal grid. This loading process constrains them to a daily 0.25 degree grid that matches the OISST data before finally converting them to a single time series per region.

```{r ERA5-prep, eval=FALSE}
# See the code/workflow script for the code used for ERA5 data prep
# There is too much code to run from an RMarkdown document as each variable must loaded individually
```

## MHWs per region

We will be using the SST values from GLORYS for calculating the MHWs and will use the standard Hobday definition with a base period of 1993-01-01 to 2018-12-25. We are using an uneven length year as the data do not quite extend to the end of December. It was decided that the increased accuracy of the climatology from the 2018 year outweighed the negative consideration of having a clim period that excludes a few days of winter.

```{r MHW-calc, eval=FALSE}
# Load the data
OISST_all_ts <- readRDS("data/OISST_all_ts.Rda")

# Calculate the MHWs
OISST_region_MHW <- GLORYS_all_ts %>%
  dplyr::select(region:temp) %>% 
  group_by(region) %>%
  nest() %>%
  mutate(clims = map(data, ts2clm,
                     climatologyPeriod = c("1993-01-01", "2018-12-25")),
         events = map(clims, detect_event),
         cats = map(events, category, S = FALSE)) %>%
  select(-data, -clims)

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

### MHW results

As this is not really the focus of the paper, I will quickly go over the MHW results here. Specifically we are looking at the summary stats for MHWs per region and season, and how they may differ between those classifying groups. The stats created here are used in the results section of the manuscript and some of the musings will be used for the discussion.

```{r MHW-results}
MHW_annual_count <- OISST_MHW_event %>%
  mutate(year = year(date_peak)) %>% 
  group_by(year) %>% 
  summarise(count = n())
```

## Clims + anoms per variable

The analyses to come are going to be performed on anomaly values, not the original time series. In order to calculate the anomalies we are first going to need the climatologies for each variable. We will use the Hobday definition of climatology creation and then subtract the expected climatology from the observed values. We are again using the 1993-01-01 to 2018-12-25 base period for these calculations to ensure consistency throughout the project.

```{r clims, eval=FALSE}
# Load the data
GLORYS_all_ts <- readRDS("data/GLORYS_all_ts.Rda")
ERA5_all_ts <- readRDS("data/ERA5_all_ts.Rda")
OISST_all_ts <- readRDS("../MHWNWA/data/OISST_region.Rda")
ALL_ts <- left_join(ERA5_all_ts, GLORYS_all_ts, by = c("region", "t")) %>% 
  left_join(OISST_all_ts, by = c("region", "t")) %>% 
  filter(t <= "2018-12-31")

# Calculate GLORYS clims and anoms
  # Also give better names to the variables
  # Also convert the Qx terms from seconds to days by multiplying by 86,400
ALL_anom <- ALL_ts %>%
  dplyr::rename(lwr = msnlwrf, swr = msnswrf, lhf = mslhf, 
                shf = msshf, mslp = msl, sst = temp.y) %>% 
  dplyr::select(-wind_dir, -cur_dir, -temp.x) %>% 
  mutate(qnet_mld = (qnet*86400)/(mld*1024*4000),
         lwr_mld = (lwr*86400)/(mld*1024*4000),
         swr_mld = (swr*86400)/(mld*1024*4000),
         lhf_mld = (lhf*86400)/(mld*1024*4000),
         shf_mld = (shf*86400)/(mld*1024*4000),
         mld_1 = 1/mld) %>% 
  pivot_longer(cols = c(-region, -t), names_to = "var", values_to = "val") %>% 
  group_by(region, var) %>%
  nest() %>%
  mutate(clims = map(data, ts2clm, y = val, roundClm = 10,
                     climatologyPeriod = c("1993-01-01", "2018-12-25"))) %>% 
  dplyr::select(-data) %>% 
  unnest(cols = clims) %>%
  mutate(anom = val-seas) %>% 
  ungroup()

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

## Cumulative heat flux terms

We also need to create cumulative heatflux terms as well as a few other choice variables. This is done by taking the first day during the MHW and adding the daily values together cumulatively until the end of the event. The daily values are first divided by the MLD on that day as seen above. The MLD value used to divide the daily variables accounts for the water density and specific heat constant: Q/(rho x Cp x hmld), where rho = 1042 and Cp ~= 4000. Th Qnet term calculated this way approximates the air-sea flux term.

The movement terms aren't very useful and may not be worth including as they don't really show advection. So rather one can say that the parts of the heating that aren't explained by anything else could be attributed to advection through the process of elimination. For the moment they are still left in here.

```{r cum-heat-flux, eval=FALSE}
# We're going to switch over to the NOAA OISST data for MHWs now
ALL_anom_cum <- ALL_anom %>%
  dplyr::select(region, var, t, anom) %>% 
  pivot_wider(id_cols = c("region", "t"), names_from = var, values_from = anom) %>% 
  dplyr::select(region:tcc, mslp, qnet, p_e, mld, mld_1, qnet_mld:shf_mld) %>% 
  # Change this line depending on GLORYS or NOAA use
  left_join(OISST_MHW_clim[,c("region", "t", "event_no")], by = c("region", "t")) %>% 
  filter(event_no > 0) %>% 
  group_by(region, event_no) %>% 
  mutate_if(is.numeric, cumsum) %>% 
  ungroup() %>% 
  dplyr::select(region, event_no, t, everything()) %>% 
  pivot_longer(cols = c(-region, -event_no, -t), names_to = "var", values_to = "anom") %>% 
  mutate(var = paste0(var,"_cum")) %>% 
  dplyr::select(region, var, event_no, t, anom)

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

In the [MHWs vs. heat flux](https://robwschlegel.github.io/MHWflux/mhw-flux.html) vignette we will take the periods of time over which MHWs occurred per region and pair those up with the GLORYS and ERA5 data. This will be used to investigate which drivers are best related to the onset and decline of MHWs. In the [SOM](https://robwschlegel.github.io/MHWNWA/som.html) vignette we will look for how 

## References