---
title: "Exposure through time"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Exposure through time}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

library(climniche)

case_path <- system.file("extdata/mediterranean_anchovy", package = "climniche")
range_summary <- read.csv(
  file.path(case_path, "anchovy_climniche_time_range_summary.csv")
)
departure_summary <- read.csv(
  file.path(case_path, "anchovy_climniche_time_departure_summary.csv")
)
change_rate <- read.csv(
  file.path(case_path, "anchovy_climniche_time_change_rate.csv")
)
```

A single future projection cannot distinguish whether Niche Boundary
Exceedance spreads to additional cells or grows within cells already beyond
the boundary. This example follows European anchovy (*Engraulis encrasicolus*)
in the Mediterranean Sea at 2030, 2050, 2070 and 2090 under Bio-ORACLE
SSP2-4.5.

At each projection, Niche Boundary Exceedance is a future state relative to
the fixed empirical radial boundary. It is not a before-and-after difference
or an observed date of boundary crossing.

## Fit one climatic reference

The analysis uses the six Bio-ORACLE variables and the continuous anchovy SDM
suitability layer prepared in the
[Mediterranean example](climniche-examples.html). The current realised
climatic niche is fitted once. Its centre, climatic weighting matrix and
empirical radial boundary are then held fixed across all four projections.

```{r series-fit, eval = FALSE}
projection_years <- c(2030, 2050, 2070, 2090)

future_series <- list(
  ssp245_2030 = climate_2030,
  ssp245_2050 = climate_2050,
  ssp245_2070 = climate_2070,
  ssp245_2090 = climate_2090
)

# Project each period into the same fitted climatic niche.
series <- fit_climniche_series(
  current = current_climate,
  future = future_series,
  time = projection_years,
  scenario = "SSP2-4.5",
  occupied = anchovy_suitability,
  occupied_threshold = sdm_threshold,
  domain = mediterranean_sea_mask,
  sensitivity = climatic_weights
)
```

The order of the future list must match `time`. Current and future rasters use
the same variables, names, resolution and cell alignment.

## Range-level exposure

The temporal summaries are derived from Niche Boundary Exceedance; they are not
additional cell-level exposure quantities. Let $E_{it}$ denote Niche Boundary
Exceedance for cell $i$ at time $t$, $B$ the fitted radial boundary
distance and $w_i$ the product of current SDM suitability and cell area.
With `boundary_exceedance_tolerance = 0`,

$$
F_t =
\frac{\sum_i w_i I(E_{it} > 0)}
     {\sum_i w_i}
$$

is the Weighted Niche Boundary Exceedance Fraction,

$$
S_t =
\frac{\sum_i w_i (E_{it}/B) I(E_{it} > 0)}
     {\sum_i w_i I(E_{it} > 0)}
$$

is Conditional Relative Niche Boundary Exceedance, and

$$
M_t =
\frac{\sum_i w_i (E_{it}/B) I(E_{it} > 0)}
     {\sum_i w_i}
= F_t S_t
$$

is Range Mean Relative Niche Boundary Exceedance. Here, $I(\cdot)$ is the
indicator function.

```{r range-code, eval = FALSE}
range_summary <- climniche_range_summary(
  series,
  scope = "current",
  area_weight = TRUE
)

subset(
  range_summary,
  select = c(
    time,
    exposed_fraction,
    conditional_relative_exceedance,
    range_wide_relative_exceedance
  )
)
```

```{r range-table, echo = FALSE}
range_table <- data.frame(
  `Projection year` = range_summary[["time"]],
  `Weighted boundary exceedance fraction` = paste0(
    round(100 * range_summary[["exposed_fraction"]], 1),
    "%"
  ),
  `Conditional relative Niche Boundary Exceedance` = paste0(
    round(100 * range_summary[["conditional_relative_exceedance"]], 1),
    "%"
  ),
  `Range mean relative Niche Boundary Exceedance` = paste0(
    round(100 * range_summary[["range_wide_relative_exceedance"]], 1),
    "%"
  ),
  check.names = FALSE
)

knitr::kable(range_table, row.names = FALSE)
```

```{r range-figure-code, eval = FALSE}
extent_plot <- plot_climniche_time(
  series,
  metric = "exposed_fraction",
  scope = "current",
  area_weight = TRUE,
  show_models = FALSE
)

