7TK.qmd

# Frequency-based metric {#sec-7TK}

Here I will exemplify the frequency based metric using a method similar to the `7TK` variable in NiN.

```{r loadData}
#| include: false
load("data/study-area-files.RData")
library(sf)
library(tidyverse)
library(units)
library(stars)
library(eaTools)
library(ggridges)
```

Making a 100x100m grid over the stydy area (the smaller zoomed in area from @fig-map2).

```{r fig-makeGrid}
#| code-fold: true
#| code-summary: "Make grid"
#| fig-cap: "100x100m grid over the stydu area bbox, with random grid cell IDs."

grid100 <- sf::st_make_grid(lev_crop,
                 cellsize = 100) |>
  st_sf()
  
id <- sample(1:nrow(grid100))
grid100 <- grid100 |>
  add_column(id = id)

plot(grid100)

```

Then I intersect the trenches data with the grid.

<!--# https://vialab.mit.edu/tutorials/module/mapping-in-r-street-trees-in-camberville/ -->

```{r fig-trenchInCel}
#| code-fold: true
#| code-summary: "Intersect trenches and grid cells"
#| fig-cap: |
#|   "Study area gridded with 100x100 m cells and trenches labled according to grid ID"
#| warning: false

trench_in_cell <- st_intersection(lev_crop, grid100)
#trench_in_cell <- st_join(lev_crop, grid100, join=st_within)
#saveRDS(trench_in_cell, "data/trench_in_cell.RDS")
trench_in_cell |>
  ggplot()+
  geom_sf(data = grid100)+
  geom_sf(aes(color = factor(id),
              fill = factor(id)))+
  guides(fill = "none",
         color = "none") +
  theme_bw()
```

Then I collapse the dataframe into unique grid cell IDs and mark these as 1's (presenses).

I then join this data with `grid100` and plot it.

```{r fig-trencharea}
#| code-fold: true
#| code-summary: "Show code"
#| fig-cap: |
#|   "Map of study area showing 100x100m grid cells classed as either trenched
#|   (blue) or not trenched (red)."

grid100_pa <- trench_in_cell |>
  st_set_geometry(NULL) |>
  group_by(id) |>
  summarise(n = n()) |>
  add_column(trenched = 1) |>
  select(-n) |>
  right_join(grid100, by = "id") |>
  mutate(trenched = case_when(
    is.na(trenched) ~ 0,
    .default = trenched
  )) |>
  st_as_sf()

grid100_pa |>
  ggplot() +
  geom_sf(aes(fill = factor(trenched)),
          na.rm = TRUE)+
  theme_bw() +
  geom_sf(data = trench_in_cell)
```

The example in @fig-trencharea is done on 100x100 grid cell just for illustration. I will now do the same, but on 10x10m grid cells, and then calculate the proportion of trenched 10x10 grid cells inside each 100x100m grid cell. The proportion of 10x10 grid cells could be aggregated to any geometry, such as a region, a municipality or a project area.

```{r fig-trenched}
#| code-fold: true
#| code-summary: "Show code"
#| warning: false
#| fig-cap: |
#|   "Map of study area showing the 100x100 grid cells that are colored
#|   based on the percentage of 10x10m grid cells that have a trench in them."


grid10 <- sf::st_make_grid(grid100,
                cellsize = 10) |>
  st_sf()

# align perfectly with grid100 and carry the id from grid100 to grid10
grid10 <- grid10 |> 
  st_intersection(grid100) |>
  mutate(
    area = units::drop_units(
      st_area(sf..st_make_grid.grid100..cellsize...10.)),
    exclude = case_when(
      area < 50 ~ 1,
      .default = 0
         )) |>
  filter(exclude != 1)


id10 <- sample(1:nrow(grid10))
grid10 <- grid10 |>
  add_column(id10 = id10) |>
  select(id, id10)

#plot(grid10)

trench_in_cell10 <- st_intersection(lev_crop, grid10)

grid100_freq <- trench_in_cell10 |>
  st_set_geometry(NULL) |>
  group_by(id10) |>
  summarise(n = n()) |>
  ungroup() |>
  add_column(trenched = 1) |>
  select(-n) |>
  right_join(grid10, by = "id10") |>
  mutate(trenched = case_when(
    is.na(trenched) ~ 0,
    .default = trenched
  )) |>
  group_by(id) |>
  summarise(freq = mean(trenched)*100) |>
  #st_as_sf() |>
  #st_intersection(grid100) |>
  right_join(grid100, by = "id") |>
  st_as_sf()

grid100_freq |>
  ggplot() +
  geom_sf(aes(fill = freq))+
  theme_bw() +
  geom_sf(data = trench_in_cell)
```

