README.Rmd

---
output: github_document
---

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


[![Build Status](https://travis-ci.org/ices-tools-dev/sfdSAR.svg?branch=master)](https://travis-ci.org/ices-tools-dev/sfdSAR)

[<img align="right" alt="ICES Logo" width="17%" height="17%" src="http://ices.dk/_layouts/15/1033/images/icesimg/iceslogo.png">](http://ices.dk)

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  out.width = "100%"
)
```
# sfdSAR

The goal of sfdSAR is to make it easy to follow the procedure of calculating swept
area ratio of an area of seabed by a fishing gear.

## Installation

You can install the sfdSAR using

```r
install.packages("sfdSAR", repos = "https://ices-tools-prod.r-universe.dev")
```

## Usage

For a summary of the package:

```r
library(sfdSAR)
?sfdSAR
```

## References


ICES 2015. Report of the Working Group on Spatial Fisheries Data (WGSFD), 8-12
June 2015, ICES Headquarters, Copenhagen, Denmark. ICES CM 2015/SSGEPI:18. 150pp

ICES 2016. Interim Report of the Working Group on Spatial Fisheries Data (WGSFD),
17-20 May 2016, Brest, France. ICES CM 2016/SSGEPI:18. 244 pp

Eigaard OR, Bastardie F, Breen M, et al. (2016) Estimating seabed pressure from
demersal trawls, seines, and dredges based on gear design and dimensions. ICES
Journal of Marine Science, 73:27-43

Church N.J., Carter A.J., Tobin D., Edwards D., Eassom A., Cameron A., Johnson
G.E., Robson, L.M. & Webb K.E. (2016) JNCC Recommended Pressure Mapping
Methodology 1. Abrasion: Methods paper for creating a geo-data layer for the
pressure 'Physical Damage (Reversible Change) - Penetration and/or disturbance
of the substrate below the surface of the seabed, including abrasion'. JNCC
report No. 515, JNCC, Peterborough

Development
-----------

sfdSAR is developed openly on
[GitHub](https://github.com/ices-tools-dev/sfdSAR).

Feel free to open an
[issue](https://github.com/ices-tools-dev/sfdSAR/issues) there if you
encounter problems or have suggestions for future versions.


## Example

The functions in this package are intended for one purpose: to compute the swept
area ratio (SAR) and the subsurface SAR of a fishing gear, which can then be
summarised over years and gear groupings.

Swept Area Ratio (SAR) is computed using the algorithm described below. The main
steps in the data processing are

1. Determine the gear width of the VMS record according to:
    - Where average gear widths are supplied these are used.
    - For VMS records with missing gear widths but that have supplied average
      vessel characteristics (i.e. average overall vessel length or average KW
      engine power): use the model described in (Eigaard et al., 2016) to
      provide an estimate of gear width
    - For VMS records with missing gear widths and missing vessel
      characteristics use a fill-in value provided by ICES (2015) based on a
      review by the JNCC or on the BENTHIS survey (Eigaard et al. 2016).
1. Estimate swept area based on gear type, fishing hours (hours), fishing speed
    (speed) and gear width (width) for each record (ICES, 2016, p 69), note here
    speed is in knots and requires to be converted to km per hour:
    - Trawl          : hours x width x speed x 1.82
    - Danish seine   : hours / 2.591234 x (width<sup>2</sup>) / (4 &#960;)
    - Scottish seine : hours / 1.9125 x  (1.5 x width<sup>2</sup>) / (4 &#960;)
1. Accumulate across gears for each year to produce annual totals of SA by
    c-square and gear category, and finally average over years within gear
    category and c-square.
1. Calculate SAR values by scaling by the area of the c-squares

The code below shows how the sfdSAR functions can be used to calculate swept
area ratio (SAR)

In the following examples the `dplyr` package is used to simplify the data
processing and a made up vms toy vms dataset (`test_vms`) will be used
```{r packages, echo=FALSE, message=FALSE, warning=FALSE}
library(dplyr)
library(sfdSAR)
data(test_vms)
```
```
library(dplyr)
library(sfdSAR)
## load sample vms data
data(test_vms)
```

### 1. Determine gear widths

The calculation of gear with is done using the data in the benthis model parameters table:

```{r gear_info, results="asis"}
library(icesVMS)
gear_widths <- get_benthis_parameters()
kableExtra::kable(gear_widths)
```

This table comes from Eigaard et al. (2016), with additions from ICES (2015).
And contains, for each benthis gear group, the proportion of the gear
contact that also affects the subsurface, the estimated average gear width,
and the coeffients and covariates of the surface contact model which relates
the gear width, properties of the vessel (kw or overall length) to
bottom contact.

In order to use this data a lookup table is required linking Metier level 6
codes to the benthis gear groupings listed above.  The lookup table is given in
the `get_metier_lookup()` function from the `icesVMS` package and contains
other gear groupings used in ICES outputs and was initially developed by
ICES (2015).

```{r metier_lookup, results='asis'}
metier_lookup <- get_metier_lookup()
kableExtra::kable(head(metier_lookup))
```


Linking the gearwidths and contact model information is done with the following
two lines


```{r join}
# join widths and lookup
aux_lookup <-
  gear_widths %>%
  right_join(metier_lookup, by = c("benthisMet" = "benthisMetiers"))

# add aux data to vms
vms <-
  aux_lookup %>%
  right_join(test_vms, by = c("leMetLevel6" = "LE_MET_level6"))
```

and the gear width model is applied using the helper function `predict_gear_width`

```{r gearwidth}
# calculate the gear width model
vms$gearWidth_model <-
  predict_gear_width(vms$gearModel, vms$gearCoefficient, vms)
```

In general, if gearwdth is available, it is used.  If average overall vessel
length (oal) or average vessel power (kW) is available then the gear width
model is used.  Finally if none of these are avaiable an average gear width is
applied.  The following code implements this

```{r fillin}
# do the fillin for gear width:
# select provided average gear width, then modelled gear with, then benthis
# average if no kw or aol supplied
vms$gearWidth_filled <-
  with(vms,
    ifelse(!is.na(avg_gearWidth), avg_gearWidth / 1000,
      ifelse(!is.na(gearWidth_model), gearWidth_model / 1000,
        gearWidth)
    ))

```

### Predicting surface contact

finaly, surface contact is computed using the appropriate surface contact
model, given by the `contact_model` feild, defined as:

```{r contact_model}
sapply(unique(gear_widths$contactModel), function(x) body(get(x)))
```

The helper function `predict_surface_contact` computes the surface contact
(usage shown below).  The feild `subsurface_prop` which has come from the
`gear_width` dataset can be used to compute subsurface contact from the surface
contact.

```{r surface}
# calculate surface contact
vms$surface <-
  predict_surface_contact(vms$contactModel,
                          vms$fishingHours,
                          vms$gearWidth_filled,
                          vms$ICES_avg_fishing_speed)
# calculate subsurface contact
vms$subsurface <- vms$surface * as.numeric(vms$subsurfaceProp) * .01

```

### Summarising accross months etc.

Normally it is required to summarise the surface quantities, which can be done
like this

```{r summary}
# compute summaries of swept area over groups
sa <-
  vms %>%
    mutate(
      mw_fishinghours = kw_fishinghours / 1000
    ) %>%
    group_by(year, c_square, fishingCategoryFo) %>%
    summarise(
      mw_fishinghours = sum(mw_fishinghours, na.rm = TRUE),
      subsurface = sum(subsurface, na.rm = TRUE),
      surface = sum(surface, na.rm = TRUE)
    ) %>%
  ungroup %>%
  mutate(
    lat = csquare_lat(c_square),
    lon = csquare_lon(c_square)
  )

```

```{r sa}
sa
```

### Computing Swept Area Ratio (SAR)

In the code below, SAR is calculated for each year, then averaged over years,
resulting in a dataset of averarage SAR per c_square.  Note that, because
grouping is taking place over `c_square` the summation in the first
`group_by` section is equivalent to `sum(surface) / area`.  The second grouping
section computes averages for each `c_square` over all years in the dataset.

```{r sar}
# compute swept area ratio per year and c_square then average over years
sar <-
  sa %>%
    mutate(
      area = csquare_area(c_square)
    ) %>%
    group_by(c_square, year) %>%
      summarise(
        surface_sar = sum(surface / area, na.rm = TRUE),
        subsurface_sar = sum(subsurface / area, na.rm = TRUE)
      ) %>%
    ungroup() %>%
    group_by(c_square) %>%
    summarise(
      surface_sar = mean(surface_sar, na.rm = TRUE),
      subsurface_sar = mean(subsurface_sar, na.rm = TRUE)
    )
sar
```


### All in one

The steps described above are combined into one code block for convienience.
This code can be applied to a larger dataset to covering a range of years,
fishing gears and c_squares.


```{r ref.label=c("join", "gearwidth", "fillin", "surface", "summary", "sar")}

```

### Notes

This document was created using:
```{r, eval=FALSE}
rmarkdown::render("README.Rmd")
```