conditional_plot <- plot_climniche_time(
  series,
  metric = "conditional_relative_exceedance",
  scope = "current",
  area_weight = TRUE,
  show_models = FALSE
)

range_plot <- plot_climniche_time(
  series,
  metric = "range_wide_relative_exceedance",
  scope = "current",
  area_weight = TRUE,
  show_models = FALSE
)

patchwork::wrap_plots(
  extent_plot,
  conditional_plot,
  range_plot,
  nrow = 1
)
```

```{r range-figure-output, echo = FALSE, out.width = "100%"}
knitr::include_graphics("figures/anchovy-climniche-through-time.png")
```

The suitability- and area-weighted fraction beyond the boundary increases from
7.7% in 2030 to 78.9% in 2090. Conditional Relative Niche Boundary Exceedance
falls from 27.6% to 12.5% over the same period. This combination is consistent
with the addition of cells with relatively small exceedance. Range Mean
Relative Niche Boundary Exceedance rises from 2.1% to 9.9%.

## Timing within the current distribution

The first exceedance year identifies the first supplied projection in which a
cell lies beyond the current niche boundary. The projection fraction records
how many of the four supplied future periods exceed that boundary. Both
summaries refer to the sampled projections and do not estimate an intervening
crossing date.

```{r departure-code, eval = FALSE}
departure <- climniche_departure(
  series,
  scope = "current"
)

rate <- climniche_change_rate(
  series,
  metric = "range_wide_relative_exceedance",
  scope = "current",
  area_weight = TRUE
)
```

```{r departure-summary, echo = FALSE}
departure_table <- data.frame(
  `Boundary exceedance fraction` = paste0(
    round(
      100 * departure_summary[["proportion_with_boundary_exceedance"]],
      1
    ),
    "%"
  ),
  `Median first exceedance` =
    departure_summary[["median_first_boundary_exceedance"]],
  `Mean projection fraction` = paste0(
    round(
      100 * departure_summary[[
        "mean_boundary_exceedance_projection_fraction"
      ]],
      1
    ),
    "%"
  ),
  check.names = FALSE
)

knitr::kable(departure_table, row.names = FALSE)
```

```{r departure-map-code, eval = FALSE}
onset_map <- plot_climniche_departure_map(
  series,
  metric = "first_boundary_exceedance",
  scope = "current",
  scenario = "SSP2-4.5",
  study_region = mediterranean_boundary,
  degree_labels = "hemisphere"
)

frequency_map <- plot_climniche_departure_map(
  series,
  metric = "boundary_exceedance_projection_fraction",
  scope = "current",
  scenario = "SSP2-4.5",
  study_region = mediterranean_boundary,
  degree_labels = "hemisphere"
)

patchwork::wrap_plots(onset_map, frequency_map, nrow = 1)
```

```{r departure-map-output, echo = FALSE, out.width = "100%"}
knitr::include_graphics("figures/anchovy-climniche-time-maps.png")
```

The map distinguishes early exceedance from cells whose future climate first
appears beyond the boundary in 2090. The accompanying frequency map shows
whether exceedance is confined to one projection or recurs across the series.

The largest increase in Range Mean Relative Niche Boundary Exceedance occurs
between 2070 and 2090:

```{r rate-table, echo = FALSE}
rate_row <- subset(
  change_rate,
  metric == "range_wide_relative_exceedance"
)
rate_table <- data.frame(
  `Interval start` = rate_row[["interval_start"]],
  `Interval end` = rate_row[["interval_end"]],
  `Increase in range mean relative exceedance` = paste0(
    round(100 * rate_row[["maximum_interval_increase"]], 1),
    " percentage points"
  ),
  check.names = FALSE
)

knitr::kable(rate_table, row.names = FALSE)
```

## Report

`climniche_series_report()` collects the fixed niche reference, projections,
range summaries, temporal Niche Boundary Exceedance summaries and interval
changes.

```{r series-report, eval = FALSE}
series_report <- climniche_series_report(
  series,
  species = "European anchovy",
  scope = "current",
  area_weight = TRUE
)

series_report
write_climniche_series_report(
  series_report,
  "anchovy-exposure-through-time.md"
)
```