@fig-trenched shows the trenching variable in its raw units (%) over the entire study area. The next step is then to exclude areas that are not part of the target ecosystem, i.e. not part of the accounting area. This can be done on the 100x100 grid level, put probably more correct to do it at the 10x10m grid level. To start I mask againt open mires only. Later I add AR5 mires.

```{r fig-trenched_mask}
#| code-fold: true
#| code-summary: "Show code"
#| warning: false
#| fig-cap: |
#|   "Map of study area (masked against mire extent) 
#|   showing 100x100m grid cells colored by the proportion of 
#|   10x10 grid cells inside them which are trenched. "

# make the ecosystem delineation map into polygons
mire <- st_as_sf(mire10_crop, merge = T)

grid100_freq2 <- trench_in_cell10 |>
  st_intersection(mire) |>
  st_drop_geometry() |>
  group_by(id10) |>
  summarise(n = n()) |>
  ungroup() |>
  add_column(trenched = 1) |>
  select(-n) |>
  right_join(grid10, by = "id10") |>
  mutate(trenched = case_when(
    is.na(trenched) ~ 0,
    .default = trenched
  )) |>
  st_as_sf() |>
  st_intersection(mire) |>
  st_drop_geometry() |>
  group_by(id) |>
  summarise(freq = mean(trenched)*100) |>
  #st_as_sf() |>
  #st_intersection(grid100) |>
  right_join(grid100, by = "id") |>
  st_as_sf() |>
  st_intersection(mire)

high <- "red"
low <- "green"


grid100_freq2 |>
  ggplot() +
  geom_sf(aes(fill = freq))+
  theme_bw() +
  geom_sf(data = trench_in_cell)+
  scale_fill_gradient("Frequency of trenches", low = low, 
        high = high)
  
```

```{r}
#| include: false
highest <- grid100_freq2 |>
  arrange(desc(freq)) |>
  slice_head(n = 3) |>
  pull(freq) |>
  round()
```

The three grid cells with the highest proportion of trenched 10x10m grid cells have
`r highest[1]`, `r highest[2]`, and `r highest[3]` percent.
The three cells visible in @fig-trenched_mask with the darkest red color
all seem to be relatively densely trenched.
Probably we need to adjust the threshold value so that at least these three all 
get quite low indicator values.
I will for now set 10% as the threshold value (X~60~) and 50% as the low reference point (X~0~). These values are updated later (see @sec-ref and @sec-ref2).
I will fit the data using a sigmoid function to soften arbitrary break points.

```{r fig-normalise}
#| code-fold: true
#| warning: false
#| code-summary: "Show code"
#| fig-cap: |
#|   "Normalisation of the variable 'Porportion of small grid cells with 
#|   trenching' based on a sigmoid transformation function and thee numerical 
#|   reference points."

ea_normalise(
  data = grid100_freq2,
  vector = "freq",
  upper_reference_level = 50,
  lower_reference_level = 0,
  break_point = 10,
  scaling_function = "sigmoid",
  reverse=T,
  plot = T)
```

```{r fig-normalisemap}
#| code-fold: true
#| code-summary: "Adding the normalised data"
#| fig-cap: |
#|   "Same figure as the previous map, but with the trenching variable 
#|   normalised into indicator values."

grid100_freq2 <- grid100_freq2 |>
  mutate(indicator = 
           ea_normalise(
  data = grid100_freq2,
  vector = "freq",
  upper_reference_level = 50,
  lower_reference_level = 0,
  break_point = 10,
  scaling_function = "sigmoid",
  reverse=T,
  plot = F))

grid100_freq2 |>
  ggplot() +
  geom_sf(aes(fill = indicator))+
  theme_bw() +
  geom_sf(data = trench_in_cell)+
  scale_fill_gradient(
    "Indicator value", 
    low = high, 
    high = low)
```

