report-QC_format2022.Rmd

---
title: "ICES VMS datacall quality check report"
date: "`r lubridate::now()`"
output: 
  html_document: 
    fig_width: 9
    fig_height: 4
    toc: yes
    toc_float: true
---

<!-- QCTEMPLATE: header -->


```{r setup, echo = FALSE}
knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE,
                      comment = NA)
```

```{r setup2}
library(tidyverse)
library(knitr)
library(icesVocab)
library(DT)

options(digits = 4)

#Specify data path
veDataPath <- "~/prj2/vms/ices_vms_data_call/delivery/iceland_annex1_2009_2022_2023-08-14.csv"

leDataPath <- "~/prj2/vms/ices_vms_data_call/delivery/iceland_annex2_2009_2022_2023-08-14.csv"



# Functions --------------------------------------------------------------------
grade <-function (x, dx) {
  if (dx > 1) 
    warning("Not tested for grids larger than one")
  brks <- seq(floor(min(x)), ceiling(max(x)), dx)
  ints <- findInterval(x, brks, all.inside = TRUE)
  x <- (brks[ints] + brks[ints + 1])/2
  return(x)
}

`an` <-
function(x){return(as.numeric(x))}

csq2pos <- function(d, dx = 0.05) {
  bind_cols(d,
            CSquare2LonLat(d$csq, degrees = dx) %>% 
              select(lon = SI_LONG,
                     lat = SI_LATI)) %>% 
    as_tibble()
}


# borrowed from the geo-package
d2ir <- function (lon, lat = NULL, useI = FALSE) {
  
  # if lon is a dataframe
  if (is.null(lat)) {
    lon <- lat$lon
    lat <- lat$lat
  }
  lat <- lat + 1e-06
  lon <- lon + 1e-06
  outside <- lat < 36 | lat >= 85.5 | lon <= -44 | lon > 68.5
  if (any(outside)) 
    warning("Positions outside of ICES statistical area")
  lat <- floor(lat * 2) - 71
  lat <- ifelse(lat < 10, paste("0", lat, sep = ""), lat)
  if (useI) 
    lettersUsed <- LETTERS[1:12]
  else lettersUsed <- LETTERS[c(1:8, 10:13)]
  lon1 <- lettersUsed[(lon + 60)%/%10]
  lon2 <- ifelse(lon1 == "A", floor(lon%%4), floor(lon%%10))
  ir <- paste(lat, lon1, lon2, sep = "")
  ir[outside] <- NA
  ir
}
ir2d <- function (ir, useI = FALSE) {
  lat <- substring(ir, 1, 2)
  lat <- as.numeric(lat)
  lat <- (lat + 71)/2 + 0.25
  lon1 <- substring(ir, 3, 3)
  lon1 <- toupper(lon1)
  lon1 <- match(lon1, LETTERS)
  if (!useI) 
    lon1 <- ifelse(lon1 > 8, lon1 - 1, lon1)
  lon1 <- lon1 - 2
  lon2 <- substring(ir, 4)
  lon2 <- as.numeric(lon2)
  lon <- ifelse(lon1 < 0, -44 + lon2 + 0.5, -40 + 10 * lon1 + 
                  lon2 + 0.5)
  data.frame(lat = lat, lon = lon)
}

CSquare2LonLat <- function(csqr,degrees){

  ra          <- 1e-6 #Artificial number to add for rounding of 5'ves (round(0.5) = 0, but in Excel (where onversion comes from, it is 1)onversion comes from, it is 1)  
  chars       <- an(nchar(csqr))
  tensqd      <- an(substr(csqr,1,1))+ra;      gqlat     <- (round(abs(tensqd-4)*2/10)*10/5)-1
  tenslatd    <- an(substr(csqr,2,2))+ra;      gqlon     <- (2*round(tensqd/10)-1)*-1
  tenslond    <- an(substr(csqr,3,4))+ra
  unitsqd     <- an(substr(csqr,6,6))+ra;      iqulat    <- round(unitsqd*2/10)
  unitslatd   <- an(substr(csqr,7,7))+ra;      iqulon    <- (round((unitsqd-1)/2,1) - floor((unitsqd-1)/2))*2
  unitslond   <- an(substr(csqr,8,8))+ra
  tenthsqd    <- an(substr(csqr,10,10))+ra;    iqtlat    <- round(tenthsqd*2/10)
  tenthslatd  <- an(substr(csqr,11,11))+ra;    iqtlon    <- (round((tenthsqd-1)/2,1) - floor((tenthsqd-1)/2))*2
  tenthslond  <- an(substr(csqr,12,12))+ra
  hundqd      <- an(substr(csqr,14,14))+ra;    iqhlat    <- round(hundqd*2/10)
  hundlatd    <- an(substr(csqr,15,15))+ra;    iqhlon    <- (round((hundqd-1)/2,1) - floor((hundqd-1)/2))*2
  hundlond    <- an(substr(csqr,16,16))+ra
  reso        <- 10^(1-floor((chars-4)/4))-((round((chars-4)/4,1)-floor((chars-4)/4))*10^(1-floor((chars-4)/4)))

  if(degrees < reso[1]) stop("Returning degrees is smaller than format of C-square")
  if(degrees == 10){   
    lat <- ((tenslatd*10)+5)*gqlat-ra                              
    lon <- ((tenslond*10)+5)*gqlon-ra              
  }
  if(degrees == 5){    
    lat <- ((tenslatd*10)+(iqulat*5)+2.5)*gqlat-ra
    lon <- ((tenslond*10)+(iqulon*5)+2.5)*gqlon-ra 
  }
  if(degrees == 1){    
    lat <- ((tenslatd*10)+ unitslatd+0.5)*gqlat-ra
    lon <- ((tenslond*10)+ unitslond+0.5)*gqlon-ra 
  }
  if(degrees == 0.5){  
    lat <- ((tenslatd*10)+ unitslatd + (iqtlat*0.5)+0.25)*gqlat-ra
    lon <- ((tenslond*10)+ unitslond + (iqtlon*0.5)+0.25)*gqlon-ra
  }         
  if(degrees == 0.1){  
    lat <- ((tenslatd*10)+ unitslatd + (tenthslatd*0.1)+0.05)*gqlat-ra
    lon <- ((tenslond*10)+ unitslond + (tenthslond*0.1)+0.05)*gqlon-ra
  }
  if(degrees == 0.05){ 
    lat <- ((tenslatd*10)+ unitslatd + (tenthslatd*0.1)+(iqhlat*0.05)+0.025)*gqlat-ra
    lon <- ((tenslond*10)+ unitslond + (tenthslond*0.1)+(iqhlon*0.05)+0.025)*gqlon-ra
  }
  if(degrees == 0.01){ 
    lat <- ((tenslatd*10)+ unitslatd + (tenthslatd*0.1)+(hundlatd*0.01)+0.005)*gqlat-ra
    lon <- ((tenslond*10)+ unitslond + (tenthslond*0.1)+(hundlond*0.01)+0.005)*gqlon-ra
}
return(data.frame(SI_LATI=lat,SI_LONG=lon))}

read_le <- function(file) {
  read.csv(file,
           header = FALSE,
           na.strings = "NULL") %>% 
    tidyr::as_tibble() %>% 
    dplyr::rename(type = 1,
                  country = 2,
                  year = 3,
                  month = 4,
                  n_vessel = 5,
                  vids = 6,
                  isq = 7,
                  gear = 8,
                  target = 9,
                  metier6 = 10,
                  lclass = 11,
                  vms = 12,
                  effort = 13,
                  kwd = 14,
                  catch = 15,
                  value = 16) %>% 
    dplyr::mutate(value = as.numeric(value),
                  lon = ir2d(isq)$lon,
                  lat = ir2d(isq)$lat)
}
read_ve <- function(file) {
  
 read.csv(file,
           header = FALSE,
           na.strings = "NULL") %>% 
    tidyr::as_tibble() %>% 
    dplyr::rename(type = 1,
                  country = 2,
                  year = 3,
                  month = 4,
                  n_vessel = 5,
                  vids = 6,
                  csq = 7,
                  gear = 8,
                  target = 9,
                  metier6 = 10,
                  lclass = 11,
                  speed = 12,
                  effort = 13,             # units in hours
                  int = 14,                # average interval
                  length = 15,             # average vessel length
                  kw = 16,                 # average kw
                  kwh = 17,                # kw x effort
                  catch = 18,
                  value = 19,
                  gearwidth = 20) %>%
#                  spread = 21) %>% 
    dplyr::mutate(value = as.numeric(value))
}

# support tables ---------------------------------------------------------------
lclass <- 
  tribble(~lclass, ~length_class,
          "A", "<8",
          "B", "08-10",
          "C", "10-12",
          "D", "12-15",
          "E", ">= 15")

# Get data ---------------------------------------------------------------------
ve <- read_ve(veDataPath) 
ve <- bind_cols(ve,
            CSquare2LonLat(ve$csq, degrees = 0.05) %>% 
              select(lon = SI_LONG,
                     lat = SI_LATI)) %>% 
    as_tibble()

ve <- ve %>% 
  filter(!is.na(lon)) %>% 
  mutate(isq = d2ir(lat, lon))

le <- read_le(leDataPath)

#file <- veDataPath 
# some bounds ------------------------------------------------------------------
country <- unique(ve$country)
spatBound <- list(xrange = range(ve$lon, na.rm = TRUE),
                  yrange = range(ve$lat, na.rm = TRUE))
spatCore  <- list(xrange = quantile(ve$lon, c(0.025, 0.975)),
                  yrange = quantile(ve$lat, c(0.025, 0.975)))
tempBound <- range(ve$year, na.rm = TRUE)
# some maps --------------------------------------------------------------------
m <- map_data("world", xlim = spatBound$xrange, ylim = spatBound$yrange)
map0 <- 
  ggplot() +
  theme_void() +
  geom_polygon(data = m, aes(x = long, y = lat, group = group),
               fill = "lightgray")
map1 <- 
  map0 +
  coord_quickmap(xlim = spatBound$xrange, ylim = spatBound$yrange)
map2 <- 
  map0 +
  coord_quickmap(xlim = spatCore$xrange, ylim = spatCore$yrange)
```

