README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  echo = F,
  warning = F,
  message = F,
  error = F
)
library(sf)
library(terra)
library(geodata)
library(tidyterra)
library(ggplot2)
library(dplyr)
library(SDMWorkflows)
library(rgbif)
library(tidyverse)
library(sars)
```

# Klimaskovfond

<!-- badges: start -->
<!-- badges: end -->

The goal of Klimaskovfond is to generate the needed data sets to get analize the Klimaskovfond score system


```{r, cache=T}
Forest <- terra::rast("o:/Nat_BDR-data/Arealanalyse/2023/CLEAN/Rasterized/Rast_AllForerst_Croped.tif") |> as.numeric()

Desciduous <- terra::rast("o:/Nat_Ecoinformatics-tmp/JakobAssmann_au634851/treetype/treetype_bjer_dec.tif") 

Conifer <- terra::rast("o:/Nat_Ecoinformatics-tmp/JakobAssmann_au634851/treetype/treetype_bjer_con.tif") 


```
```{r asPolygons, cache=T}
Forest[Conifer== 1] <- NA

Forest_Poly <-  as.polygons(Forest) |> aggregate(by = "C_02") |> 
  disagg()
Forest_Poly$ID <- 1:nrow(Forest_Poly)
terra::writeVector(Forest_Poly, "Forest.shp", overwrite = T)
sf::write_sf(st_as_sf(Forest_Poly), "Forest.geojson")
```


```{r, echo = F, eval = F}
library(doParallel)
library(foreach)
require(doSNOW)
library(data.table)

cl <- makeSOCKcluster(50)

Forest_Poly <- terra::vect("Forest_Area.shp")

MoreThan2 <- Forest_Poly[Forest_Poly$Ha > 2,]

MoreThan2DF <- as.data.frame(MoreThan2)
MoreThan2DF$Richness <- NA

terra::writeVector(MoreThan2, "MoreThan2.shp")

doSNOW::registerDoSNOW(cl)

pb <- txtProgressBar(min=1, max=nrow(MoreThan2DF), style=3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress=progress)
#v
Species <- foreach(i=1:nrow(MoreThan2DF), .packages = c("terra", "sf", "rgbif", "SDMWorkflows", "dplyr"), .options.snow=opts)%dopar% {
  try({
    MoreThan2 <- terra::vect("MoreThan2.shp")
    WKT <- MoreThan2[i,]|>
      terra::convHull() |>
      terra::project("epsg:4326") |>
      st_as_sf() |>
      st_as_sfc() |>
      st_as_text()

    Test <- rgbif::occ_count(hasCoordinate = T,
                             geometry = WKT,
                             year = '1999,2023',
                             facet="scientificName", facetLimit=10000)

    Temp <- SDMWorkflows::Clean_Taxa(Taxons = Test$scientificName) |>
      dplyr::select(kingdom, phylum, class, family, genus, species) |>
      dplyr::distinct()

    Temp$ID <- i
    MoreThan2DF[i,]$Richness <- nrow(Temp)
    message(paste(i, "of", nrow(MoreThan2DF) ,"ready", Sys.time()))
    Temp

  })
}

readr::write_csv(MoreThan2DF, "RichnessParallel.csv")
stopCluster(cl)

Species <- Species |> purrr::keep(is.data.frame) |> purrr::map(data.table::as.data.table) |> data.table::rbindlist()

saveRDS(Species, "Species.rds")

MoreThan2DF <- as.data.frame(MoreThan2) |> dplyr::select(Ha)
MoreThan2DF$ID <- 1:nrow(MoreThan2DF)

MoreThan2DF <- MoreThan2DF[MoreThan2DF$ID %in% unique(Species$ID),]

Species_Animals_Plants <- Species[kingdom %chin% c("Animalia", "Plantae")]

Species_Animals_Plants_N <- Species_Animals_Plants[, .N, by = ID]

Species_Animals_Plants_N <- Species_Animals_Plants_N |> 
  dplyr::left_join(MoreThan2DF)

write_csv(Species_Animals_Plants_N, "Animals_plants_richness.csv")

Species_Plants <- Species[kingdom %chin% c("Plantae")]

Species_Plants_N <- Species_Plants[, .N, by = ID]

Species_Plants_N <- Species_Plants_N |> 
  dplyr::left_join(MoreThan2DF)

