---
title: "FjordLight Example"
author: "Bernard Gentili & Robert Schlegel"
date: "2023-10-19"
description: "This vignette walks the user through a viable workflow for using the functions within this package to work with the PAR dataset."
output: 
  rmarkdown::html_vignette:
    fig_caption: yes
vignette: >
  %\VignetteIndexEntry{FjordLight Example}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: bibliography.bib
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
  )
```

# Overview

The __`FjordLight`__ package has been built explicitly to facilitate the use of the Arctic fjord PAR dataset developed via the [FACE-IT project](https://www.face-it-project.eu/) as an extension of the work by @gattuso2020global. FACE-IT has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 869154.

Below one will find a range of examples that should help to clarify the useful order of these functions in a workflow. All of these examples are also available within the help files for the given functions.

```{r setup}
library(FjordLight)
library(raster)
```

# Caveat emptor

With 64 GB of RAM (a typical desktop) your computer should have no problem loading the data, even with the most extensive fjords. With 4 GB (a typical laptop) you can only handle the smaller fjords: 'kong', 'young', and 'nuup'. However, if you do not need the monthly `PARbottom` or `KPAR` time series, any computer should be able to load all fjords. In this case, when calling the `fl_LoadFjord()` function, set the `TS` argument to `FALSE` (this is the default value).

# Workflow

## List of available fjords
```{r eval=FALSE}
fl_ListFjords()
```

## Download a fjord

To save a file to your computer:

```{r eval=FALSE}
# Choose a fjord from the list above
fjord <- "kong"

# Note: this will require that a folder named 'data' exists in your current working directory
# One may see one's working directory written in small text next to the version of R in the console pane
# If the file has already been downloaded a message will be shown
fl_DownloadFjord(fjord, dirdata = "data")
```

For the purposes of this vignette we will use the small toy dataset that is included within __`FJordLight`__:

```{r}
# Chose to load all of the monthly bottom PAR values or not
WANT_TIME_SERIES <- TRUE

# Load the data
fjord <- "test"
fjorddata <- fl_LoadFjord(fjord, dirdata = system.file("extdata", package = "FjordLight"), TS = WANT_TIME_SERIES)
str(fjorddata, list.len = 15)
```


## Get geographic parameters

```{r}
flget_geoparameters(fjorddata)
```

## Get bathymetry

```{r fig.width = 6, fig.height=6}
# All depths (what = "o" ; o for Ocean), as raster
all_bathy <- flget_bathymetry(fjorddata, what = "o", mode = "raster", PLOT = TRUE)

# Coastal zone [0-200m] (what = "c" ; c for Coastal), as raster
coastal_bathy <- flget_bathymetry(fjorddata, what = "c", mode = "raster", PLOT = TRUE)

# Shallow zone [0-50m] (what = "sl" ; s for Shallow, l to add land), as raster
shallow_land <- flget_bathymetry(fjorddata, what = "sl", mode = "raster", PLOT = TRUE)

# Just land; note the difference in colour palette
just_land <- flget_bathymetry(fjorddata, what = "l", mode = "raster", PLOT = TRUE)

# As a data.frame (mode = "df" : longitude, latitude, depth)
sea <- flget_bathymetry(fjorddata, what = "s", mode = "df", PLOT = FALSE)
cz <- flget_bathymetry(fjorddata, what = "c", mode = "df", PLOT = FALSE)

# you may add letter "l" if you want land elevation
sealand <- flget_bathymetry(fjorddata, what = "sl", mode = "df")

# Example of structure
str(sealand)
```

## Get climatologies of PAR0m, Kpar, and PARbottom

```{r fig.width = 6, fig.height=6}
# PAR0m and PARbottom for July
P07 <- flget_climatology(fjorddata, optics = "PAR0m", period = "Clim", month = 7, PLOT = TRUE)
print(P07)
Pk7 <- flget_climatology(fjorddata, optics = "Kpar", period = "Clim", month = 7, PLOT = TRUE)
print(Pk7)
Pb7 <- flget_climatology(fjorddata, optics = "PARbottom", period = "Clim", month = 7, PLOT = TRUE)
print(Pb7)