<!--------------------------
Checks for vms/ais data
------------------------ -->

```{r logbook-checks}
spatBoundLog <- list(xrange = range(le$lon, na.rm = TRUE),
                     yrange = range(le$lat, na.rm = TRUE))
tempBoundLog <- range(le$year, na.rm = TRUE)
```

# VMS

## Overall checks (fishing hours) {.tabset}


### By month


```{r by_month}
ve %>% 
  group_by(year, month) %>% 
  summarise(FH = sum(effort),
            .groups = "drop") %>% 
  spread(year, FH) %>% 
  kable(caption = "Effort(fishing hours)")

ve %>% 
  group_by(year) %>% 
  summarise(FH = sum(effort)) %>% 
  ggplot(aes(year, FH)) +
  geom_line(lwd = 1) +
  labs(x = NULL, y = "Fishing hours",
       title = "Effort (fishing hours) by year")

 
```


```{r pings-by-month}
ve %>% 
  group_by(year, month) %>% 
  summarise(FH=sum(effort)) %>%
  ggplot(aes(year, FH)) +
  geom_point() +
  geom_line() +
  facet_wrap(~ month) +
  labs(x = NULL, y = NULL, title = "Fishing hours by year & month")

#  complete(year, month, fill = list(n = 0)) %>% 

```

### By length class

```{r by_lclass}

ve %>% 
  group_by(year, lclass) %>% 
  summarise(FH = sum(effort),
            .groups = "drop") %>% 
  spread(year, FH) %>% 
  kable(caption = "Effort(fishing hours) by length class and year")

ve %>% 
  group_by(year, lclass) %>% 
  summarise(FH=sum(effort)) %>%
  left_join(lclass) %>% 
  unite(lclass, length_class, col = "lclass", sep = ": ") %>% 
  ggplot(aes(year, FH)) +
  geom_point() +
  geom_line() +
  facet_wrap(~ lclass, scales = "free_y") +
  labs(x = NULL, y = NULL, title = "Fishing hours by year & length class")

#complete(year, lclass, fill = list(n = 0)) %>% 
```

### By Gear

```{r by_dcf4}
ve %>% 
  group_by(year) %>% 
  summarise(gear = n_distinct(gear)) %>% 
  spread(year, gear) %>% 
  kable(caption = "Number of unique vms Gear by year")
ve %>% 
  count(year, gear) %>% 
  spread(year, n) %>% 
  kable(caption = "Number of vms records by Gear")
ve %>% 
  count(year, gear) %>% 
  ggplot(aes(year, n)) +
  geom_point() +
  geom_line() +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = NULL, title = "Number of vms records by Gear") 
```


### By Target

```{r by_dcf5}
ve %>% 
  group_by(year) %>% 
  summarise(target = n_distinct(target)) %>% 
  spread(year, target) %>% 
  kable(caption = "Number of unique vms Target by year")
ve %>% 
  count(year, target) %>% 
  spread(year, n) %>% 
  kable(caption = "Number of vms records by Target")

ve %>% 
  count(year, target) %>% 
  ggplot(aes(year, n)) +
  geom_point() +
  geom_line() +
  facet_wrap(~ target, scales = "free_y") +
  labs(x = NULL, y = NULL, title = "Number of vms records by Target") 
```


### By Metier6

```{r by_dcf6, fig.height = 6}
ve %>% 
  group_by(year) %>% 
  summarise(metier6 = n_distinct(metier6)) %>% 
  spread(year, metier6) %>% 
  kable(caption = "Number of unique vms Metier6 by year")
ve %>% 
  count(metier6, year) %>% 
  spread(year, n) %>% 
  kable(caption = "Number of vms records by Metier6 by year")
# add here labels on the metier6 level
ve %>% 
  count(year, gear, target, metier6) %>% 
  ggplot(aes(year, n, colour = target, group = metier6)) +
  geom_point() +
  geom_line() +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = NULL,
       title = "Number of vms records by Metier6",
       subtitle = "needs further work")
```


## Distribution {.tabset}

### Fishing speed

```{r table-average-fishing-speed, results='asis'}
kable(do.call(rbind, tapply(ve$speed, ve$year, summary)), booktabs = TRUE)
```

```{r average-fishing-speed, fig.height = 6}
ve %>% 
  mutate(speed = grade(speed, 0.2)) %>% 
  count(year, gear, speed) %>% 
  group_by(year, gear) %>% 
  mutate(p = n / sum(n)) %>% 
  ggplot(aes(speed, p, colour = gear)) +
  geom_line(lwd = 1) +
  labs(x = NULL,
       y = "Proportion",
       title = "Average fishing speed by Gear") +
  facet_wrap( ~ year, ncol = 2) +
  scale_colour_viridis_d()
```

### Fishing hours

```{r table-fishing-hours, results='asis'}
kable(do.call(rbind, tapply(ve$effort, ve$year, summary)), booktabs = TRUE,
      caption = "Summary statistics of average fishing hours")
```


