gfwrPackageStraightforwardApproach.Rmd

---
title: "gfwr package tool"
author: "Kilian GRIHAULT BARREIRO"
date: "2024-01-14"
output:
  html_document: default
  pdf_document: default
editor_options:
  markdown:
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE)
```

# How to use : Global Fishing Watch R Package (`gfwr`)

The `gfwr` package provides convenient functions to pull GFW data
directly into R into usable formats.\
It contains three main functions, including `get_vessel_info()`,
`get_event()` and `get_raster()`. The two first being devoted to
retrieving information and features on one ore several specific vessels.
The last is of particular interest to us because it allows us to gather
information from global fishing watch raster on the fishing effort
(further details in the function appropriate section).

The time spent fishing is computed using Automatic Identification System
(AIS) data, which is transmitted by most industrial fishing vessels. The
AIS data provides information on the location, speed, and direction of
the vessel, which can be used to identify when the vessel is actively
fishing.

#### AIS Caveats and limitations

The AIS coverage of vessels has several limitations such as :<br> 1. The
number of vessels that are captured (AIS provides approximately 70'000
of the 2.8 million identified fishing vessels). <br> 2. The size of the
vessels (52-85% for vessels larger than 24 meters against 1% for vessels
under 12 meters). <br> *Good to know: IMO mandates AIS for most vessels
larger than 36 meters.* <br> 3. AIS interference with each other in
areas of high vessel density. <br> 4. Some terrestrial satellites only
receive messages near shore. <br>

## Installation

```{r, eval=FALSE}
remotes::install_github("GlobalFishingWatch/gfwr")
```

```{r, message=FALSE}
library(gfwr)
library(ggplot2)
library(tidyverse)
library(sf)
```

## API

To access GFW APIs, you need to :<br> 1. register for a GFW account here

:   \<[gfw account creation link](https://tinyurl.com/2pr5nxxk)\> <br> 2. Request
    API key here : <https://globalfishingwatch.org/our-apis/tokens>

Once you have your token, add it to your .Renviron file (by executing
the chunk below), by writing (GFW_TOKEN = "YOUR_TOKEN") in the file.
<br> *(You could be asked to restart R for changes to take effect.)*

```{r, eval=FALSE}
usethis::edit_r_environ()
```

We save the key in an object that will be used in gfwr functions.

```{r}
key <- gfwr::gfw_auth()
```

## Fishing effort visualization

A region_id is necessary to use the `get_raster` function.

```{r, results='hide'}
region_id <- get_region_id(
  region_name = "Australia",
  region_source = "eez",
  key = key
)$id
```

The `get_raster` function gets a raster of fishing effort from the API
and converts the response to a data frame which contains occurrences for
each vessel and for each grid cell (data is binned into grid cells of
different resolution), the `Vessel IDs`, `Flag`, `Geartype` and
`Apparent fishing Hours` which are basically the amount of fishing hours
of each vessel per grid cell (`geometry`).

Data can be provided through :<br> - `daily`, `monthly` and `yearly`
temporal resolutions. <br> - `low` (0.1 deg) and `high` (0.01 deg)
spatial resolutions. <br> - `vessel_id`, `flag`, `gearType`,
`flagAndGearType`.

```{r, eval=FALSE, message=FALSE}
data <- get_raster(
  spatial_resolution = "low",
  temporal_resolution = "monthly",
  group_by = "flagAndGearType",
  date_range = "2022-01-01,2023-01-01",
  region = region_id,
  region_source = "eez",
  key = key
)
```

*(You can remove the option* `message = FALSE` *to see the columns
types.)*

#### `get_raster` caveats and limitations.

Date range is limited to 1-year. Nevertheless, with some modifications,
we can get round these problems through `get_gfwFishingEffort`.

#### `get_gfwFishingEffort` function.

The `get_gfwFishingEffort` function recover the data of Global Fishing
Watch and returns it as a sf object. <br> We have the same parameters
than the `get_raster` function, plus `n_crs` which is the crs for the
sf_modification <br> Different possible values can be combined and are :
<br> - `Time Range`, `Flag`, `Geartype`. <br> *(A combination can be :
c('Time Range','Geartype'), if you want to get the sum of fishing hours
per date and geartype, for example you want to display the drifting
longline fishing in a specific year)* <br> <br> **Notes :** <br> 1. For
the moment we are limited to the EEZs of each region, but we can
potentially restrict the working area to specific MPAs (further details
in the gfwr package). <br> 2. Days indicated in the\_\_ `start_date`
**and** `end_date` \_\_variables are included in the data recovery.

```{r}
get_gfwData <- function(region, start_date, end_date, temp_res,
                        spat_res = "low",
                        key = gfwr::gfw_auth(),
                        cCRS = "EPSG:4326",
                        compress = FALSE) {
  region_id <- gfwr::get_region_id(
    region_name = region,
    region_source = "eez",
    key = key
  )$id[1]

  # Convert dates into Date objects
  start_date <- as.Date(start_date, format = "%Y-%m-%d")
  end_date <- as.Date(end_date, format = "%Y-%m-%d")

  # Function to obtain data for a specific date range
  get_data_for_range <- function(start_date, end_date) {
    date_range <- paste(start_date, end_date, sep = ",")

    data <- gfwr::get_raster(
      spatial_resolution = spat_res,
      temporal_resolution = temp_res,
      group_by = "flagAndGearType",
      date_range = date_range,
      region = region_id,
      region_source = "eez",
      key = key
    )

    return(data)
  }

  # Check whether the date range is less than or equal to 366 days
  if (as.numeric(difftime(end_date, start_date, units = "days")) <= 366) {
    # If yes, obtain data for the entire date range
    data_df <- get_data_for_range(start_date, end_date)
  } else {
    # If not, divide the date range into 366-day chunks and obtain the data for each chunk.
    date_chunks <- seq(start_date, end_date, by = "366 days")
    data_df <- purrr::map_dfr(
      date_chunks,
      ~ get_data_for_range(.x, min(.x + 365, end_date))
    )
  }

  if (isTRUE(compress)) {
    # GFW data will always be "EPSG:4326". No need to have CRS as an option here

    data_df <- data_df %>%
      dplyr::select("Lon", "Lat", "Apparent Fishing Hours") %>%
      dplyr::group_by(Lon, Lat) %>%
      dplyr::summarise("Apparent Fishing Hours" = sum(`Apparent Fishing Hours`,
        na.rm = TRUE
      )) %>%
      dplyr::ungroup()

    data_sf <- data_df %>%
      terra::rast(type = "xyz", crs = "EPSG:4326") %>% # Convert to polygons for easier use
      terra::as.polygons(
        trunc = FALSE,
        dissolve = FALSE,
        na.rm = TRUE,
        round = FALSE
      ) %>%
      sf::st_as_sf()

    if (dim(data_df)[1] != dim(data_sf)[1]) {
      stop("Data dimensions of df and sf data do not match after conversion")
    }
  } else if (isFALSE(compress)) {
    # Combine data frames in the list into one data frame
    data_df <- bind_rows(data_df)

    # Separate the "Time Range" column based on the specified temp_res
    if (temp_res == "yearly") {
      data_sf <- data_df %>%
        dplyr::mutate(Year = `Time Range`) %>%
        sf::st_as_sf(
          coords = c("Lon", "Lat"),
          crs = "EPSG:4326"
        )
    } else {
      # Sinon, séparer la colonne "Time Range" selon le temp_res spécifié
      if (temp_res == "monthly") {
        data_sf <- data_df %>%
          tidyr::separate("Time Range",
            into = c("Year", "Month"),
            sep = "-",
            remove = FALSE
          ) %>%
          sf::st_as_sf(
            coords = c("Lon", "Lat"),
            crs = "EPSG:4326"
          )
      } else if (temp_res == "daily") {
        data_sf <- data_df %>%
          tidyr::separate("Time Range",
            into = c("Year", "Month", "Day"),
            sep = "-",
            remove = FALSE
          ) %>%
          sf::st_as_sf(
            coords = c("Lon", "Lat"),
            crs = "EPSG:4326"
          )
      }
    }
  }

  # But you may wish to return the data in a different CRS. For this you need transform
  if (isFALSE(cCRS == "EPSG:4326")) {
    data_sf <- data_sf %>%
      sf::st_transform(crs = cCRS)
  }

  return(data_sf)
}
```

To get the raw values given by GFW we do not need to specify any
variable of gathering in `group_by_vars` parameter.

```{r, message=FALSE}
data_sf_combined <- get_gfwData(
  region = "Australie",
  start_date = "2021-01-01",
  end_date = "2023-12-31",
  temp_res = "yearly",
  spat_res = "low",
  key = key,
  compress = FALSE
)