# PARbottom Global
PbG <- flget_climatology(fjorddata, optics = "PARbottom", period = "Global", PLOT = TRUE)
print(PbG)

# PAR0m, Kpar, and PARbottom for year 2012 as 3 columns data.frame
P02012 <- flget_climatology(fjorddata, optics = "PAR0m", period = "Yearly", year = 2012, mode = "df")
k2012 <- flget_climatology(fjorddata, optics = "Kpar", period = "Yearly", year = 2012, mode = "df")
Pb2012 <- flget_climatology(fjorddata, optics = "PARbottom", period = "Yearly", year = 2012, mode = "df")
head(Pb2012)
```

If you want the data as a data.frame:

```{r}
# first get pixels area
area <- flget_area(fjorddata, mode = "df")

# Then bind the data frames and remove rows with missing values
PAR_area <- cbind(sea, area[3], P02012[3], k2012[3], Pb2012[3])
PAR_area <- PAR_area[complete.cases(PAR_area),]
head(PAR_area)
```

## Get monthly time series of PARbottom

```{r fig.width = 6, fig.height=6}
# Years 2003 to 2004 - months July to August
mts <- flget_PARbottomMonthlyTS(fjorddata, month = 7:8, year = 2003:2004, PLOT = TRUE)
print(mts)

# Or as a data.frame
mts_2003 <- flget_PARbottomMonthlyTS(fjorddata, year = 2003, PLOT = FALSE, mode = "df")
head(mts_2003)
```

```{r eval=FALSE}
# All months - all years - as data.frame: 
  # columns = months (8, March to October) * years (20, 2003 to 2022) + 2 (lon lat) = 162
# NB: This may be too large for some laptops, proceed with caution
mts_full <- flget_PARbottomMonthlyTS(fjorddata, mode = "df", PLOT = FALSE)
```

## Get monthly time series of KPAR

As an addendum to the original monthly bottom PAR data, data for K~PAR~ at the same spatial and temporal resolution were made available as a separate set of files. To download and work with these files one follows a similar workflow to that shown above:

```{r eval=FALSE}
# Choose a fjord from the possible options
fjord <- "kong"

# Note: this will require that a folder named 'data' exists in your current working directory
# One may see one's working directory written in small text next to the version of R in the console pane
# If the file has already been downloaded a message will be shown
fl_DownloadFjord(fjord, layer = "K_PAR" dirdata = "data")
```

For the purposes of this vignette we will use the `KPAR` test file included with the package.

```{r}
# Chose to load all of the monthly bottom PAR values or not
WANT_TIME_SERIES <- TRUE

# Load the data
fjord <- "test"
fjorddata_KPAR <- fl_LoadFjord(fjord, layer = "K_PAR", dirdata = system.file("extdata", package = "FjordLight"), TS = WANT_TIME_SERIES)
str(fjorddata_KPAR, list.len = 15)
```

Note that these files do not contain global attributes because they are a supplement to the original bottom PAR files.

Once the data have been loaded they may be worked with in the same way as the bottom PAR monthly time series data by changing the use of the function `flget_PARbottomMonthlyTS()` to `flget_KPARMonthlyTS()`.

```{r}
# Years 2003 to 2004 - months July to August
mts_KPAR <- flget_KPARMonthlyTS(fjorddata_KPAR, month = 7:8, year = 2003:2004, PLOT = TRUE)
print(mts)

# Or as a data.frame
mts_KPAR_2003 <- flget_KPARMonthlyTS(fjorddata_KPAR, year = 2003, PLOT = FALSE, mode = "df")
head(mts_2003)
```

## Standard deviations for climatologies

Also appended to the original bottom PAR dataset have been files that document the standard deviation (SD) for the monthly and annual climatologies. These may be downloaded/loaded by using the argument `layer = "ClimSD"` or `layer = "YearlySD"` in the corresponding functions as seen below:

```{r eval=FALSE}
# Choose a fjord from the possible options
fjord <- "kong"