```{r fishing-hours, fig.height = 6}
#d <- 
#  ve %>% 
#  group_by(year) %>% 
#  summarise(effort = max(effort))
ggplot(ve, aes(x = effort)) +
  geom_histogram() +
  scale_x_log10() +
  #geom_vline(data = d, aes(xintercept = effort)) +
  labs(x = "Average fishing hours", y = "Count") +
  facet_wrap( ~ year, ncol = 2)
```


### VMS Intervals

```{r table-ping-interval, results='asis'}
kable(do.call(rbind, tapply(ve$int, ve$year, summary)), booktabs = TRUE,
      caption = "Summary statistics of average ping interval")
```


```{r ping-interval, fig.height = 6}
#d <- 
#  ve %>% 
#  group_by(year) %>% 
#  summarise(effort = max(effort))
ggplot(ve, aes(x = int)) +
  geom_histogram() +
  scale_x_log10() +
  #geom_vline(data = d, aes(xintercept = effort)) +
  labs(x = "Average ping interval", y = "Count") +
  facet_wrap( ~ year, ncol = 2)
```


### Vessel length

```{r by_length, , fig.height = 6}
kable(do.call(rbind, tapply(ve$length, ve$year, summary)),
      caption = "Quantile distribution of mean vessel length")
ggplot(ve, aes(x = length))+
  geom_histogram() +
  labs(x = NULL, y = NULL,
       title = "Distribution of average vessel length") +
  facet_wrap( ~ year, ncol = 2)
```

### kW

```{r average_kW, fig.height = 6}
kable(do.call(rbind, tapply(ve$kw, ve$year, summary)),
      caption = "Quantile distribution of average kW")

ggplot(ve, aes(x = kw))+
  geom_histogram() +
  xlab("Average kW") + ylab ("Count") +
  facet_wrap( ~ year, ncol = 2)
```

### kW-hours

```{r kwh, fig.height = 6}
kable(do.call(rbind, tapply(ve$kwh, ve$year, summary)),
      caption = "Quantile distribution of average kW")
ggplot(ve, aes(x = kwh))+
  geom_histogram() +
  labs(x = NULL, y = "Count",
       title = "Average kw-hours") +
  facet_wrap(~ year) +
  scale_x_log10()
```

## Sums and medians {.tabset}

### Landings by Gear

```{r landings-by-target-year}
ve %>% 
  group_by(year, gear) %>% 
  summarise(catch = sum(catch) / 1000,
            .groups = "drop") %>% 
  spread(year, catch) %>% 
  kable(caption = "Landings [kt]")

ve %>% 
  group_by(year, gear) %>% 
  summarise(catch = sum(catch) / 1000,
            .groups = "drop") %>% 
  ggplot(aes(year, catch)) +
  geom_line(lwd = 1) +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = "Landings [kt]",
       title = "Catch by target")
```

### Landings by Gear - comparison with logbooks

```{r landings-by-target-year_lgs-comparision}
bind_rows(ve %>% 
            group_by(year, gear) %>% 
            summarise(catch = sum(catch) / 1000,
                      .groups = "drop") %>% 
            mutate(source = "vms"),
          le %>% 
            group_by(year, gear) %>% 
            summarise(catch = sum(catch) / 1000,
                      .groups = "drop") %>% 
            mutate(source = "logbook")) %>% 
  ggplot(aes(year, catch, colour = source)) +
  geom_line(lwd = 1) +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = "Landings [kt]",
       title = "Catch by target - comparison with logbooks") +
  expand_limits(y = 0)
```

### Effort by Gear

```{r effort-by-target-year}
ve %>% 
  group_by(year, gear) %>% 
  # NOTE: EFFORT SHOULD NOT BE NA, FIX UPSTREAM
  summarise(effort = sum(effort) / 24,
            .groups = "drop") %>% 
  spread(year, effort) %>% 
  kable(caption = "Effort [days]")

ve %>% 
  group_by(year, gear) %>% 
  summarise(effort = sum(effort) / 24,
            .groups = "drop") %>% 
  ggplot(aes(year, effort)) +
  geom_line(lwd = 1) +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = "Effort [days]",
       title = "Effort by target")
```

### Effort by Gear - comparison with logbooks

```{r effort-by-target-year_lgs-comparision}
bind_rows(ve %>% 
            group_by(year, gear) %>% 
            summarise(effort = sum(effort) / 24,
                      .groups = "drop") %>% 
            mutate(source = "vms"),
          le %>% 
            group_by(year, gear) %>% 
            summarise(effort = sum(effort),
                      .groups = "drop") %>% 
            mutate(source = "logbooks")) %>% 
  ggplot(aes(year, effort, colour = source)) +
  geom_line(lwd = 1) +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = "Effort [days]",
       title = "Effort [days] by target - comparison with logbooks",
       caption = "vms: Fishing activity in days, lgs: Days at sea") +
  expand_limits(y = 0) +
  scale_colour_brewer(palette = "Set1")
```

NOTE: The logbook (lgs) efforts is effectively days at sea (in the Icelandic case, days at sea when some a fishing activity was reported). Thus, logbook efforts is expected to be always higher than the vms effort that estimated the duration of a fishing activity.

### Median landings per hour

```{r landings-by-gear0}
ve %>% 
  group_by(year, gear) %>% 
  summarise(cph = median(catch/effort),
            .groups = "drop") %>% 
  spread(year, cph) %>% 
  kable(caption = "Median landing in per hour")
ve %>% 
  group_by(year, gear) %>% 
  summarise(cph = median(catch/effort),
            .groups = "drop") %>% 
  ggplot(aes(year, cph)) +
  geom_line() +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = NULL,
       title = "Median landing [kg] per hour")
```

### Median landing per kWh

```{r landings-by-target}
ve %>% 
  group_by(year, gear) %>% 
  summarise(cpkwh = median(catch/kwh),
            .groups = "drop") %>% 
  spread(year, cpkwh) %>% 
  kable(caption = "Median landing in kilograms per kilowhathour")
ve %>% 
  group_by(year, gear) %>% 
  summarise(cpkwh = median(catch/kwh),
            .groups = "drop") %>% 
  ggplot(aes(year, cpkwh)) +
  geom_line() +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = NULL,
       title = "Median landing [kg] per kilowhathour")
```

### Value

```{r value-by-target-year}
ve %>% 
  group_by(year, gear) %>% 
  summarise(value = sum(value),
            .groups = "drop") %>% 
  spread(year, value) %>% 
  kable(caption = "Value [EUR]")
ve %>% 
  group_by(year, gear) %>% 
  summarise(value = sum(value),
            .groups = "drop") %>% 
  ggplot(aes(year, value)) +
  geom_line(lwd = 1) +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = "[EUR]",
       title = "Value by target")
```


### Median value per kWh

```{r mean-value-by-KWhours-year}
ve %>% 
  group_by(year, gear) %>% 
  summarise(vpkwh = median(value/kwh, na.rm=T),
            .groups = "drop") %>% 
  spread(year, vpkwh) %>% 
  kable(caption = "Median value [EUR] per kilowhathour")
ve %>% 
  group_by(year, gear) %>% 
  summarise(vpkwh = median(value/kwh, na.rm=T),
            .groups = "drop") %>% 
  ggplot(aes(year, vpkwh)) +
  geom_line() +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = NULL,
       title = "Median value [EUR] per kilowhathour")
```

###  Median price per kg

```{r price}
ve %>% 
  group_by(year, gear) %>% 
  summarise(ppkg = median(value/catch),
            .groups = "drop") %>% 
  spread(year, ppkg) %>% 
  kable(caption = "Median value per kg")
ve %>% 
  group_by(year, gear) %>% 
  summarise(ppkg = median(value/catch),
            .groups = "drop") %>% 
  ggplot(aes(year, ppkg)) +
  geom_line() +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = NULL,
       title = "Median value per kg")
```


