README-extended.Rmd

---
output: 
  github_document:
    pandoc_args: --webtex
    # --webtex=https://latex.codecogs.com/png.latex?%5Cdpi{150}
  pdf_document:
    extra_dependencies: ["float"]
header-includes:
   - \usepackage[default]{sourcesanspro}
   - \usepackage[T1]{fontenc}
bibliography: documentation/codos/inst/references/references.bib
cls: documentation/codos/inst/references/proceedings-of-the-royal-society-a.csl
---
<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%",
  dpi = 500,
  dev = "cairo_pdf"
)
```

# CO<sub>dos</sub>: CO<sub>2</sub> Correction Tools 
<!-- <img src="documentation/codos/inst/images/logo.png" alt="logo" align="right" height=200px/> -->

<!-- badges: start -->
`r if(!knitr::is_latex_output()) badger::badge_devel("special-uor/codos", "yellow")`
`r if(!knitr::is_latex_output()) badger::badge_cran_release("codos", "red")`
`r if(!knitr::is_latex_output()) badger::badge_github_actions("special-uor/codos")`
`r if(!knitr::is_latex_output()) badger::badge_doi("10.5281/zenodo.5083309", "blue")`
<!-- badges: end -->

## Installation

You can(not) install the released version of codos from [CRAN](https://CRAN.R-project.org) with:

``` r
install.packages("codos")
```

And the development version from [GitHub](https://github.com/) with:

``` r
# install.packages("devtools")
devtools::install_github("special-uor/codos", "dev")
```
<!-- ## Example -->

<!-- - CRU TS 4.04: [inst/extdocs/cru-ts-4.04.md](inst/extdocs/cru-ts-4.04.md) -->


```{r, include = !knitr::is_latex_output(), results = 'asis'}
cat(paste0("
-----

Note: Some of the equations on this document are not displayed properly (due to a server issue), check out the [README-extended.pdf](README-extended.pdf).

-----
"))
```

## Background:

### Vapour-pressure deficit (`vpd`)
`vpd` is given by mean daily growing season temperature, `tmp` [°C] and moisture index, `mi` [-]. Using the CRU TS 4.04 dataset [@cru404] we found the following relation:

<center>
$$\text{vpd} = 4.6123 \times \exp(0.0609 \times \text{tmp}-0.8726 \times \text{mi})$$
</center>

The steps performed were:

1. Generate a monthly climatology for the period between 1961 and 1990 (inclusive). Variables used: `cld`, `pre`, `tmn`, `tmx`, `vap`.
  
  ```{r, eval = FALSE, echo = TRUE}
  # Monthly climatology for `tmn`
  codos::monthly_clim("cru_ts4.04.1901.2019.tmn.dat.nc", "tmn", 1961, 1990)
  ```

Output file:

  ```bash
  "cru_ts4.04.1901.2019.tmn.dat-clim-1961-1990.nc"
  ```
  
2. Interpolate the monthly data to daily. Variables used: `cld`, `pre`, `tmn`, `tmx`, `vap`.

  ```{r, eval = FALSE, echo = TRUE}
  # Monthly to daily interpolation for `tmn`
  codos::nc_int("cru_ts4.04.1901.2019.tmn.dat-clim-1961-1990.nc", "tmn")
  ```
  
Output file:

  ```bash
  "cru_ts4.04.1901.2019.tmn.dat-clim-1961-1990-int.nc"
  ```
  
3. Calculate daily temperature, `tmp`. Variables used: `tmn` and `tmx`.

  ```{r, eval = FALSE, echo = TRUE}
  codos::daily_temp(tmin = list(filename = "cru_ts4.04.1901.2019.tmn.dat-clim-1961-1990-int.nc",
                                id = "tmn"),
                    tmax = list(filename = "cru_ts4.04.1901.2019.tmx.dat-clim-1961-1990-int.nc", 
                                id = "tmx"),
                    output_filename = "cru_ts4.04-clim-1961-1990-daily.tmp.nc")
  ```
  
4. Calculate mean growing season for daily temperature

  ```{r, eval = FALSE, echo = TRUE}
  codos::nc_gs("cru_ts4.04-clim-1961-1990-daily.tmp.nc", "tmp", thr = 0)
  ```

Output file:

  ```bash
  "cru_ts4.04-clim-1961-1990-daily.tmp-gs.nc"
  ```
  
5. Calculate potential evapotranspiration (`pet`)

Install `SPLASH` (unofficial R package) as follows:

  ```{r, eval = FALSE}
  remotes::install_github("villegar/splash", "dev")
  ```

Or, download from the official source: https://bitbucket.org/labprentice/splash.

  ```{r, eval = FALSE, echo = TRUE}
  elv <- codos:::nc_var_get("halfdeg.elv.nc", "elv")$data
  tmp <- codos:::nc_var_get("cru_ts4.04.1901.2019.daily.tmp.nc", "tmp")$data
  cld <- codos:::nc_var_get("cru_ts4.04.1901.2019.cld.dat-clim-1961-1990-int.nc", 
                            "cld")$data
  
  codos::splash_evap(output_filename = "cru_ts4.04-clim-1961-1990-pet.nc", 
                     elv, # Elevation, 720x360 grid 
                     sf = 1 - cld / 100, 
                     tmp, 
                     year = 1961, # Reference year 
                     lat = codos::lat, 
                     lon = codos::lon)
  ```

Output file:

  ```bash
  "cru_ts4.04-clim-1961-1990-pet.nc"
  ```
  
6. Calculate moisture index (`mi`)

<center>
$$MI_{i,j} = \frac{\text{Total precipitation}}{\text{Total PET}}$$
</center>

  ```{r, eval = FALSE, echo = TRUE}
  pet <- codos:::nc_var_get("cru_ts4.04-clim-1961-1990-pet.nc", 
                            "pet")$data
  pre <- codos:::nc_var_get("cru_ts4.04.1901.2019.pre.dat-new-clim-1961-1990-int.nc", 
                            "pre")$data
  codos::nc_mi(output_filename = "cru_ts4.04-clim-1961-1990-mi.nc", 
               pet, # potential evapotranspiration
               pre) # precipitation
  ```
  
Output file:

  ```bash
  "cru_ts4.04-clim-1961-1990-mi.nc"
  ```

7. Approximate `vpd`

  ```{r, eval = FALSE, echo = TRUE}
  tmp <- codos:::nc_var_get("cru_ts4.04-clim-1961-1990-daily.tmp.nc", 
                            "tmp")$data
  vap <- codos:::nc_var_get("cru_ts4.04.1901.2019.vap.dat-clim-1961-1990-int.nc", 
                            "vap")$data
  output_filename <- file.path(path, "cru_ts4.04-clim-1961-1990-vpd-tmp.nc")
  codos::nc_vpd(output_filename, tmp, vap)
  ```

Output file:

  ```bash
  "cru_ts4.04-clim-1961-1990-vpd-tmp.nc"
  ```
  
8. Find the coeffients for the following equation

<center>
$$\text{vpd} = a \times \exp(\text{kTmp} \times \text{tmp}-\text{kMI} \times \text{mi})$$
</center>

```{r, eval = FALSE, echo = TRUE}
mi <- codos:::nc_var_get("cru_ts4.04-clim-1961-1990-mi.nc", "mi")$data
Tmp <- codos:::nc_var_get("cru_ts4.04-clim-1961-1990-daily.tmp-gs.nc", "tmp")$data
vpd <- codos:::nc_var_get("cru_ts4.04-clim-1961-1990-vpd-tmp-gs.nc", "vpd")$data

# Apply ice mask
mi[codos:::ice_mask] <- NA
Tmp[codos:::ice_mask] <- NA
vpd[codos:::ice_mask] <- NA

# Filter low temperatures, Tmp < 5
mi[Tmp < 5] <- NA
Tmp[Tmp < 5] <- NA

# Create data frame
df <- tibble::tibble(Tmp = c(Tmp), 
                     vpd = c(vpd),
                     MI = c(mi))

# Filter grid cells with missing Tmp, vpd, or MI
df <- df[!is.na(df$Tmp) & !is.na(df$vpd) & !is.na(df$MI), ]

# Linear approximation
lmod <- lm(log(vpd) ~ Tmp + MI, data = df)
# Non-linear model
exp_mod <- nls(vpd ~ a * exp(kTmp * Tmp - kMI * MI),
               df,
               start = list(a = exp(coef(lmod)[1]),
                            kTmp = coef(lmod)[2],
                            kMI = coef(lmod)[3]),
               control = list(maxiter = 200))
```

Summary statistics:
```{r, echo = FALSE, eval = FALSE}
df <- read.csv("documentation/codos/inst/extdata/df.csv")
# Linear approximation
lmod <- lm(log(vpd) ~ Tmp + MI, data = df)
# Non-linear model
exp_mod <- nls(vpd ~ a * exp(kTmp * Tmp - kMI * MI),
               df,
               start = list(a = exp(coef(lmod)[1]),
                            kTmp = coef(lmod)[2],
                            kMI = coef(lmod)[3]),
               control = list(maxiter = 200))
```

```{r, eval = FALSE}
summary(exp_mod)
coefficients(exp_mod)
```

```{r, include=knitr::is_latex_output(), echo = FALSE}
knitr::asis_output('\n\n\\newpage')
```

### Corrected `mi` from reconstructed `mi`

The following equations were used:

<center>$$f = e/1.6 = D/[\text{c}_\text{a}(1-\chi)]$$</center>
<!-- OLD: -->
<!-- <center>$$\chi  =  [\xi/(\xi + \text{vpd}^ {1/2})] [1 - \Gamma^*/\text{c}_\text{a}] + \Gamma^*/\text{c}_\text{a}$$</center> -->
<!-- NEW: -->

<center>$$\chi = \frac{\xi \times \text{vpd}^{1/2} + \text{vpd}}{\text{c}_\text{a}/\left(\text{c}_\text{a}+9.7\right)}$$</center>

<center>$$\xi = [\beta (K + \Gamma^*) / (1.6 \eta^*)] ^ {1/2}$$</center>

where:

- $e$	ratio of water lost to carbon fixed [–]
- $\text{vpd}$	vapour pressure deficit [Pa]
- $\text{c}_\text{a}$	ambient CO2 partial pressure [Pa]
- $\chi$	ratio of leaf-internal to ambient CO2 partial pressures [–]
- $\xi$	stomatal sensitivity factor [Pa1/2]
- $\Gamma^*$	photorespiratory compensation point [Pa]: a function of temperature and elevation
- $\beta$	ratio of cost factors for carboxylation and transpiration = 146 [–]
- $K$	effective Michaelis constant of Rubisco [Pa]: a function of temperature and elevation
- $\eta^*$ viscosity of water relative to its value at 25°C [–]

And the equilibrium relation:

<center>$$f(\text{T}_\text{c1}, \text{MI}_\text{1}, \text{c}_\text{a,1}) =   f(\text{T}_\text{c0}, \text{MI}_\text{0}, \text{c}_\text{a,0})$$</center>
where:

- $\text{T}_\text{c1}$	past temperature (assume equal to reconstructed value) [K]
- $\text{MI}_\text{1}$	past MI (unknown) [–]
- $\text{c}_\text{a,1}$ past ambient CO2 partial pressure [Pa], adjusted for elevation
- $\text{T}_\text{c0}$	present temperature [K]
- $\text{MI}_\text{0}$	reconstructed MI [–]
- $\text{c}_\text{a,1}$	'recent' ambient CO2 partial pressure [Pa], adjusted for elevation

Steps in the solution:

1. Evaluate $f(\text{T}_\text{c0}, \text{MI}_\text{0}, \text{c}_\text{a,0})$
2. Equate this to:

<!-- OLD: -->
<!-- $$[\xi(\text{T}_\text{c1}, z) \text{vpd}_1^{1/2} + \text{vpd}_1] / [\text{c}_\text{a,1}(z) - \Gamma^*(\text{T}_\text{c1}, z)] $$ -->
<!-- NEW: -->
$$[\xi(\text{T}_\text{c1}, z) \times \text{vpd}_1 ^ {1/2} + \text{vpd}_1] / [\text{c}_\text{a,1}(z) / (\text{c}_\text{a,1}(z) + 9.7)] $$

where:

- $z$ is elevation 
- $\text{vpd}_1$ is past vapour pressure deficit

And solve for $\text{vpd}_1$.

3. Convert $\text{vpd}_1$ back to MI (at temperature $\text{T}_\text{c1}$), to yield an estimate of $\text{MI}_\text{1}$.

Using `codos`, all the steps translate to a simple function call

```{r, eval = FALSE}
corrected_mi <- codos::corrected_mi(present_t,
                                    past_temp,
                                    recon_mi,
                                    modern_co2,
                                    past_co2)
```

Note that this function takes temperatures in [°C] and ambient CO$_2$ partial pressures in 
[$\mu\text{mol}/\text{mol}$] (unless, `scale_factor` is overwritten, e.g. `scale_factor = 1` to use ambient CO$_2$ partial pressures in [Pa]).

More details:

```{r, eval = FALSE}
?codos::corrected_mi
```

# References

<!-- [1] University of East Anglia Climatic Research Unit; Harris, I.C.; Jones, P.D.;  -->
<!-- Osborn, T. (2020): CRU TS4.04: Climatic Research Unit (CRU) Time-Series (TS) -->
<!-- version 4.04 of high-resolution gridded data of month-by-month variation in -->
<!-- climate (Jan. 1901- Dec. 2019). Centre for Environmental Data Analysis. -->
<!-- <https://catalogue.ceda.ac.uk/uuid/89e1e34ec3554dc98594a5732622bce9> -->