# Note: this will require that a folder named 'data' exists in your current working directory
# One may see one's working directory written in small text next to the version of R in the console pane
# If the file has already been downloaded a message will be shown
fl_DownloadFjord(fjord, layer = "ClimSD", dirdata = "data")
fl_DownloadFjord(fjord, layer = "YearlySD", dirdata = "data")
```

As with the other examples above, for this vignette we will proceed with the pre-packaged example files. Note however that there are no time series data in these files, so we do not need to specificy this argument (default behaviour is `FALSE`).

```{r}
# Monthly climatology SD
fjorddata_ClimSD <- fl_LoadFjord(fjord, layer = "ClimSD", dirdata = system.file("extdata", package = "FjordLight"))
str(fjorddata_ClimSD, list.len = 15)
fjorddata_YearlySD <- fl_LoadFjord(fjord, layer = "YearlySD", dirdata = system.file("extdata", package = "FjordLight"))
str(fjorddata_YearlySD, list.len = 15)
```

To extract and work with the monthly climatology SD data one may then follow these steps:

```{r}
# Determine coordinates
lon <- fjorddata$longitude
lat <- fjorddata$latitude
str(lon); str(lat)

# Select a month
month <- 8 # August
im <- which(fjorddata$Months == month)
print(im)

# PAR0m Standard Deviation
PAR0mSD <- raster::raster(list(x = lon, y = lat, z = fjorddata_ClimSD$ClimPAR0mSD[, , im]))
plot(PAR0mSD, main = paste("PAR0m StDev", month.abb[month]))

# kdpar Standard Deviation
KparSD <- raster::raster(list(x = lon, y = lat, z = fjorddata_ClimSD$ClimKparSD[, , im]))
plot(KparSD, main = paste("Kpar StDev", month.abb[month]))

# PARbottom Standard Deviation
PARbottomSD <- raster::raster(list(x = lon, y = lat, z = fjorddata_ClimSD$ClimPARbottomSD[, , im]))
plot(PARbottomSD, main = paste("PARbottom StDev", month.abb[month]))

# PARbottom
PAR0m <- flget_climatology(fjorddata, optics = "PAR0m", period = "Clim", month = month, PLOT = TRUE)

VarCoef <- PAR0mSD / PAR0m
plot(VarCoef, main = "PAR0m Coefficient of Variation")
```

Likewise, for the yearly climatology SD data:

```{r}
# Determine coordinates
lon <- fjorddata$longitude
lat <- fjorddata$latitude
str(lon); str(lat)

# Choose a year
year <- 2007
iy <- which(fjorddata$Years == year)
print(iy)

# PAR0m Standard Deviation
PAR0mSD <- raster::raster(list(x = lon, y = lat, z = fjorddata_YearlySD$YearlyPAR0mSD[, , iy]))
plot(PAR0mSD, main = paste("PAR0m StDev", year))

# kdpar Standard Deviation
KparSD <- raster::raster(list(x = lon, y = lat, z = fjorddata_YearlySD$YearlyKparSD[, , iy]))
plot(KparSD, main = paste("Kpar StDev", year))

# PARbottom Standard Deviation
PARbottomSD <- raster::raster(list(x = lon, y = lat, z = fjorddata_YearlySD$YearlyPARbottomSD[, , iy]))
plot(PARbottomSD, main = paste("PARbottom StDev", year))

# PARbottom
PAR0m <- flget_climatology(fjorddata, optics = "PAR0m", period = "Yearly", year = year, PLOT = TRUE)

VarCoef <- PAR0mSD / PAR0m
plot(VarCoef, main = "PAR0m Coefficient of Variation")
```


## P-functions

```{r fig.width = 6, fig.height=6}
# One may create their own functions
fG <- flget_Pfunction(fjorddata, type = "coastal", period = "Global", plot = FALSE)
# Then you can use it; for instance :
irradiance_levels <- c(0.1, 1, 10)
fG(irradiance_levels)

# Or load the pre-calculated values as a 2 column data.frame
f2012 <- flget_Pfunction(fjorddata, type = "coastal", period = "Yearly", year = 2012, mode = "df")
str(f2012)

# Plot P-functions
fGlob <- flget_Pfunction(fjorddata, type = "coastal", period = "Global", PLOT = TRUE, lty = 1, col = 1, lwd = 2, 
                         Main = paste(fjord, "coastal P-functions"), ylim = c(0, 50))
```

# References