study-area.qmd

# Study area and trenching data {#sec-area}

In this project we will use data from Levanger in central Norway. This area was part of the pilot study to identify mire trenches in Norway based on LiDAR [@jansson2024] (see also @sec-background).

```{r setup}
#| warning: false
#| error: false
#| include: false
library(sf)
library(tidyverse)
library(basemaps)
library(stars)
library(ggpubr)

options(scipen = 999)
dir <- substr(getwd(), 1,2)
# conditional path
path1 <- ifelse(dir == "C:", "R:/", "/data/R/")
  
```

```{r data}
#| code-fold: true
#| code-summary: "Read shape file and convert to gpkg" 
#| eval: false

path <- "41201785_okologisk_tilstand_2022_2023/data/grofter/Kolstad"
lev <- sf::read_sf(paste0(path1, path, "/LevangerGroftPolygon.shp")) |>
  sf::st_make_valid()
lev |> sf::st_write("data/lev.gpkg")
```

```{r data2}
#| code-fold: true
#| code-summary: "Read data"
lev <- sf::read_sf("data/lev.gpkg")
```

```{r stats}
#| code-fold: true
#| code-summary: "Calculate study extent"

lev_bb <- sf::st_bbox(lev)
myx <- unname(round(lev_bb["xmax"] - lev_bb["xmin"]))
myy <- unname(round(lev_bb["ymax"] - lev_bb["ymin"]))
```

The study area is `r myx` 𝖷 `r myy` meters, or `r round(myx*myy/10^6)` km^2^.

```{r basemapSetup}
#| code-fold: true
#| code-summary: "Set up basemap"

basemaps::set_defaults(map_service = "osm", map_type = "topographic")
lev_t <- lev |>
  st_transform(3857)
```

```{r fig-map1}
#| code-fold: true
#| code-summary: "Make map"
#| fig-cap: |
#|   "Map of study area in Levanger, central Norway, showing trenched mires i purple. 
#|   Note that the colored polygon borders exagerrate the area of the trenches."
#| cache: true
#| message: false
#| warning: false

ggplot() + 
  basemap_gglayer(sf::st_bbox(lev_t)) +
  scale_fill_identity() + 
  coord_sf() +
  geom_sf(data = lev_t,
          fill = "purple",
          color = "purple") +
  theme_bw()
```

Let's zoom in on an area just north of Tomtvatnet, and overlay the trenches data with a map of open mires [@bakkestuen2023].

```{r zoombbox}
#| code-fold: true
#| code-summary: "Creat new bbox"
#| warning: false

bbox_new <- st_bbox(lev) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

bbox_new[1] <- bbox_new[1] + (0.3 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] - (0.5 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] + (0.4 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] - (0.5 * yrange) # ymax - top

lev_crop <- lev |>
  sf::st_crop(bbox_new)
```

```{r readMire}
#| code-fold: true
#| code-summary: "Read mire data"
#| eval: false

# I took the gdb file in the same folder and opened it in arcgis pro
# and exported a section of it as geoTIFF
# Raster GBD is not readable in R with GDAL <3.7
path <- "41201785_okologisk_tilstand_2022_2023/data/Myrmodell/myrmodell_levanger.tif"
mire90 <- stars::read_stars(paste0(path1, path)) |>
  st_crop(lev_bb) |> 
  setNames("prob") |>
  mutate(prob = case_when(
    prob > 900 ~ 1,
    TRUE ~ NA
  ))
stars::write_stars(mire90, "data/mire_levanger_reclassified90.tif")
saveRDS(mire90, "data/mire_levanger_reclassified90.rds")

mire70 <- stars::read_stars(path) |>
  st_crop(lev_bb) |> 
  setNames("prob") |>
  mutate(prob = case_when(
    prob > 700 ~ 1,
    TRUE ~ NA
  ))
stars::write_stars(mire70, "data/mire_levanger_reclassified70.tif")
saveRDS(mire70, "data/mire_levanger_reclassified70.rds")

mire50 <- stars::read_stars(path) |>
  st_crop(lev_bb) |> 
  setNames("prob") |>
  mutate(prob = case_when(
    prob > 500 ~ 1,
    TRUE ~ NA
  ))
stars::write_stars(mire50, "data/mire_levanger_reclassified50.tif")
saveRDS(mire50, "data/mire_levanger_reclassified50.rds")

mire10 <- stars::read_stars(path) |>
  st_crop(lev_bb) |> 
  setNames("prob") |>
  mutate(prob = case_when(
    prob > 100 ~ 1,
    TRUE ~ NA
  ))
stars::write_stars(mire10, "data/mire_levanger_reclassified10.tif")
saveRDS(mire10, "data/mire_levanger_reclassified10.rds")

```