write_csv(Species_Plants_N, "Plants_richness.csv")

Species_Animals <- Species[kingdom %chin% c("Animalia")]

Species_Animals_N <- Species_Animals[, .N, by = ID]

Species_Animals_N <- Species_Animals_N |> 
  dplyr::left_join(MoreThan2DF)

write_csv(Species_Animals_N, "Animals_richness.csv")
```


```{r HaArea, cache=T}
Forest_Poly <- terra::vect("Forest.shp")
Forest_Poly$Ha <- terra::expanse(Forest_Poly, unit= "ha") |> round(3)
terra::writeVector(Forest_Poly, "Forest_Area.shp", overwrite = T)
sf::write_sf(st_as_sf(Forest_Poly), "Forest_Area.geojson")

Forest_DF <- as.data.frame(Forest_Poly) |> dplyr::mutate(LogHa = log(Ha))
```

## Overview of Forest Patches

We have a total of `r prettyNum(nrow(Forest_DF), big.mark = ",")` forest patches, with sizes ranging from `r prettyNum(round(min(Forest_DF$Ha), 3), big.mark = ",")` to `r prettyNum(max(Forest_DF$Ha), big.mark = ",")` hectares. The distribution of forest patch areas is visualized below:

```{r histogramplot}
ggplot(Forest_DF, aes(x = Ha)) + geom_histogram(bins = nclass.Sturges(Forest_DF$LogHa)) + scale_x_log10(labels = scales::comma) + scale_y_continuous(labels = scales::comma)
```

Next, let's focus on forest patches larger than 2 hectares.

```{r, echo = F, cache=T}
Forest_Poly <- terra::vect("Forest_Area.shp")

MoreThan2 <- Forest_Poly[Forest_Poly$Ha > 2,]

MoreThan2DF <- as.data.frame(MoreThan2)

MoreThan2DF <- MoreThan2DF |> dplyr::mutate(LogHa = log(Ha))
```


Out of the total forest patches, only `r prettyNum(nrow(MoreThan2DF), big.mark = ",")` correspond to patches larger than 2 hectares. The distribution of these larger patches is visualized below:

```{r histogramplot2}
ggplot(MoreThan2DF, aes(x = Ha)) + geom_histogram(bins = nclass.Sturges(Forest_DF$LogHa)) + scale_x_log10(labels = scales::comma) + scale_y_continuous(labels = scales::comma)
```

## Species Area Relationship Analysis

To determine the optimal forest patch size for biodiversity preservation, we conducted a Species Area Relationship analysis using the sars R package. The analysis utilized GBIF data on species presences within forest polygons from 1999 to 2023, resolving synonyms.

The distribution of species across kingdoms is summarized below: 

```{r tablePercentage}
Species <- readRDS("Species.rds")

Porp <- Species$kingdom |> table() |> prop.table()
DF <- data.frame(Kingdom = names(Porp), Percentage = round(as.numeric(Porp)*100, 2)) |> dplyr::arrange(desc(Percentage))

knitr::kable(DF)
```

Additionally, we present the top 10 classes with the highest proportions:

```{r tablePercentage2, echo = F}
DF2 <- Species |> 
  group_by(kingdom, phylum, class) |> 
  summarise(n = n()) |> 
  arrange(desc(n))

Total <- sum(DF2$n)

DF2 <- DF2 |> 
  dplyr::mutate(Percentage = round(100*(n/Total), 2)) |> 
  ungroup() |> 
  dplyr::slice_max(order_by = Percentage, n = 10) |> 
  dplyr::select(-n) |> 
  dplyr::mutate(Cumulative_Percentage = cumsum(Percentage))