I think the indicator values are reasonably well representative of what I 
expect the ecological effect of the different trenching densities to be.
These reference points are, however, not informed with any data.

```{r indicatorValue}
#| code-fold: true
#| code-summary: "Calculate indicator values"

grid100_freq2 <- grid100_freq2 |>
  mutate(area = st_area(sf..st_make_grid.lev_crop..cellsize...100.))

ind_dist_open <- NULL
for(i in 1:1000){
  ind_dist_open[i] <- 
    mean(sample(grid100_freq2$indicator, 
       nrow(grid100_freq2),
       replace = TRUE,
       prob = grid100_freq2$area))
}

#ind_dist_open |>
#  as_tibble() |>
#  ggplot()+
#  geom_density(aes(x = value),
#    fill = "burlywood") +
#  labs(x = "Indicator value")
```



## Including forest peatlands

Now let's calculate the indicator value based on the orange areas (the third map)
and the green areas (the last map) in @fig-mireExtents,
and see if it looks very different what we had.

```{r fig-indicatorNonOpen}
#| code-fold: true
#| code-summary: "Calculate indicator values for non-open mires in AR5"
#| warning: false
#| fig-cap: |
#|   "Indicator values over non-open AR5 mires only"

grid100_freq3 <- trench_in_cell10 |>
  st_intersection(ar5_nonOpen) |>
  st_drop_geometry() |>
  group_by(id10) |>
  summarise(n = n()) |>
  ungroup() |>
  add_column(trenched = 1) |>
  select(-n) |>
  right_join(grid10, by = "id10") |>
  mutate(trenched = case_when(
    is.na(trenched) ~ 0,
    .default = trenched
  )) |>
  st_as_sf() |>
  st_intersection(ar5_nonOpen) |>
  st_drop_geometry() |>
  group_by(id) |>
  summarise(freq = mean(trenched)*100) |>
  right_join(grid100, by = "id") |>
  st_as_sf() |>
  st_intersection(ar5_nonOpen)

grid100_freq3 <- grid100_freq3 |>
  mutate(indicator = 
           ea_normalise(
  data = grid100_freq3,
  vector = "freq",
  upper_reference_level = 50,
  lower_reference_level = 0,
  break_point = 10,
  scaling_function = "sigmoid",
  reverse=T,
  plot = F))

grid100_freq3 |>
  ggplot() +
  geom_sf(aes(fill = indicator))+
  theme_bw() +
  geom_sf(data = trench_in_cell)+
  scale_fill_gradient(
    "Indicator value", 
    low = high, 
    high = low)
```

@fig-indicatorNonOpen is a lot more red than @fig-normalisemap!

```{r indicatorValue_nonopen}
#| code-fold: true
#| code-summary: "Calculate indicator values for non-open mire"

grid100_freq3 <- grid100_freq3 |>
  mutate(area = st_area(sf..st_make_grid.lev_crop..cellsize...100.))

ind_dist_nonopen <- NULL
for(i in 1:1000){
  ind_dist_nonopen[i] <- 
    mean(sample(grid100_freq3$indicator, 
       nrow(grid100_freq3),
       replace = TRUE,
       prob = grid100_freq3$area))
}

#ind_dist_nonopen |>
#  as_tibble() |>
#  ggplot()+
#  geom_density(aes(x = value),
#    fill = "burlywood") +
#  labs(x = "Indicator value")
```

Let's do the same, but include all AR5 mires, and the little bit extra areas 
from the open mire map.