```{r}
#| code-fold: true
#| code-summary: "Read mire data"

#90% probability
mire90 <- readRDS("data/mire_levanger_reclassified90.rds")
mire90_crop <- mire90 |>
  st_crop(st_bbox(lev_crop))

#10% probability
mire10 <- readRDS("data/mire_levanger_reclassified10.rds")
mire10_crop <- mire10 |>
  st_crop(st_bbox(lev_crop))


```

```{r fig-map2}
#| code-fold: true
#| code-summary: "Create map"
#| cache: true
#| warning: false
#| fig-cap: |
#|   "Zoomed in version of the previous map showing systematic trenching,
#|   but also more sporadic trenches that might be errors in the model.
#|   The red areas are open mires with 10% probability.
#|   Yellow areas are mires with 90% probability. From areal photos, the 10% probability seem to match best.
#|   Note in any case that most of the trenches are located in 
#|   other nature types, mosty forests."


ggplot() + 
  geom_sf(data = lev_crop,
          fill = "purple",
          color = "purple") +
  geom_stars(data = mire10_crop,
             downsample = 0,
             fill = "red",
             na.action = na.omit,
             alpha=.7)+
  geom_stars(data = mire90_crop,
             downsample = 0,
             fill = "yellow",
             na.action = na.omit,
             alpha=.7)+
  theme_bw()

```

The trenches model originally results in a raster with a continuous probability value for the occurrence of a trench in each 1 𝖷 1 m cell. What I show above (@fig-map1, @fig-map2) is a version of this data where a threshold value of 90% for the probabilities turn it into presence-absence data in the shape of polygons.


## Adding AR5

The open wetlands, shown as red areas in @fig-map2, are largely laying outside of the systematically trenched areas. This can be because the trenched areas, which may have been open mire to begin with, now are overgrown with trees. The map of open mires is in way trained to find areas that are not trenched, so no wonder if the ecosystem condition looks good in these areas. Instead I think we should include forested mires in the Ecosystem Type (ET). We can use the mire class in the AR5 maps in addition to the open mires map. AR5 underestimates the total mire area, but it probably includes more forested mires. Another alternative would be to use the new Ecosystem Map of Norway [@strand2023], but this uses AR50, i.e. as coarser generalization of AR5.

```{r ar5}
#| code-fold: true
#| code-summary: "Prepare AR5 data"
#| eval: false
#| include: false

ar5_path <- "GeoSpatialData/Topography/Norway_FKB/Original/FKB-AR5 FGDB-format/Basisdata_50_Trondelag_25832_FKB-AR5_FGDB.gdb"
st_layers(paste0(path1, ar5_path))
# this takes a while and should be done on the server
ar5 <- sf::read_sf(paste0(path1, ar5_path),
                   layer = "fkb_ar5_omrade")
ar5_crop <- ar5 |>
  st_crop(bbox_new)

ar5_crop <- ar5_crop |>
  mutate(Arealtype = case_match(
    arealtype,
    "60" ~ "Myr",
    "12" ~ "Samferdsel",
    "21" ~ "Fulldyrka mark",
    "23" ~ "Innmarksbeite",
    "30" ~ "Barskog",
    "50" ~ "Åpen mark",
    "81" ~ "Ferskvann",
    .default = arealtype
    
  ))

ar5_union <- ar5_crop |>
  group_by(Arealtype) |>
  summarise()

saveRDS(ar5_union, "data/ar5.rds")
```

```{r fig-ar5}
#| fig-cap: |
#|   "AR5 overlayed with the mire trenches data. 
#|   AR5 classifies most of the systamtically trenched peatlands as mire,
#|   and some as forest.
#| code-fold: true
#| code-summary: "Read AR5 and plot"


ar5 <- readRDS("data/ar5.rds")
ggplot() + 
  geom_sf(data = ar5,
          aes(fill = Arealtype)) +
  geom_sf(data = lev_crop,
          fill = "yellow",
          color = "yellow")
  
```

Let's separate out the areas that are mire in AR5, but not in the open mire map,
and see if adding these areas make a difference. Of coarse, we are only looking at
a very small spatial extent here, and it is not representative for all the 
differences between AR5 and the opne mire map, but still.