knitr::kable(DF2)
```

# Species Area Relationship for Different Groups

Now, we perform a species area relationship analysis for three groups: plants and animals together, plants only, and animals only.

Plants and Animals


```{r mysarsfuncc, echo = F, cache=T}
MySarsPlot <- function(x, type = "multi", allCurves = FALSE,
                       pch = 16, cex = 1.2, pcol = "dodgerblue2", ModTitle = NULL,
                       TiAdj = 0, TiLine = 0.5, cex.main = 1.5, cex.lab = 1.3, cex.axis = 1,
                       yRange = NULL, lwd = 2, lcol = "dodgerblue2", mmSep = FALSE,
                       lwd.Sep = 6, col.Sep = "black", pLeg = TRUE, modNames = NULL,
                       cex.names = 0.88, subset_weights = NULL, confInt = FALSE)
{
  if (confInt) {
    if (length(x$details$confInt) == 1)
      stop("No confidence interval information in the fit object")
    CI <- x$details$confInt
  }
  ic <- x[[2]]$ic
  dat <- x$details$fits
  dat2 <- dat
  df <- dat[[1]]$data
  xx <- df$A
  yy <- df$S
  nams <- vapply(dat2, function(x) x$model$name, FUN.VALUE = character(1))
  xx2 <- seq(min(xx), max(xx), length.out = 1000)
  mf <- lapply(dat2, function(y) {
    if (y$model$name == "Linear model") {
      c1 <- y$par[1]
      m <- y$par[2]
      ff <- c1 + (m * xx2)
    }
    else {
      ff <- y$model$mod.fun(xx2, y$par)
      if (anyNA(ff)) {
        stop("Error in plotting, contact package owner")
      }
    }
    ff
  })
  mf2 <- matrix(unlist(mf), ncol = length(dat2), byrow = FALSE)
  mf2 <- as.data.frame(mf2)
  colnames(mf2) <- nams
  icv <- vapply(dat2, function(x) unlist(x[[ic]]), FUN.VALUE = double(1))
  delt <- icv - min(icv)
  akaikesum <- sum(exp(-0.5 * (delt)))
  aw <- exp(-0.5 * delt)/akaikesum
  if (round(sum(aw), 0) != 1)
    stop("IC weights do not sum to 1")
  mf3 <- matrix(nrow = nrow(mf2), ncol = ncol(mf2))
  for (i in seq_along(aw)) mf3[, i] <- mf2[, i] * aw[i]
  wfv <- rowSums(mf3)

  nams2 <- nams
  DF2 <- data.frame(Richness = wfv, Area = xx2)
  CI$x <- df$A
  G <- ggplot(DF2, aes(x = Area, y = Richness)) + geom_path() + theme_bw()
  if(confInt){
    G <- G + geom_ribbon(data = CI, aes(x = x, ymax = U, ymin = L))
  }
  print(G)
  return(DF2)
}
```


```{r fitPlantAndAnimals, echo = F, cache=T}

Richness <- read.csv("Animals_plants_richness.csv") |>
  dplyr::filter(!is.na(N)) |>
  dplyr::select(Ha, N)

mm_Richness <- sar_average(data = Richness, verb = FALSE, confInt = F, display = F)
```

Here we can see the table of selected models:

```{r table_Plants_Animals, echo = F}
Sum <-  summary(mm_Richness)$Model_table

Sum$CumWeight <- cumsum(Sum$Weight)
knitr::kable(Sum)
```

and the plot of the relationship

```{r plotAll}
plot(mm_Richness, pLeg = TRUE, mmSep = TRUE, allCurves = F)
```


## Plants only

Now we fit the same model but for plants only

```{r fitPlant, echo = F, cache=T}

Richness <- read.csv("Plants_richness.csv") |>
  dplyr::filter(!is.na(N)) |>
  dplyr::select(Ha, N)

mm_Richness_plants <- sar_average(data = Richness, verb = FALSE, confInt = F, display = F)
```

Here we can see the table of selected models:

```{r table_Plants, echo = F}
Sum <-  summary(mm_Richness_plants)$Model_table

Sum$CumWeight <- cumsum(Sum$Weight)
knitr::kable(Sum)
```

and the plot of the model

```{r plotPlants}
plot(mm_Richness_plants, pLeg = TRUE, mmSep = TRUE, allCurves = F)
```


## Animals only

Now we fit the same model but for Animals only


```{r fitAnimals, echo = F, cache=T}

Richness <- read.csv("Animals_richness.csv") |>
  dplyr::filter(!is.na(N)) |>
  dplyr::select(Ha, N)

mm_Richness_animals <- sar_average(data = Richness, verb = FALSE, confInt = F,  display = F)
```

Here we can see the table of selected models:

```{r table_animals, echo = F}
Sum <-  summary(mm_Richness_animals)$Model_table