data_sf_combined <- get_gfwData("Australie",
  "2019-01-01",
  "2023-12-31",
  "yearly",
  compress = FALSE
)
```

## Visualization

To display the data, we load : <br> - The coastline from `rnaturalearth`
package and modify it to get an sf object, and we constrain it to the
boundaries of the given data. <br> - EEZ Polygons from \_Flanders Marine
Institute (2023). Maritime Boundaries Geodatabase: Maritime Boundaries
and Exclusive Economic Zones (200NM), version 12. Available online at
<https://www.marineregions.org/>. <https://doi.org/10.14284/632_>

```{r, message = FALSE, results='hide'}
# Check and modify if necessary the spatial reference of data_sf_combined
data_sf_combined <- sf::st_set_crs(
  data_sf_combined,
  sf::st_crs(rnaturalearth::ne_coastline(scale = "large"))
)

coast_clipped <- rnaturalearth::ne_coastline(scale = "large") %>%
  sf::st_as_sf() %>%
  sf::st_intersection(sf::st_as_sfc(sf::st_bbox(data_sf_combined)))

# Load EEZ polygons
## Downloaded data
# eezs <- sf::st_read(dsn = '~/2024_Research_A_UQ/World_EEZ_v12_20231025/') %>%
#   sf::st_transform(crs = sf::st_crs(data_sf_combined)) %>%
#   sf::st_make_valid() %>%
#   sf::st_intersection(sf::st_as_sfc(sf::st_bbox(data_sf_combined)))