## Invalid ICES rectangles {.tabset}

```{r invalid stat-sqrs}
if (any(is.na(ve$lon))) {
  kable(table(`ICES Rectangle` = ve$isq[is.na(ve$lon)],
              Year = ve$year[is.na(ve$lon)]), booktabs = TRUE)
} else {
  x <- "There were no invalid Statistical Rectangles reported"
  attr(x, "format") <- "markdown"
  attr(x, "class") <- "knit_asis"
  x
}
```



## Space {.tabset}

### Overall

```{r table-exent}
ve %>% 
  group_by(year) %>% 
  summarise(min.lon = min(lon),
            max.lon = max(lon),
            min.lat = min(lat),
            max.lat = max(lat)) %>% 
  kable(caption = "Extent of vms data by year")
```

```{r area-of-data0}
d <- 
  ve %>% 
  select(year, lon, lat) %>% 
  distinct() %>% 
  count(lon, lat)
map1 +
  geom_tile(data = d, aes(lon, lat, fill = n)) +
  scale_fill_viridis_c(guide = "legend",
                       breaks = 1:50) +
  labs(fill = "Years",
       title = "Number of reported years per csquare")
```

```{r area-of-data,  fig.height = 9}
d <- 
  ve %>% 
  select(year, lon, lat) %>% 
  distinct()
map1 +
  geom_tile(data = d, aes(lon, lat), fill = "darkred") +
  scale_fill_viridis_c() +
  facet_wrap(~ year, ncol = 3) +
  labs(title = "Reported csquares per year")
```

```{r area-of-data2,  fig.height = 9}
map2 +
  geom_tile(data = d, aes(lon, lat), fill = "darkred") +
  scale_fill_viridis_c() +
  facet_wrap(~ year, ncol = 2) +
  labs(title = "Core area: Reported csquares per year")
```

### Effort

```{r spatial-effort-year, results='asis', fig.height = 9}
# do plots
core  <- list(xlim = quantile(ve$lon, c(0.025, 0.975)),
              ylim = quantile(ve$lat, c(0.025, 0.975)))
d <- 
  ve %>% 
  group_by(year, lon, lat) %>% 
  summarise(fishing_days = sum(effort) / 24,
            .groups = "drop")
map0 +
  geom_tile(data = d,
            aes(lon, lat, fill = fishing_days)) +
  facet_wrap(~ year, ncol = 2) +
  coord_quickmap(xlim = core$xlim, ylim = core$ylim) +
  scale_fill_viridis_c(trans = "log", option = "B", guide = "legend",
                       direction = -1,
                       breaks = c(10000, 2000, 500, 100,50, 10, 1, 0.5,0.1,0),
                       labels = c(10000, 2000, 500, 100,50, 10, 1, 0.5,0.1,0)) +
  labs(title = paste("Days@Sea"))
```

### Effort by Target

```{r spatial-extent-dominant-gears, fig.height = 9}
top4 <- 
  ve %>% 
  group_by(target) %>% 
  summarise(effort = sum(effort)) %>% 
  arrange(-effort) %>% 
  slice(1:6) %>% 
  pull(target) 
d <- 
  ve %>% 
  filter(target %in% top4) %>% 
  group_by(target, year, lon, lat) %>% 
  summarise(fishing_days = sum(effort) / 24,
            .groups = "drop")
# do plots
for (target in top4) {
  
  tmp <- d %>% filter(target == target)
  core  <- list(xlim = quantile(tmp$lon, c(0.025, 0.975)),
                ylim = quantile(tmp$lat, c(0.025, 0.975)))
  
  p <-
    map0 +
    geom_tile(data = d %>% filter(target == target),
              aes(lon, lat, fill = fishing_days)) +
    facet_wrap(~ year, ncol = 2) +
    coord_quickmap(xlim = core$xlim, ylim = core$ylim) +
    scale_fill_viridis_c(trans = "log", option = "B", guide = "legend",
                         direction = -1,
                         breaks = c(10000, 2000, 500, 100,50, 10, 1, 0.5,0.1,0),
                         labels = c(10000, 2000, 500, 100,50, 10, 1, 0.5,0.1,0)) +
    labs(title = paste("Days@Sea: ", target))
  print(p)
}
```

### Effort difference `r tempBound[1]`:`r (tempBound[2]-1)` vs `r tempBound[2]`

```{r, fig.height = 9}
base <-
  ve %>% 
  group_by(year, csq, lon, lat) %>% 
  summarise(effort = sum(effort),
            .groups = "drop")
recent <- 
  base %>% 
  filter(year == max(year)) %>% 
  select(-year)
d <- 
  base %>% 
  filter(year < max(year)) %>% 
  group_by(csq, lon, lat) %>% 
  summarise(median = median(effort),
            .groups = "drop") %>% 
  full_join(recent,
            by = c("csq", "lon", "lat")) %>% 
  mutate(effort = replace_na(effort, 0),
         median = replace_na(median, 0)) %>% 
  group_by(csq, lon, lat) %>% 
  mutate(r = 1 / (pmax(effort, 1e-9) / pmax(median, 1e-9))) %>% 
  ungroup() %>% 
  mutate(r = ifelse(r > 2, 2, r),
         r = ifelse(r < 1/2, 1/2, r))

ggplot() +
  geom_tile(data = d,
            aes(lon, lat, fill = r)) +
  coord_quickmap() +
  scale_fill_viridis_c(guide = "legend",
                       breaks = c(-Inf, 1/2, 1/1.5, 1/1.25, 1/1.05, 1, 1.05, 1.25, 1.5, 2, Inf),
                       labels = c("historic >>","historic> +100%","historic> +50%","historic> +25%",
                                  "+/-5%",
                                  "recent> +5%","recent> +25%","recent> +50%","recent> +100%","recent >>", "dummy"))
```


```{r spatial-effort-difference1, fig.height = 9}
base <- with(ve,
             aggregate(effort,
                       by = list(c_square = csq, year = year),
                       FUN = sum, na.rm = TRUE))
base <- dplyr::rename(base, fishing_hours = x)

# calculate total fishing hours for recent year
recent <- base[base$year == tempBound[2],]

# calculate median of the total fishing hours per square for historical years
#Josefine: only one year so doesnt work.
#base <- with(base[base$year < tempBound[2],],
#             aggregate(fishing_hours,
#                       by = list(c_square = c_square),
#                       FUN = median, na.rm = TRUE))

base <- with(base,
             aggregate(fishing_hours,
                       by = list(c_square = c_square),
                       FUN = median, na.rm = TRUE))

base <- dplyr::rename(base, fishing_hours_median = x)

# join
dat2plot <- dplyr::full_join(base,
                             recent[,c("c_square","fishing_hours")])

# set NAs to zero
dat2plot$fishing_hours_median[is.na(dat2plot$fishing_hours_median)] <- 0
dat2plot$fishing_hours[is.na(dat2plot$fishing_hours)] <- 0

# calculate ratio (with exceptions for zeros)
dat2plot$ratio <- 1/with(dat2plot, pmax(fishing_hours, 1e-9) / pmax(fishing_hours_median, 1e-9))

# add back in lat and long
dat2plot <- cbind(dat2plot,
                  CSquare2LonLat(dat2plot$c_square, degrees = 0.05))

breaks <- c(-Inf, 1/2, 1/1.5, 1/1.25, 1/1.05, 1, 1.05, 1.25, 1.5, 2, Inf)
legval <- c("historic >>","historic> +100%","historic> +50%","historic> +25%",
            "+/-5%",
            "recent> +5%","recent> +25%","recent> +50%","recent> +100%","recent >>")
palette <- RColorBrewer::brewer.pal(length(breaks)-1,"RdYlBu")  #colour for fishing intensity
dat2plot$colgrp <- as.numeric(cut(dat2plot$ratio, breaks = breaks))
dat2plot$cols <- palette[dat2plot$colgrp]
grps <- sort(unique(dat2plot$colgrp))
map1 +
  geom_tile(data = dat2plot,
            aes(SI_LONG, SI_LATI, fill = cols)) +
  #coord_quickmap() +
  scale_fill_manual(values = rev(palette)[grps], labels = legval[grps]) +
  labs(fill = "Days@Sea",
       subtitle = "Whole area",
       title = paste0("Effort difference: ", tempBound[1], ":", tempBound[2]-1, " vs ", tempBound[2]))
map2 +
  geom_tile(data = dat2plot,
            aes(SI_LONG, SI_LATI, fill = cols)) +
  #coord_quickmap() +
  scale_fill_manual(values = rev(palette)[grps], labels = legval[grps]) +
  labs(fill = "Days@Sea",
       subtitle = "Core area",
       title = paste0("Effort difference: ", tempBound[1], ":", tempBound[2]-1, " vs ", tempBound[2]))

```