Sum$CumWeight <- cumsum(Sum$Weight)
knitr::kable(Sum)
```

and the plot of the model

```{r plotanimals}
plot(mm_Richness_animals, pLeg = TRUE, mmSep = TRUE, allCurves = F)
```

# Buffer and intersection generation

## Tasks needed

* Existing nature (A), add a 5 ha buffer, add jakob's Forest dataset
* A = Existing nature
* B = Jakobs Forest or Fredskov
* Buffer ~5 Ha (225m) around B butcannot be part of A, and has to be a part of agriculture
* For every buffer patch we calculate the contiguity with B and A patches and calculate total area and proportion of forest within A and B and add Lavbund proportion


## Read A and plot

```{r loadPackages}
library(terra)
library(geodata)
library(ggplot2)
library(tidyterra)
library(BDRUtils)
```

## Existing Nature (A) and Forest (B) Datasets

We begin by loading the existing nature dataset (A) and the dataset representing deciduous forest areas from Jakob Assmann's work, which we'll refer to as B.

```{r ReadNatureandForest}
# Read in Nature
A <- terra::rast("o:/Nat_Sustain-proj/_user/derekCorcoran_au687614/biodiversitetsradet.github.io/Issue_20_BDR_KR/Class_A.tif")
# Read in forest
B <- terra::rast("o:/Nat_Ecoinformatics-tmp/JakobAssmann_au634851/treetype/treetype_bjer_dec.tif")
```

## Agricultural Land Dataset

Next, we acquire the dataset that represents agricultural land in Denmark. This information is crucial, as we want to ensure that the newly created forested areas are restricted to the existing agricultural land.

```{r Agriculture}
# Read in the dataset representing agricultural land in Denmark
Agriculture <- terra::rast("o:/Nat_BDR-data/Arealanalyse/2023/CLEAN/Rasterized/Rast_markblokkort_Croped.tif")
```


## Generating a Buffer around Deciduous Forest

We generate a 225-meter buffer around the deciduous forest areas, approximating a 5-hectare squared region.


```{r buffer225All, cache = TRUE}
# Create a copy of the deciduous forest dataset (B)
B_All <- B
# Set all zero values to 1 to simplify buffer creation
B_All[B_All == 0] <- 1
# Generate a 225-meter buffer around the deciduous forest
Buffer <- terra::buffer(B_All, 225)
# Create a temporary copy of the buffer to handle modifications
Temp <- Buffer
# Set all values to NA in the temporary buffer
values(Temp) <- NA
# Exclude existing natural areas from the buffer
Temp[Buffer] <- 1
A <- terra::extend(A, terra::ext(Temp))
Temp[!is.na(A)] <- NA
# Exclude forested areas from the buffer
Temp[!is.na(B_All)] <- NA
# Ensure only agricultural areas remain in the buffer
Temp[is.na(Agriculture)] <- NA
# Save the resulting buffer as a COG (Cloud-Optimized GeoTIFF)
BDRUtils::write_cog(Temp, "Buffer_all_225.tif")
```

This script creates a buffer around deciduous forest areas, ensuring it conforms to the specified conditions and constraints. The resulting buffer is then saved as a Cloud-Optimized GeoTIFF (COG) file named "Buffer_all_225.tif."

We can now visualize all this categories

```{r, echo = F, cache=TRUE}
Buffer <- terra::rast("SpatialData/Buffer_all_225.tif")

A <- terra::rast("o:/Nat_Sustain-proj/_user/derekCorcoran_au687614/biodiversitetsradet.github.io/Issue_20_BDR_KR/Class_A.tif")
A <- as.numeric(A)
A[A == 1] <- 2
A <- terra::extend(A, Buffer)
# Read in forest
B <- terra::rast("o:/Nat_Ecoinformatics-tmp/JakobAssmann_au634851/treetype/treetype_bjer_dec.tif")
B[!is.na(B)] <- 3

A[is.na(A)] <- 0
B[is.na(B)] <- 0
Buffer[is.na(Buffer)] <- 0


Final <- A+B+Buffer

Final[Final == 0] <- NA

df <- data.frame(id = c(1,2,3,5), levels = c("Potential Forest", "Non-forest existing Nature", "Forest", "Forest in existing Nature"))


levels(Final) <- df

BDRUtils::write_cog(Final, "AllCategories.tif")