```{r fig-indicatorAll}
#| code-fold: true
#| code-summary: "Calculate indicator values for open and non-open mires"
#| warning: false
#| fig-cap: |
#|   "Indicator values over all peatlands (AR5 + the open mire model)"
#| cache: true

grid100_freq4 <- trench_in_cell10 |>
  st_intersection(peatlands) |>
  st_drop_geometry() |>
  group_by(id10) |>
  summarise(n = n()) |>
  ungroup() |>
  add_column(trenched = 1) |>
  select(-n) |>
  right_join(grid10, by = "id10") |>
  mutate(trenched = case_when(
    is.na(trenched) ~ 0,
    .default = trenched
  )) |>
  st_as_sf() |>
  st_intersection(peatlands) |>
  st_drop_geometry() |>
  group_by(id) |>
  summarise(freq = mean(trenched)*100) |>
  right_join(grid100, by = "id") |>
  st_as_sf() |>
  st_intersection(peatlands)

grid100_freq4 <- grid100_freq4 |>
  mutate(indicator = 
           ea_normalise(
  data = grid100_freq4,
  vector = "freq",
  upper_reference_level = 50,
  lower_reference_level = 0,
  break_point = 10,
  scaling_function = "sigmoid",
  reverse=T,
  plot = F))

# Option to have uniform class:
grid100_freq4 <- grid100_freq4 |>
  st_cast("MULTIPOLYGON") |>
  st_cast("POLYGON")

#saveRDS(grid100_freq4, "data/grid100_freq4.rds")

grid100_freq4 |>
  ggplot() +
  geom_sf(aes(fill = indicator))+
  theme_bw() +
  geom_sf(data = trench_in_cell)+
  scale_fill_gradient(
    "Indicator value", 
    low = high, 
    high = low)
```

```{r indicatorValue_all}
#| code-fold: true
#| code-summary: "Calculate indicator values for all mires"

grid100_freq4 <- grid100_freq4 |>
  mutate(area = st_area(sf..st_make_grid.lev_crop..cellsize...100.))

ind_dist_all <- NULL
for(i in 1:1000){
  ind_dist_all[i] <- 
    mean(sample(grid100_freq4$indicator, 
       nrow(grid100_freq4),
       replace = TRUE,
       prob = grid100_freq4$area))
}

ind_dist_comb <- tibble(
  open = ind_dist_open,
  forested = ind_dist_nonopen,
  all = ind_dist_all
) |>
  pivot_longer(
    cols = everything(),
    names_to = "extent",
    values_to = "indicator_value")
```

```{r fig-ridge}
#| code-fold: true
#| code-summary: "Make ridge plot"
#| fig-cap: |
#|   "Indicator values (distributions from bootstrapping)
#|   for the mire trenching indicator cacualted for a small test area
#|   and using three different ecosystem delineation maps. 
#|   Vertical lines are 2.5, 50 and 97.5 percentiles."
#| warning: false


ind_dist_comb |>
  ggplot(
    aes(
      x = indicator_value, 
      y = extent,
      fill = after_stat(x))
  ) +
  geom_density_ridges_gradient(
    bandwidth=.005,
    quantile_lines = TRUE, 
    quantiles = c(0.025, .5, 0.975)
  ) +
  scale_fill_distiller(palette = "RdYlGn",
    direction = 1) +
  theme_bw(base_size = 15) +
  guides(fill = "none") +
  geom_vline(xintercept = .6,
    size=1.2,
    lty=2) +
  theme(axis.text.x = element_text(angle=45, hjust=1, vjust=1))+
  labs(y = "", x = "Indicator values")+
  xlim(c(0,1)
)


```

@fig-ridge shows us that the ecosystem delineation makes a huge difference
on the indicator values. Combining AR5 and the open mire model gives a mean indicator 
value of `r round(mean(ind_dist_all), 2)`, reflection quite predictably the 
rather binary situation in @fig-indicatorAll where half the area is dark red
and the other half is green.

I think it makes most sense to combine AR5 and
the open mire model for this indicator. All indicators used in the
same assessment must use the same overall ecosystem delineation, but it is still possible that one indicator related only to a smaller part of that delineation. 
If using the AR5+open mire combination is not agreeable, then using AR5, alternatively the new Ecosystem Map of Norway [@strand2023], would be preferable to using just the 
open mire model.