eezs <- mregions::mr_shp(key = "Australia:eez") %>%
  dplyr::filter(geoname == "Australian Exclusive Economic Zone") %>%
  sf::st_transform(crs = sf::st_crs(data_sf_combined)) %>%
  sf::st_make_valid() %>%
  sf::st_intersection(sf::st_as_sfc(sf::st_bbox(data_sf_combined)))
```

#### Here we display the Fishing Effort in Australia from 2019 to 2023.

##### Raw Fishing Effort

```{r, echo=FALSE}
main_plot <- ggplot(data_sf_combined) +
  geom_sf(aes(color = log10(`Apparent Fishing Hours`))) +
  geom_sf(data = coast_clipped, color = "black", fill = NA) + # Add coastline
  geom_sf(data = eezs, fill = NA, color = "red") + # Add the EEZ with hatching
  scale_color_viridis_c(guide = "legend") +
  theme_minimal() +
  labs(
    title = "2022 Vessel Activity Map",
    subtitle = "Fishing Hours recorded by GFW in Australia",
    color = "Fishing Hours (log10)"
  ) +
  theme(
    legend.position = "bottom",
    legend.text = element_text(size = 8),
    legend.title = element_text(size = 10)
  ) +
  guides(color = guide_colorbar(
    title.position = "top",
    title.vjust = 0.5,
    title.hjust = -0.5,
    label.theme = element_text(size = 8),
    barwidth = 5,
    barheight = 0.5
  ))

# The display and writing in this section is for information purposes only, to
# understand how the information on the grid is translated.
overlay_plot <- ggplot(data_sf_combined) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf), fill = "white") +
  geom_sf(aes(color = log10(`Apparent Fishing Hours`))) +
  geom_sf(data = coast_clipped, color = "black", fill = NA) + # Add coastline
  geom_sf(data = eezs, fill = NA, color = "red") +
  scale_color_viridis_c(guide = "legend") +
  labs(
    title = "Vessel Activity Map in Australia between 2019 and 2023",
    subtitle = "Fishing Hours data recorded by GFW",
    color = "Fishing Hours \n (log10)"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    title = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.border = element_rect(color = "black", fill = NA, size = 1)
  ) +
  coord_sf(xlim = c(152, 155), ylim = c(-27, -29))

main_plot +
  annotation_custom(
    ggplotGrob(overlay_plot),
    xmin = 130,
    xmax = 170,
    ymin = -20,
    ymax = -36
  )
```

##### By years

```{r, echo=FALSE}
ggplot(data_sf_combined) +
  geom_sf(aes(color = as.factor(Year))) +
  geom_sf(data = coast_clipped, color = "black", fill = NA) + # Add coastline
  geom_sf(data = eezs, color = "red", fill = NA) + # Ajouter la EEZ avec hachures
  theme_minimal() +
  scale_color_viridis_d(guide = "legend") +
  labs(
    title = "Vessel Activity Map in Australia between 2019 and 2024",
    subtitle = "Fishing Hours data recorded by GFW",
    color = "Years"
  ) +
  theme(
    legend.position = "bottom",
    legend.text = element_text(size = 8),
    legend.title = element_text(size = 10)
  )
```

##### Year-on-year comparison

We may need to compare different timeframes, such as seasons, to see if
there are any patterns. <br> **Note :** As more vessels have adopted AIS
(mainly in economically developed countries) since the deployment of
these technologies, the rise in activities must be seen in the context
of this increase and not necessarily of more intense fishing activity.

```{r, message=FALSE}
# We need to change the temporal range according to our need group by it to
# display the total fishing hours. <br>
data_sf_combined <- get_gfwData(
  "Australie",
  "2019-01-01",
  "2023-12-31",
  "monthly"
) %>%
  dplyr::group_by(Year, Month) %>%
  dplyr::summarize(Total_Fishing_Hours = sum(`Apparent Fishing Hours`))
```

```{r, echo=FALSE}
ggplot(
  data_sf_combined,
  aes(x = Month, y = Total_Fishing_Hours, color = Year, group = Year)
) +
  geom_line() +
  geom_point() +
  labs(
    title = "Total Fishing Hours per month (2014-2023)",
    x = "Month", y = "Total Fishing Hours"
  ) +
  theme_minimal()