ggplot() + geom_spatraster(data = Final, maxcell = 5e+06)
```


Now in order to calculate areas and adjacencies the raster will be transformed into polygons

```{r ToShapefiles, cache=TRUE, eval = T}
Buffer_all_225 <- terra::rast("AllCategories.tif")

Buffer_all_225_SF <- as.polygons(Buffer_all_225) |> aggregate(by = "levels") |> 
  disagg()
Buffer_all_225_SF$ID <- 1:nrow(Buffer_all_225_SF)
terra::writeVector(Buffer_all_225_SF, "Buffer_all_225_SF.shp", overwrite = T)
```

Now to actually calculate the values we will first unite resolve for A and B join them in the largest possible polygons

```{r uniteAB, cache = T}
B <- terra::rast("o:/Nat_Ecoinformatics-tmp/JakobAssmann_au634851/treetype/treetype_bjer_dec.tif")
B[!is.na(B)] <- 1
A <- terra::rast("o:/Nat_Sustain-proj/_user/derekCorcoran_au687614/biodiversitetsradet.github.io/Issue_20_BDR_KR/Class_A.tif")
A <- as.numeric(A)
A <- terra::extend(A, B)

A[is.na(A)] <- 0
B[is.na(B)] <- 0

AB <- A+B
AB[AB > 0] <- 1
AB[AB == 0] <- NA

names(AB) <- "level"

AB_SF <- as.polygons(AB) |> aggregate(by = "level") |> disagg()
AB_SF$ID <- 1:nrow(AB_SF)
terra::writeVector(AB_SF, "AB_SF.shp", overwrite = T)
```
Now we filter form the joint polygons only the ones with areas higher than 200 ha, 100 Ha, 50 Ha and 25 Ha

```{r filter200ha, cache=TRUE}
AB <- terra::vect("AB_SF.shp")
AB$Ha <- terra::expanse(AB, unit = "ha")
AB_25 <- AB[AB$Ha > 25,]

terra::writeVector(AB_25, "AB_25.shp", overwrite = T)

AB_50 <- AB_25[AB_25$Ha > 50,]
terra::writeVector(AB_50, "AB_50.shp", overwrite = T)


AB_100 <- AB_50[AB_50$Ha > 100,]
terra::writeVector(AB_100, "AB_100.shp", overwrite = T)


AB_200 <- AB_100[AB_100$Ha > 200,]
terra::writeVector(AB_200, "AB_200.shp", overwrite = T)

```


Now we go one by one and we generate the potential forest content

### 200 ha

```{r potentialForest200ha, cache=TRUE}
A <- terra::rast("o:/Nat_Sustain-proj/_user/derekCorcoran_au687614/biodiversitetsradet.github.io/Issue_20_BDR_KR/Class_A.tif")
# Read in forest
B <- terra::rast("o:/Nat_Ecoinformatics-tmp/JakobAssmann_au634851/treetype/treetype_bjer_dec.tif")

B[!is.na(B)] <- 1
Agriculture <- terra::rast("o:/Nat_BDR-data/Arealanalyse/2023/CLEAN/Rasterized/Rast_markblokkort_Croped.tif")

AB_200 <- terra::vect("AB_200.shp")

B_200 <- B |> 
  terra::crop(AB_200) |> 
  terra::mask(AB_200)

Buffer_200_ha <- terra::buffer(B_200, 225)
# Create a temporary copy of the buffer to handle modifications
Temp <- Buffer_200_ha
# Set all values to NA in the temporary buffer
values(Temp) <- NA
# Exclude existing natural areas from the buffer
Temp[Buffer_200_ha] <- 1
A <- terra::crop(A, terra::ext(Temp))
Temp[!is.na(A)] <- NA
# Exclude forested areas from the buffer
B <- terra::crop(B, terra::ext(Temp))
Temp[!is.na(B)] <- NA
# Ensure only agricultural areas remain in the buffer
Agriculture <- terra::crop(Agriculture, terra::ext(Temp))
Temp[is.na(Agriculture)] <- NA
names(Temp) <- "Patches"
# Save the resulting buffer as a COG (Cloud-Optimized GeoTIFF)
BDRUtils::write_cog(Temp, "potentialForest200ha.tif")
```

And now we transform this into polygons

```{r ToShapefiles200ha, cache=TRUE, eval = T}
potentialForest200ha <- terra::rast("potentialForest200ha.tif")