### Effort difference `r tempBound[2]-1` vs `r tempBound[2]`

```{r, fig.height = 9}
base <- with(ve,
             aggregate(effort,
                       by = list(c_square = csq, year = year),
                       FUN = sum, na.rm = TRUE))
base <- dplyr::rename(base, fishing_hours = x)

# calculate total fishing hours for recent year
recent <- base[base$year == tempBound[2],]

# previous year
base <- base[base$year == tempBound[2]-1,]
base <- dplyr::rename(base, fishing_hours_median = fishing_hours)

# join
dat2plot <- dplyr::full_join(base,
                             recent[,c("c_square","fishing_hours")])


# set NAs to zero
dat2plot$fishing_hours_median[is.na(dat2plot$fishing_hours_median)] <- 0
dat2plot$fishing_hours[is.na(dat2plot$fishing_hours)] <- 0

# calculate ratio (with exceptions for zeros)
dat2plot$ratio <- 1/with(dat2plot, pmax(fishing_hours, 1e-9) / pmax(fishing_hours_median, 1e-9))

# add back in lat and long
dat2plot <- cbind(dat2plot,
                  CSquare2LonLat(dat2plot$c_square, degrees = 0.05))

breaks <- c(-Inf, 1/2, 1/1.5, 1/1.25, 1/1.05, 1, 1.05, 1.25, 1.5, 2, Inf)
legval <- c("historic >>","historic> +100%","historic> +50%","historic> +25%",
            "+/-5%",
            "recent> +5%","recent> +25%","recent> +50%","recent> +100%","recent >>")
palette <- RColorBrewer::brewer.pal(length(breaks)-1,"RdYlBu")  #colour for fishing intensity
dat2plot$colgrp <- as.numeric(cut(dat2plot$ratio, breaks = breaks))
dat2plot$cols <- palette[dat2plot$colgrp]
grps <- sort(unique(dat2plot$colgrp))
map1 +
  geom_tile(data = dat2plot,
            aes(SI_LONG, SI_LATI, fill = cols)) +
  scale_fill_manual(values = rev(palette)[grps], labels = legval[grps]) +
  labs(fill = "Days@Sea",
       subtitle = "Whole area",
       title = paste0("Effort difference: ", tempBound[2]-1, " vs ", tempBound[2]))
map2 +
  geom_tile(data = dat2plot,
            aes(SI_LONG, SI_LATI, fill = cols)) +
  scale_fill_manual(values = rev(palette)[grps], labels = legval[grps]) +
  labs(fill = "Days@Sea",
       subtitle = "Core area",
       title = paste0("Effort difference: ", tempBound[2]-1, " vs ", tempBound[2]))
```

## Check codes in VMS data {.tabset}

### Country code

Check VMS CountryCode with ICES vocabulary //vocab.ices.dk/?ref=337
```{r}
#Check Logbook CountryCode
CountryCode <- getCodeList("ISO_3166")
CountryCode$VocabCountry <- CountryCode$Key

VECountry <- left_join(ve, CountryCode, by = c("country" = "Key"))%>%
  select(country, VocabCountry) %>%
  distinct()

VECountryErr <- VECountry[is.na(VECountry$VocabCountry),]
datatable(VECountryErr)
```

### Gear type

Check VMS GearType with ICES vocabulary //vocab.ices.dk/?ref=1498
```{r}
#Check VMS GearType
GearType <- getCodeList("GearTypeL4")
GearType$VocabGear <- GearType$Key

VEGear <- left_join(ve, GearType, by = c("gear" = "Key"))%>%
  select(gear, VocabGear) %>%
  distinct()

VEGearErr <- VEGear[is.na(VEGear$VocabGear),]
datatable(VEGearErr) 
```

### Target assemblage

Check VMS TargetAssemblage with ICES vocabulary //vocab.ices.dk/?ref=1499
```{r}
#Check VMS TargetAssemblage
TargetAssemblage <- getCodeList("TargetAssemblage")
TargetAssemblage$VocabTargetAssemblage <- TargetAssemblage$Key

VETargetAssemblage <- left_join(ve, TargetAssemblage, by = c("target" = "Key"))%>%
  select(target, VocabTargetAssemblage) %>%
  distinct()

VETargetAssemblageErr <- VETargetAssemblage[is.na(VETargetAssemblage$VocabTargetAssemblage),]
datatable(VETargetAssemblageErr) 
```

### Vessel length

Check VMS VesselLength with ICES vocabulary //vocab.ices.dk/?ref=1502
```{r}
#Check VMS VesselLength
VesselLength <- getCodeList("BYC_VesselLRange")
VesselLength$VocabVesselLength <- VesselLength$Key

VEVesselLength <- left_join(ve, VesselLength, by = c("lclass" = "Key"))%>%
  select(lclass, VocabVesselLength) %>%
  distinct()

VEVesselLengthErr <- VEVesselLength[is.na(VEVesselLength$VocabVesselLength),]
datatable(VEVesselLengthErr) 
```


# Logbooks
___

## Counts {.tabset}

### By time


```{r logbook-records}
le %>% 
  count(year, month) %>% 
  mutate(month = as.character(str_pad(month, width = 2, pad = "0"))) %>% 
  add_row(le %>% 
            count(year) %>% 
            mutate(month = "Total")) %>% 
  spread(year, n) %>% 
  kable(caption = "Number of logbook records by year & months")
le %>% 
  count(year) %>% 
  ggplot(aes(year, n)) +
  geom_point() +
  scale_x_continuous(breaks = 2000:2030) +
  labs(x = NULL, y = NULL, title = "Number of logbook records by year") 
```

```{r logbook-entries-month}
le %>% 
  count(year, month) %>% 
  complete(year, month, fill = list(n = 0)) %>% 
  ggplot(aes(year, n)) +
  geom_point() +
  geom_line() +
  facet_wrap(~ month) +
  labs(x = NULL, y = NULL, title = "Number of logbook records by year & month")
```

### By length class