```

##### Fishing gear type

Here we display the Vessel activity in 'Micronesia' in 2022 according to
the fishing gear type.

```{r, message=FALSE}
data_sf_combined <- get_gfwData(
  "Micronesia",
  "2022-01-01",
  "2022-12-31",
  "monthly",
  "low",
  key
)
```

```{r, echo=FALSE, message=FALSE, results='hide'}
# Check and modify if necessary the spatial reference of data_sf_combined
data_sf_combined <- sf::st_set_crs(
  data_sf_combined,
  sf::st_crs(rnaturalearth::ne_coastline(scale = "large"))
)

coast_clipped <- rnaturalearth::ne_coastline(scale = "large") %>%
  sf::st_as_sf() %>%
  sf::st_intersection(sf::st_as_sfc(sf::st_bbox(data_sf_combined)))

# Load EEZ polygons
## Downloaded data
# eezs <- sf::st_read(dsn = '~/2024_Research_A_UQ/World_EEZ_v12_20231025/') %>%
#   sf::st_transform(crs = sf::st_crs(data_sf_combined)) %>%
#   sf::st_make_valid() %>%
#   sf::st_intersection(sf::st_as_sfc(sf::st_bbox(data_sf_combined)))

eezs <- mregions::mr_shp(key = "Micronesia:eez") %>%
  dplyr::filter(geoname == "Micronesian Exclusive Economic Zone") %>%
  sf::st_transform(crs = sf::st_crs(data_sf_combined)) %>%
  sf::st_make_valid() %>%
  sf::st_intersection(sf::st_as_sfc(sf::st_bbox(data_sf_combined)))
```

```{r, echo=FALSE}
# Create the map
ggplot(data_sf_combined) +
  geom_sf(aes(color = `Geartype`)) +
  geom_sf(data = coast_clipped, color = "black", fill = NA) + # Add coastline
  geom_sf(data = eezs, color = "red", fill = NA) + # Add EEZ
  theme_minimal() +
  labs(
    title = "2022 Vessel Activity Map",
    subtitle = "recorded by GFW in Micronesia",
    color = "Gear types"
  ) +
  theme(legend.position = "right")
```

##### Flags

Here we display the Vessel activity in Papua New Guinea according to
Vessels flags.

```{r, echo=FALSE, message=FALSE, results='hide'}
data_sf_combined <- get_gfwData(
  "Papua New Guinea",
  "2022-01-01",
  "2022-12-31",
  "monthly",
  "low",
  key
) %>%
  sf::st_set_crs(sf::st_crs(rnaturalearth::ne_coastline(scale = "large")))

coast_clipped <- rnaturalearth::ne_coastline(scale = "large") %>%
  sf::st_as_sf() %>%
  sf::st_intersection(sf::st_as_sfc(sf::st_bbox(data_sf_combined)))

# Load EEZ polygons
# eezs <- sf::st_read(dsn = '~/2024_Research_A_UQ/World_EEZ_v12_20231025/') %>%
#   sf::st_transform(crs = sf::st_crs(data_sf_combined)) %>%
#   sf::st_make_valid() %>%
#   sf::st_intersection(sf::st_as_sfc(sf::st_bbox(data_sf_combined)))

eezs <- mregions::mr_shp(key = "Papua New Guinea:eez") %>%
  dplyr::filter(geoname == "Papua New Guinean Exclusive Economic Zone") %>%
  sf::st_transform(crs = sf::st_crs(data_sf_combined)) %>%
  sf::st_make_valid() %>%
  sf::st_intersection(sf::st_as_sfc(sf::st_bbox(data_sf_combined)))
```

```{r, echo=FALSE}
# Create the map
ggplot(data_sf_combined) +
  geom_sf(aes(color = `Flag`)) +
  geom_sf(data = coast_clipped, color = "black", fill = NA) + # Add coastline
  geom_sf(data = eezs, color = "red", fill = NA) + # Add EEZ
  scale_size_continuous(
    range = c(1, 10),
    guide = "legend",
    name = "Flag"
  ) +
  theme_minimal() +
  labs(
    title = "2022 Vessel Activity Map",
    subtitle = "recorded by GFW in Papua New Guinea",
    color = "Flag"
  ) +
  theme(legend.position = "right")
```

##### Supplementary materials.

The fishing detection model was trained on AIS data from 503 vessels and
identified fishing activity with over 90% accuracy, which means that it
can identify a fishing and non-fishing activity with high accuracy. More
details on AIS operation and limitations here : <br>
<https://globalfishingwatch.org/dataset-and-code-fishing-effort/> <br>

###### Hierarchy of vessels gear types :

<https://globalfishingwatch.org/datasets-and-code-vessel-identity/>\*