potentialForest200ha_SF <- as.polygons(potentialForest200ha) |> aggregate(by = "Patches") |> 
  disagg()
potentialForest200ha_SF$ID <- 1:nrow(potentialForest200ha_SF)

potentialForest200ha_SF$Ha <- terra::expanse(potentialForest200ha_SF, unit = "ha")

terra::writeVector(potentialForest200ha_SF, "potentialForest200ha_SF.shp", overwrite = T)
```

We now we process the potential forest to add total area considering AB (Total_Area), the area of the potential forest is included (Ha), we caclulate the area of forest considering the adjacent A and B areas, Forest_Area

```{r adjacency200, eval = T}
PotentialForest <- terra::vect("potentialForest200ha_SF.shp")

PotentialForest <-PotentialForest[PotentialForest$Ha > 5,]
PotentialForest$Total_Ha <- NA
PotentialForest$Ajdacent_AB_ID <- NA
PotentialForest$Forest_Ha <- NA

AB_200 <- terra::vect("AB_200.shp")

B <- terra::rast("o:/Nat_Ecoinformatics-tmp/JakobAssmann_au634851/treetype/treetype_bjer_dec.tif")

B[!is.na(B)] <- 1

for(i in 1:nrow(PotentialForest)){
  try({
    Temp <- terra::nearby(PotentialForest[i,], AB_200,  distance = 100, centroids = F, symmetrical = F) |> as.data.frame()
  PotentialForest$Total_Ha[i] <- PotentialForest$Ha[i] + AB_200[Temp$to_id,]$Ha
  PotentialForest$Ajdacent_AB_ID[i] <- paste(Temp$to_id, collapse = ", ")
  PotentialForest$Forest_Ha[i] <- terra::crop(B,  AB_200[Temp$to_id,]) |> terra::mask(AB_200[Temp$to_id,]) |> freq() |> dplyr::filter(value == 1) |> dplyr::pull(count) |> magrittr::multiply_by(100) |> magrittr::multiply_by(0.0001)
  })
  
}
terra::writeVector(PotentialForest, "PotentialForest200Ha_param.shp")
```



### 100 ha

```{r potentialForest100ha, cache=TRUE}
A <- terra::rast("o:/Nat_Sustain-proj/_user/derekCorcoran_au687614/biodiversitetsradet.github.io/Issue_20_BDR_KR/Class_A.tif")
# Read in forest
B <- terra::rast("o:/Nat_Ecoinformatics-tmp/JakobAssmann_au634851/treetype/treetype_bjer_dec.tif")

B[!is.na(B)] <- 1
Agriculture <- terra::rast("o:/Nat_BDR-data/Arealanalyse/2023/CLEAN/Rasterized/Rast_markblokkort_Croped.tif")

AB_100 <- terra::vect("AB_100.shp")

B_100 <- B |> 
  terra::crop(AB_100) |> 
  terra::mask(AB_100)

Buffer_100_ha <- terra::buffer(B_100, 225)
# Create a temporary copy of the buffer to handle modifications
Temp <- Buffer_100_ha
# Set all values to NA in the temporary buffer
values(Temp) <- NA
# Exclude existing natural areas from the buffer
Temp[Buffer_100_ha] <- 1
A <- terra::crop(A, terra::ext(Temp))
Temp[!is.na(A)] <- NA
# Exclude forested areas from the buffer
B <- terra::crop(B, terra::ext(Temp))
Temp[!is.na(B)] <- NA
# Ensure only agricultural areas remain in the buffer
Agriculture <- terra::crop(Agriculture, terra::ext(Temp))
Temp[is.na(Agriculture)] <- NA
names(Temp) <- "Patches"
# Save the resulting buffer as a COG (Cloud-Optimized GeoTIFF)
BDRUtils::write_cog(Temp, "potentialForest100ha.tif")
```

And now we transform this into polygons

```{r ToShapefiles100ha, cache=TRUE, eval = T}
potentialForest100ha <- terra::rast("potentialForest100ha.tif")

potentialForest100ha_SF <- as.polygons(potentialForest100ha) |> aggregate(by = "Patches") |> 
  disagg()
potentialForest100ha_SF$ID <- 1:nrow(potentialForest100ha_SF)
terra::writeVector(potentialForest100ha_SF, "potentialForest100ha_SF.shp", overwrite = T)
```