```{r makeDatasetsAndPlots}
#| code-fold: true
#| code-summary: "Filter data"
#| warning: false

  
# make the ecosystem delineation map into polygons
mire <- st_as_sf(mire10_crop, merge = T)
mire_gg <- ggplot(data = mire)+
  geom_sf(fill = "blue")+
  labs(title = "Open mire model")

# Filter out the mire part of AR5
ar5_myr <- ar5 |>
  filter(Arealtype == "Myr") |>
  st_cast("POLYGON")
ar5_myr_gg <- ggplot(data = ar5_myr)+
  geom_sf(fill = "cyan4")+
  labs(title = "AR5")

# Take out the open mire part of AR5 (leaving the forested peatlands)
st_erase = function(x, y) st_difference(x, st_union(st_combine(y)))
ar5_nonOpen <- ar5_myr |>
  st_erase(mire)

# Take out the AR5 part of the open mire model
only_open <- mire |>
  st_erase(ar5_myr) |>
  st_cast("MULTIPOLYGON") |>
  st_cast("POLYGON")

only_open_gg <- ggplot(data = only_open)+
  geom_sf(fill = "purple")+
  labs(title = "Open mires not in AR5")

ar5_nonOpen_gg <- 
  ggplot(data = ar5_nonOpen)+
  geom_sf(fill = "orange")+
  labs(title = "AR5 minus open mires")

# check that polygons dont overlap
check <- 
  ggplot(data = ar5_nonOpen)+
  geom_sf(fill = "red",
          alpha=.5)+
  geom_sf(data = mire,
          fill = "blue",
          alpha=.5)+
  labs(title = "AR5 (red) and open mires\n (blue) combined")

# combine into new mire map
# This first attempt with st_union produces way to many small polygon
#peatlands <- ar5_myr |>
#  st_union(mire) |>
#  st_cast("MULTIPOLYGON") |>
#  st_cast("POLYGON")

#but this works:
peatlands <- ar5_myr |>
  st_join(mire) |>
  st_intersection()

peatlands_gg <- ggplot(data = peatlands)+
  geom_sf(fill = "green")+
  labs(title = "AR5 + open mires")
```




```{r fig-mireExtents}
#| code-fold: true
#| code-summary: "Plot data"
#| fig-cap: |
#|   "Different combinations of the open mire model
#|   and the AR5 mire category."

ggarrange(mire_gg, 
          ar5_myr_gg, 
          ar5_nonOpen_gg, 
          only_open_gg,
          check, 
          peatlands_gg,
          ncol=2, nrow=3)

```


```{r savefiles}
#| include: false
#| eval: false
save(
  mire10_crop,
  lev_crop,
  ar5_nonOpen,
  peatlands,
  file = "data/study-area-files.RData"
)

```

## Ecosystem delineation {#sec-delineation}

We presented a study area that is practically square (@fig-map1). 
However, for the purpose of an ecosystem condition account or assessment, 
this area needs to be partitioned into unique and non-overlapping ecosystems,
and our indicator needs to address the ecosystem condition the target ecosystem explicitly.
The ecosystem that this indicator is aimed at is wetlands, or _våtmark_ in Norwegian.


We might therefore need to mask our spatial data at some stage in order to exclude 
trenches identified in other ecosystems. 
However, this is problematic for several reasons.

We do not have a very precise ecosystem delineation map for wetlands,
and it is not yet clear if the ecosystem should be defined as all types of wetlands,
or just open wetlands. In the case of the latter, we have a quite good ecosystem map that we can use (@fig-map2). But by excluding everything that is not open mire,
we automatically exclude wetlands, or previous wetlands, that 
have become forested due to the lowered water table that comes from trenching.
In other words, we exclude the areas that are in the worst condition.
From @fig-map2 we can see that this is a very common scenario.
Alternatively, we can use things like AR5, or a combination of the two (@fig-mireExtents).


When it comes to updating the indicator in the future (see @sec-update),
we will be in the case that more and more trenched wetland becomes classified as other
nature types, mainly forests, because that is what happens when you trench a mire. 
As a consequence, the indicator values will look better and better over time,
reflection the land use change and the associated loss of the areas in worst condition. This is a well known issue in ecosystem accounting (ref Jacobson).

We recommend for this indicator to delineate the nature type as open mires in the open mire model [@bakkestuen2023], in combination with mires in AR5 (@fig-mireExtents; last pane). 



# References {.unnumbered}