```{r by_lclass2}
le %>% 
  count(year, lclass) %>% 
  spread(year, n) %>% 
  kable(caption = "Number of logbook records by length class and year")
le %>% 
  count(year, lclass) %>% 
  complete(year, lclass, fill = list(n = 0)) %>% 
  left_join(lclass) %>% 
  unite(lclass, length_class, col = "lclass", sep = ": ") %>% 
  ggplot(aes(year, n)) +
  geom_point() +
  geom_line() +
  facet_wrap(~ lclass, scales = "free_y") +
  labs(x = NULL, y = NULL, title = "Number of logbook records by vessel length class")
```

### By Gear

```{r by_dcf42}
le %>% 
  group_by(year) %>% 
  summarise(gear = n_distinct(gear)) %>% 
  spread(year, gear) %>% 
  kable(caption = "Number of unique logbook Gear by year")
le %>% 
  count(year, gear) %>% 
  spread(year, n) %>% 
  kable(caption = "Number of logbook records by Gear")
le %>% 
  count(year, gear) %>% 
  ggplot(aes(year, n)) +
  geom_point() +
  geom_line() +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = NULL, title = "Number of logbook records by Gear") 
```

### By Target

```{r by_dcf52}
le %>% 
  group_by(year) %>% 
  summarise(target = n_distinct(target)) %>% 
  spread(year, target) %>% 
  kable(caption = "Number of unique logbooks Target by year")
le %>% 
  count(year, target) %>% 
  spread(year, n) %>% 
  kable(caption = "Number of logbooks records by Target")
le %>% 
  count(year, target) %>% 
  ggplot(aes(year, n)) +
  geom_point() +
  geom_line() +
  facet_wrap(~ target, scales = "free_y") +
  labs(x = NULL, y = NULL, title = "Number of logbook records by Target") 
```


### By Metier6

```{r by_dcf62, fig.height = 6}
le %>% 
  group_by(year) %>% 
  summarise(metier6 = n_distinct(metier6)) %>% 
  spread(year, metier6) %>% 
  kable(caption = "Number of unique logbook Metier6 by year")
le %>% 
  count(metier6, year) %>% 
  spread(year, n) %>% 
  kable(caption = "Number of logbook records by Metier6 by month and year")
# add here labels on the metier6 level
le %>% 
  count(year, gear, target, metier6) %>% 
  ggplot(aes(year, n, colour = target, group = metier6)) +
  geom_point() +
  geom_line() +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = NULL,
       title = "Number of logbook records by Metier6",
       subtitle = "needs further work")
```

### Unique vessels

```{r unique-vessels222}
le %>% 
  count(n_vessel) %>% 
  mutate(p = n / sum(n) * 1000,
         p = ifelse(p < 1.03, 1.03, p)) %>% 
  ggplot(aes(n_vessel, p)) +
  geom_col() +
  labs(x = "Number of unique vessels",
       title = "Unnagregated",
       subtitle = "Promill of ices-squares with given number of unique vessels") +
  scale_y_log10()
```

```{r unique-vessels22}
tmp <- le %>%
  dplyr::group_by(year, isq) %>%
  dplyr::summarise(
    new_count =
      ifelse(
        any(n_vessel > 2),
        3,
        length(unique(unlist(strsplit(vids, ";"))))
      )
  ) %>%
  dplyr::mutate(
    new_count = ifelse(new_count >= 3, "3+", paste(new_count))
  )

try(
  barplot(
    prop.table(
      table(factor(tmp$new_count))
    ),
    main = "aggregated to year level",
    ylab = "Proportion of ices-squares with given number of unique vessels",
    xlab = "number of unique vessels"
  )
)
```

```{r unique-vessels32}
tmp <- le %>%
  dplyr::group_by(isq) %>%
  dplyr::summarise(
    new_count =
      ifelse(
        any(n_vessel > 2),
        3,
        length(unique(unlist(strsplit(vids, ";"))))
      )
  ) %>%
  dplyr::mutate(
    new_count = ifelse(new_count >= 3, "3+", paste(new_count))
  )

try(
  barplot(
    prop.table(
      table(factor(tmp$new_count))
    ),
    main = "fully aggregated",
    ylab = "Proportion of ices-squares with given number of unique vessels",
    xlab = "number of unique vessels"
  )
)
```

## Distribution {.tabset}


### Fishing days

```{r table-fishing-day}
kable(do.call(rbind, tapply(le$effort, le$year, summary)), booktabs = TRUE,
      caption = "Summary statistics of average fishing days")
```

```{r fishing-days, fig.height = 6}
ggplot(le, aes(x = effort)) +
  geom_histogram() +
  scale_x_log10(breaks = c(1, 2, 3, 5, 10, 50, 100, 500)) +
  labs(x = "Average fishing days", y = "Count") +
  facet_wrap( ~ year, ncol = 2)
```


### kW-days

```{r average_kWdays, fig.height = 6}
kable(do.call(rbind, tapply(le$kwd, le$year, summary)),
      caption = "Quantile distribution of average kW days")

ggplot(le, aes(x = kwd))+
  geom_histogram() +
  scale_x_log10() +
  xlab("Average kW days") + ylab ("Count") +
  facet_wrap( ~ year, ncol = 2)
```

## Sums and medians {.tabset}

### Landings by Gear

```{r landings-by-target-year2}
le %>% 
  group_by(year, gear) %>% 
  summarise(catch = sum(catch) / 1000,
            .groups = "drop") %>% 
  spread(year, catch) %>% 
  kable(caption = "Landings [kt]")

le %>% 
  group_by(year, gear) %>% 
  summarise(catch = sum(catch) / 1000,
            .groups = "drop") %>% 
  ggplot(aes(year, catch)) +
  geom_line(lwd = 1) +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = "[kt]",
       title = "Landings by target")
```

### Median landing per kW-days

```{r landings-by-gear2}
le %>% 
  group_by(year, gear) %>% 
  summarise(cpkwh = median(catch/kwd),
            .groups = "drop") %>% 
  spread(year, cpkwh) %>% 
  kable(caption = "Median landing in kilograms per kilowhathour")
le %>% 
  group_by(year, gear) %>% 
  summarise(cpkwh = median(catch/kwd),
            .groups = "drop") %>% 
  ggplot(aes(year, cpkwh)) +
  geom_line() +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = NULL,
       title = "Median landing [kg] per kilowattdays")
```

### Value

```{r value-by-target-year2}
le %>% 
  group_by(year, gear) %>% 
  summarise(value = sum(value),
            .groups = "drop") %>% 
  spread(year, value) %>% 
  kable(caption = "Value [EUR]")
le %>% 
  group_by(year, gear) %>% 
  summarise(value = sum(value),
            .groups = "drop") %>% 
  ggplot(aes(year, value)) +
  geom_line(lwd = 1) +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = "[EUR]",
       title = "Value by target")
```

### Median value per kW-days

```{r mean-value-by-KWhours-year2}
le %>% 
  group_by(year, gear) %>% 
  summarise(vpkwh = median(value/kwd),
            .groups = "drop") %>% 
  spread(year, vpkwh) %>% 
  kable(caption = "Median value [EUR] per kilowatt-days")
le %>% 
  group_by(year, gear) %>% 
  summarise(vpkwh = median(value/kwd),
            .groups = "drop") %>% 
  ggplot(aes(year, vpkwh)) +
  geom_line() +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = NULL,
       title = "Median value [EUR] per kilowatt-days")
```


###  Median price per kg

```{r price2}
le %>% 
  group_by(year, gear) %>% 
  summarise(ppkg = median(value/catch, na.rm=T),
            .groups = "drop") %>% 
  spread(year, ppkg) %>% 
  kable(caption = "Median value per kg")
le %>% 
  group_by(year, gear) %>% 
  summarise(ppkg = median(value/catch, na.rm=T),
            .groups = "drop") %>% 
  ggplot(aes(year, ppkg)) +
  geom_line() +
  facet_wrap(~ gear, scales = "free_y") +
  labs(x = NULL, y = NULL,
       title = "Median value per kg")
```


## Invalid ICES rectangles {.tabset}

```{r invalid stat-sqrs-le}
if (any(is.na(le$lon))) {
  kable(table(`ICES Rectangle` = le$isq[is.na(le$lon)],
              Year = le$year[is.na(le$lon)]), booktabs = TRUE)
} else {
  x <- "There were no invalid Statistical Rectangles reported"
  attr(x, "format") <- "markdown"
  attr(x, "class") <- "knit_asis"
  x
}
```

## Space {.tabset}

### Overall

```{r table-exent2}
le %>% 
  group_by(year) %>% 
  summarise(min.lon = min(lon),
            max.lon = max(lon),
            min.lat = min(lat),
            max.lat = max(lat)) %>% 
  kable(caption = "Extent of vms data by year")
```

```{r area-of-data02}
d <- 
  le %>% 
  select(year, lon, lat) %>% 
  distinct() %>% 
  count(lon, lat)
map1 +
  geom_tile(data = d, aes(lon, lat, fill = n)) +
  scale_fill_viridis_c(guide = "legend",
                       breaks = 1:50) +
  labs(fill = "Years",
       title = "Number of reported years per ices square")
```

```{r area-of-data3,  fig.height = 9}
d <- 
  le %>% 
  select(year, lon, lat) %>% 
  distinct()
map1 +
  geom_tile(data = d, aes(lon, lat), fill = "darkred") +
  scale_fill_viridis_c() +
  facet_wrap(~ year, ncol = 3) +
  labs(title = "Reported ices squares per year")
```

```{r area-of-data22,  fig.height = 9}
map2 +
  geom_tile(data = d, aes(lon, lat), fill = "darkred") +
  scale_fill_viridis_c() +
  facet_wrap(~ year, ncol = 2) +
  labs(title = "Core area: Reported ices squares per year")
```

### Effort

```{r spatial-effort-year2, results='asis', fig.height = 9}
# do plots
core  <- list(xlim = quantile(le$lon, c(0.025, 0.975)),
              ylim = quantile(le$lat, c(0.025, 0.975)))
d <- 
  le %>% 
  group_by(year, lon, lat) %>% 
  summarise(fishing_days = sum(effort),
            .groups = "drop")
map0 +
  geom_tile(data = d,
            aes(lon, lat, fill = fishing_days)) +
  facet_wrap(~ year, ncol = 2) +
  coord_quickmap(xlim = core$xlim, ylim = core$ylim) +
  scale_fill_viridis_c(trans = "log", option = "B", guide = "legend",
                       direction = -1,
                       breaks = c(10000, 2000, 500, 100,50, 10, 1, 0.5,0.1,0),
                       labels = c(10000, 2000, 500, 100,50, 10, 1, 0.5,0.1,0)) +
  labs(title = paste("Days@Sea"))
```

### Effort by Target

```{r spatial-extent-dominant-gears2, fig.height = 9}
top4 <- 
  le %>% 
  group_by(target) %>% 
  summarise(time = sum(effort)) %>% 
  arrange(-time) %>% 
  slice(1:6) %>% 
  pull(target) 
d <- 
  le %>% 
  filter(target %in% top4) %>% 
  group_by(target, year, lon, lat) %>% 
  summarise(fishing_days = sum(effort),
            .groups = "drop")
# do plots
for (target in top4) {
  
  tmp <- d %>% filter(target == target)
  core  <- list(xlim = quantile(tmp$lon, c(0.025, 0.975)),
                ylim = quantile(tmp$lat, c(0.025, 0.975)))
  
  p <-
    map0 +
    geom_tile(data = d %>% filter(target == target),
              aes(lon, lat, fill = fishing_days)) +
    facet_wrap(~ year, ncol = 2) +
    coord_quickmap(xlim = core$xlim, ylim = core$ylim) +
    scale_fill_viridis_c(trans = "log", option = "B", guide = "legend",
                         direction = -1,
                         breaks = c(10000, 2000, 500, 100,50, 10, 1, 0.5,0.1,0),
                         labels = c(10000, 2000, 500, 100,50, 10, 1, 0.5,0.1,0)) +
    labs(title = paste("Days@Sea: ", target))
  print(p)
}
```

### Effort difference `r tempBound[1]`:`r (tempBound[2]-1)` vs `r tempBound[2]`

```{r, eval = FALSE, fig.height = 9}
base <-
  le %>% 
  group_by(year, isq, lon, lat) %>% 
  summarise(time = sum(effort),
            .groups = "drop")
recent <- 
  base %>% 
  filter(year == max(year)) %>% 
  select(-year)
d <- 
  base %>% 
  filter(year < max(year)) %>% 
  group_by(isq, lon, lat) %>% 
  summarise(median = median(time),
            .groups = "drop") %>% 
  full_join(recent,
            by = c("isq", "lon", "lat")) %>% 
  mutate(time = replace_na(time, 0),
         median = replace_na(median, 0)) %>% 
  group_by(isq, lon, lat) %>% 
  mutate(r = 1 / (pmax(time, 1e-9) / pmax(median, 1e-9))) %>% 
  ungroup() %>% 
  mutate(r = ifelse(r > 2, 2, r),
         r = ifelse(r < 1/2, 1/2, r))

ggplot() +
  geom_tile(data = d,
            aes(lon, lat, fill = r)) +
  coord_quickmap() +
  scale_fill_viridis_c(guide = "legend",
                       breaks = c(-Inf, 1/2, 1/1.5, 1/1.25, 1/1.05, 1, 1.05, 1.25, 1.5, 2, Inf),
                       labels = c("historic >>","historic> +100%","historic> +50%","historic> +25%",
                                  "+/-5%",
                                  "recent> +5%","recent> +25%","recent> +50%","recent> +100%","recent >>", "dummy"))
```


```{r spatial-effort-difference12, fig.height = 9}
base <- with(le,
             aggregate(effort,
                       by = list(c_square = isq, year = year),
                       FUN = sum, na.rm = TRUE))
base <- dplyr::rename(base, fishing_hours = x)

# calculate total fishing hours for recent year
recent <- base[base$year == tempBound[2],]

# calculate median of the total fishing hours per square for historical years
#Josefine: Only one year
#base <- with(base[base$year < tempBound[2],],
#             aggregate(fishing_hours,
#                       by = list(c_square = c_square),
#                       FUN = median, na.rm = TRUE))
base <- with(base,
             aggregate(fishing_hours,
                       by = list(c_square = c_square),
                       FUN = median, na.rm = TRUE))


base <- dplyr::rename(base, fishing_hours_median = x)

# join
dat2plot <- dplyr::full_join(base,
                             recent[,c("c_square","fishing_hours")])

# set NAs to zero
dat2plot$fishing_hours_median[is.na(dat2plot$fishing_hours_median)] <- 0
dat2plot$fishing_hours[is.na(dat2plot$fishing_hours)] <- 0

# calculate ratio (with exceptions for zeros)
dat2plot$ratio <- 1/with(dat2plot, pmax(fishing_hours, 1e-9) / pmax(fishing_hours_median, 1e-9))

# add back in lat and long
dat2plot <- 
  dat2plot %>% 
  mutate(lon = ir2d(c_square)$lon,
         lat = ir2d(c_square)$lat)

breaks <- c(-Inf, 1/2, 1/1.5, 1/1.25, 1/1.05, 1, 1.05, 1.25, 1.5, 2, Inf)
legval <- c("historic >>","historic> +100%","historic> +50%","historic> +25%",
            "+/-5%",
            "recent> +5%","recent> +25%","recent> +50%","recent> +100%","recent >>")
palette <- RColorBrewer::brewer.pal(length(breaks)-1,"RdYlBu")  #colour for fishing intensity
dat2plot$colgrp <- as.numeric(cut(dat2plot$ratio, breaks = breaks))
dat2plot$cols <- palette[dat2plot$colgrp]
grps <- sort(unique(dat2plot$colgrp))
map1 +
  geom_tile(data = dat2plot,
            aes(lon, lat, fill = cols)) +
  #coord_quickmap() +
  scale_fill_manual(values = rev(palette)[grps], labels = legval[grps]) +
  labs(fill = "Days@Sea",
       subtitle = "Whole area",
       title = paste0("Effort difference: ", tempBound[1], ":", tempBound[2]-1, " vs ", tempBound[2]))
map2 +
  geom_tile(data = dat2plot,
            aes(lon, lat, fill = cols)) +
  #coord_quickmap() +
  scale_fill_manual(values = rev(palette)[grps], labels = legval[grps]) +
  labs(fill = "Days@Sea",
       subtitle = "Core area",
       title = paste0("Effort difference: ", tempBound[1], ":", tempBound[2]-1, " vs ", tempBound[2]))

```

### Effort difference `r tempBound[2]-1` vs `r tempBound[2]`

```{r, fig.height = 9}
base <- with(le,
             aggregate(effort,
                       by = list(c_square = isq, year = year),
                       FUN = sum, na.rm = TRUE))
base <- dplyr::rename(base, fishing_hours = x)

# calculate total fishing hours for recent year
recent <- base[base$year == tempBound[2],]

# previous year
base <- base[base$year == tempBound[2]-1,]
base <- dplyr::rename(base, fishing_hours_median = fishing_hours)

# join
dat2plot <- dplyr::full_join(base,
                             recent[,c("c_square","fishing_hours")])


# set NAs to zero
dat2plot$fishing_hours_median[is.na(dat2plot$fishing_hours_median)] <- 0
dat2plot$fishing_hours[is.na(dat2plot$fishing_hours)] <- 0

# calculate ratio (with exceptions for zeros)
dat2plot$ratio <- 1/with(dat2plot, pmax(fishing_hours, 1e-9) / pmax(fishing_hours_median, 1e-9))

# add back in lat and long
dat2plot <- 
  dat2plot %>% 
  mutate(lon = ir2d(c_square)$lon,
         lat = ir2d(c_square)$lat)
breaks <- c(-Inf, 1/2, 1/1.5, 1/1.25, 1/1.05, 1, 1.05, 1.25, 1.5, 2, Inf)
legval <- c("historic >>","historic> +100%","historic> +50%","historic> +25%",
            "+/-5%",
            "recent> +5%","recent> +25%","recent> +50%","recent> +100%","recent >>")
palette <- RColorBrewer::brewer.pal(length(breaks)-1,"RdYlBu")  #colour for fishing intensity
dat2plot$colgrp <- as.numeric(cut(dat2plot$ratio, breaks = breaks))
dat2plot$cols <- palette[dat2plot$colgrp]
grps <- sort(unique(dat2plot$colgrp))
map1 +
  geom_tile(data = dat2plot,
            aes(lon, lat, fill = cols)) +
  scale_fill_manual(values = rev(palette)[grps], labels = legval[grps]) +
  labs(fill = "Days@Sea",
       subtitle = "Whole area",
       title = paste0("Effort difference: ", tempBound[2]-1, " vs ", tempBound[2]))
map2 +
  geom_tile(data = dat2plot,
            aes(lon, lat, fill = cols)) +
  scale_fill_manual(values = rev(palette)[grps], labels = legval[grps]) +
  labs(fill = "Days@Sea",
       subtitle = "Core area",
       title = paste0("Effort difference: ", tempBound[2]-1, " vs ", tempBound[2]))
```


## Check codes in logbook data {.tabset}

### ICES rectangle

Check Logbook ICES rectangle with ICES vocabulary //vocab.ices.dk/?ref=107
```{r}
#Check Logbook ICES rectangle
icesRect <- getCodeList("StatRec")
icesRect$VocabRect <- icesRect$Key

LERect <- left_join(le, icesRect, by = c("isq" = "Key"))%>%
  select(isq, VocabRect) %>%
  distinct()

LERectErr <- LERect[is.na(LERect$VocabRect),]
datatable(LERectErr) 
```

### Country code

Check Logbook CountryCode with ICES vocabulary //vocab.ices.dk/?ref=337
```{r}
#Check Logbook CountryCode
CountryCode <- getCodeList("ISO_3166")
CountryCode$VocabCountry <- CountryCode$Key

LECountry <- left_join(le, CountryCode, by = c("country" = "Key"))%>%
  select(country, VocabCountry) %>%
  distinct()

LECountryErr <- LECountry[is.na(LECountry$VocabCountry),]
datatable(LECountryErr)
```

### Gear type

Check Logbook GearType with ICES vocabulary //vocab.ices.dk/?ref=1498
```{r}
#Check Logbook GearType
GearType <- getCodeList("GearTypeL4")
GearType$VocabGear <- GearType$Key

LEGear <- left_join(le, GearType, by = c("gear" = "Key"))%>%
  select(gear, VocabGear) %>%
  distinct()

LEGearErr <- LEGear[is.na(LEGear$VocabGear),]
datatable(LEGearErr) 
```

### Target assemblage

Check Logbook TargetAssemblage with ICES vocabulary //vocab.ices.dk/?ref=1499
```{r}
#Check Logbook TargetAssemblage
TargetAssemblage <- getCodeList("TargetAssemblage")
TargetAssemblage$VocabTargetAssemblage <- TargetAssemblage$Key

LETargetAssemblage <- left_join(le, TargetAssemblage, by = c("target" = "Key"))%>%
  select(target, VocabTargetAssemblage) %>%
  distinct()

LETargetAssemblageErr <- LETargetAssemblage[is.na(LETargetAssemblage$VocabTargetAssemblage),]
datatable(LETargetAssemblageErr) 
```

### Vessel length

Check Logbook VesselLength with ICES vocabulary //vocab.ices.dk/?ref=1502
```{r}
#Check Logbook VesselLength
VesselLength <- getCodeList("BYC_VesselLRange")
VesselLength$VocabVesselLength <- VesselLength$Key

LEVesselLength <- left_join(le, VesselLength, by = c("lclass" = "Key"))%>%
  select(lclass, VocabVesselLength) %>%
  distinct()

LEVesselLengthErr <- LEVesselLength[is.na(LEVesselLength$VocabVesselLength),]
datatable(LEVesselLengthErr) 
```

### VMSEnabled

Check Logbook VMSEnabled with ICES vocabulary //vocab.ices.dk/?ref=316
```{r}
#Check Logbook VMSEnabled
VMSEnabled <- getCodeList("YesNoFields")
VMSEnabled$VocabVMSEnabled <- VMSEnabled$Key

LEVMSEnabled <- left_join(le, VMSEnabled, by = c("vms" = "Key"))%>%
  select(vms, VocabVMSEnabled) %>%
  distinct()

LEVMSEnabledErr <- LEVMSEnabled[is.na(LEVMSEnabled$VocabVMSEnabled),]
datatable(LEVMSEnabledErr) 
```


# Comparisons
___


```{r}
le %>% count(year, metier6, name = "logbook records") %>% 
  full_join(ve %>% count(year, metier6, name = "vms records")) %>% 
  kable(caption = "Comparison of Metier level 6 reporting between logbook records and VMS records")
```