---
title: "Greater North Sea mixoplankton"
knit: (function(input_file, encoding) {
  out_dir <- '../docs';
  rmarkdown::render(input_file,
 encoding=encoding,
 output_file=file.path(dirname(input_file), out_dir, 'index.html'))})
author: "Willem Stolte"
date: "17/06/2020"
output:
  html_document:
    code_folding: hide
    df_print: paged
subtitle: A first analysis
---


```{r setup, include=FALSE}

knitr::opts_chunk$set(
  echo = TRUE, comment = F,
  # comment = FALSE, 
  warning = FALSE
)

require(sf)
require(tidyverse)
require(lubridate)
require(raster)
require(EMODnetBiologyMaps)
dataDir <- "../data/derived_data"
productDir <- "../product"
csvDir <- file.path(productDir, "csv")
pngDir <- file.path(productDir, "grid_plots")
tiffDir <- file.path(productDir, "tiff")
select <- dplyr::select
```

```{r styleDef}
EMODnetMapCleanStyle <- theme_bw() +
  theme(axis.text = element_blank(),
        axis.title = element_blank())

EMODnetMapNormalStyle <- 
  theme(axis.text = element_text(size = 12, color = "blue"),
        axis.title = element_text(size = 14, color = "blue"),
        legend.text = element_text(size = 10),
        legend.position = "right",
        legend.key.width = unit(10,"mm"),
        legend.key.height = unit(10, "mm"),
        plot.background = element_rect(fill = "grey90", colour = "white", size = 0),
        panel.border = element_blank()
)
```


# Introduction

## Documentation

Github project: https://github.com/wstolte/mixoplanktonFractions

# Data extraction

## Temporal scale

Data from 1995 until now are considered for the current product.

## Geographical scale

The regions that were selected were assembled from the intersection of the IHO regions and the EEZ from the different countries. These subregions have ID's that can be used in the WFS query to the EMODnet Biology database.

```{r mapGreaterNorthSea, message=F, warning=F, comment=F, out.width="50%"}
regions <- sf::st_read(quiet = T, dsn = "../data/derived_data/greater_north_sea-selection_from_eez-iho_v4.geojson")
regions %>% ggplot() +
  geom_sf(fill = "blue", color = "white") +
  geom_sf_text(aes(label = mrgid), size = 2, color = "white") +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())
```

## Strategy

For each of these subregions, data extraction was done in a two-step approach. 

1. Extract phytoplankton data by selecting only observations linkted to the trait "phytoplankton".
1. The unique datasets obtained in the first steps were inspected, and suitable datasets, which were expected to contain phytoplankton, were selected. The suitable datasets were then downloaded completely (without selection for traits)


## Extraction of complete relevant datasets

Data were extracted from EMODnet Biology via WFS, using the requestData.R script. These data are stored locally and taken into this script here. 


```{r, message=F, warning=F}
# all2Data <- read_delim(file.path(dataDir, "all2Data.csv"), delim = ";")
load(file.path(dataDir, "all2Data.Rdata"))
```

In total, this resulted in `r length(all2Data$id)` observations

Per subregion, the number of observations/km^2^ looks like this:

```{r, message=F, warning=F}
regionN <- all2Data %>%
  group_by(mrgid) %>% summarize(n = n()) %>% ungroup() %>%
  mutate(mrgid = as.numeric(mrgid))

regions %>% right_join(regionN, by = c(mrgid = "mrgid")) %>%
  mutate(n_km2 = n/area_km2) %>%
  ggplot() +
  geom_sf(aes(fill = n_km2)) +
  scale_fill_viridis_c()
```


## Join with trophic mode table

Schneider et al. (submitted to Journal of Biological Data, data submitted to WoRMS) has defined trophic mode for about 1500 species of plankton. A first subdivision is in 

* phytoplankton - phototrophic, no phagotrophy known
* mixoplankton - capable of phototrophy and phagotrophy
* zooplankton - phagotrophic, no phototrophy known.

Mixoplankton is subdivided in different classes, but at this moment, we will not use this. 


```{r}
# read information available in WoRMS traits. 
# This is a local copy

trophy <- read_delim("../data/raw_data/traits.csv", delim = ";", guess_max = 10000)

```



```{r}
trophyData <- all2Data %>% 
  mutate(AphiaID = sub("http://marinespecies.org/aphia.php[?]p=taxdetails&id=", "", aphiaid, )) %>% 
  full_join(trophy, by = c(AphiaID = "AphiaID")) %>% 
  select(
    mrgid, datasetID, datecollected, decimallongitude, decimallatitude,
    scientificname, AphiaID, Trophy, typeMX
  ) %>%
  drop_na(Trophy) %>%
  mutate(
    season = case_when(
      month(datecollected) %in% c(6,7,8) ~ "summer", 
      month(datecollected) %in% c(3,4,5) ~ "spring", 
      month(datecollected) %in% c(9,10,11) ~ "autumn", 
      month(datecollected) %in% c(12, 1, 2) ~ "winter"
    )
  ) %>% drop_na(season) %>%
  mutate(season = factor(season, levels = c("winter", "spring", "summer", "autumn")))


```

# Fraction of mixoplankton (first results)

As a first approach, the fraction of mixoplankton/(phytoplankto + mixoplankton) is plotted on a map:


```{r trophyMap, fig.width=10, fig.height=10, out.width="100%"}
plottrophy <- trophyData %>%
  group_by(decimallongitude, decimallatitude, Trophy, season) %>%
  summarize(n = n()) %>%
  filter(Trophy %in% c("mixotroph", "phototroph")) %>%
  spread(Trophy, n) %>%
  mutate(fMix = mixotroph/(mixotroph + phototroph)) %>%
  drop_na(fMix)

write_delim(plottrophy, file.path(csvDir, "trophydata.csv"), delim = ";")

require(rworldxtra)
data(countriesHigh)
world <- st_as_sf(countriesHigh) %>% 
  filter(REGION == "Europe") %>%
  filter(ADMIN %in% c("Denmark", "Germany", "Netherlands", "Belgium", "United Kingdom", "France", "Norway", "Sweden", "Ireland")) %>%
  select(ADMIN) %>%
  st_crop(xmin = -15, xmax = 17, ymin = 45, ymax = 65)


ggplot() +
  stat_summary_2d(data = plottrophy, aes(x = decimallongitude, y = decimallatitude, z=fMix), fun=mean, binwidth = c(1,0.5)) +
  geom_sf(data = world, fill = "transparent", color = "grey", alpha = 0) +
  scale_fill_viridis_c(option = "viridis") +
  coord_sf(xlim = c(-15, 17), ylim = c(45, 65), expand = T, clip = "on") +
  facet_wrap(~ season) +
  EMODnetMapCleanStyle
```


This first analysis shows highest fraction of mixoplankton in summer. Even in summer, the fraction of mixoplankton is still low along the Dutch, Belgian and French coast. This is in line with Schneider et al. (2020a). It is also in line with a general idea that mixoplankton, at least the constitutive mixoplankton, appear in post-bloom clear water conditions, which in the North Sea are associated with offshore, possibly stratified waters. 


```{r include=FALSE}

  proWG <- CRS("EPSG:4326")
  proUTM <- CRS("+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs")
  r1<-raster(ext=extent(-16,9,46,66),ncol=100,nrow=160,crs=proWG,vals=0)
  r2<-raster(ext=extent(-16,9,46,66),ncol=25,nrow=40,crs=proWG,vals=0)
  
  seasons <- c('winter', 'spring', 'summer', 'autumn')
  
  for(jj in c(1:length(seasons))){
    
    seasonalData <- plottrophy %>% dplyr::filter(season == seasons[jj])
    
    coordinates(seasonalData)<- ~decimallongitude+decimallatitude
    projection(seasonalData)<-proWG # because data were saved in UTM. Could be changed
    # seasonalData <- sp::spTransform(seasonalData, proWG)
    sp_r_highres<-rasterize(seasonalData,r1,field="fMix",fun=mean)
    sp_r_lowres<-rasterize(seasonalData,r2,field="fMix",fun=mean)
    
    #save tiffs
    
    ec<-emodnet_colors()
    spAphId <- "fMix"
    
    tifnam <- file.path(productDir, "tiff", paste(today(), gsub(" ", "-", spAphId), seasons[jj], ".tif", sep = "_"))
    tifnam_lowres <- file.path(productDir, "tiff", paste(today(), gsub(" ", "-", spAphId), seasons[jj], "_lowres.tif", sep = "_"))
    
    writeRaster(sp_r_highres, tifnam, options=c('TFW=YES'), overwrite = T)
    writeRaster(sp_r_lowres, tifnam_lowres, options=c('TFW=YES'), overwrite = T)
    
    plot_grid_highres <- emodnet_map_plot(data=sp_r_highres,title=paste0(spAphId),subtitle=paste("probability of occurrence in", seasons[jj]),
                                  zoom=FALSE,seaColor=ec$darkgrey,landColor=ec$lightgrey,legend="",
                                  xlim = c(2351321, 4275244), ylim = c(2543530,4985495)) +
      scale_fill_viridis_c(limits = c(0,0.6))
    
    filnam <- file.path(productDir, "grid_plots", paste(today(), gsub(" ", "-", spAphId), seasons[jj], "highres.png", sep = "_"))
    
    emodnet_map_logo(plot_grid_highres,path=filnam,width=120,height=160,dpi=300,units="mm",offset="+0+0")
    
    
    plot_grid_lowres <- emodnet_map_plot(data=sp_r_lowres,title=paste0(spAphId),subtitle=paste("probability of occurrence in", seasons[jj]),
                                         zoom=FALSE,seaColor=ec$darkgrey,landColor=ec$lightgrey,legend="",
                                         xlim = c(2351321, 4275244), ylim = c(2543530,4985495)) +
      scale_fill_viridis_c(limits = c(0,0.6))
    
    filnam_lowres <- file.path(productDir, "grid_plots", paste(today(), gsub(" ", "-", spAphId), seasons[jj], "_lowres.png", sep = "_"))
    
    emodnet_map_logo(plot_grid_lowres,path=filnam_lowres,width=120,height=160,dpi=300,units="mm",offset="+0+0")
    
  }
```



# Limitations of this analysis - discussion

## Incomplete classification

At the moment, only 1500 species could be automatically classified as either phytoplankton or mixoplankton. This means that we are missing a number of species in this analysis, because the data contain more than 5000 species names. Part of those may be heterotrophic, which does not influence the maps above. If we assume that the fraction of mixoplankton is not very different in these unknown species as compared to the known species, the bias is not very strong. In the current analysis, the number of observations where species could not be classified was lower than the number of phytoplankton observations (figure ]\@ref(fig:. 

```{r classifications, fig.cap="Distribution of observations over known mixotrophs, phagotrophs, phototrophs and unknown trophy according to Schneider et al., (2020b). "}
all2Data %>% 
  mutate(AphiaID = sub("http://marinespecies.org/aphia.php[?]p=taxdetails&id=", "", aphiaid, )) %>% 
  full_join(trophy, by = c(AphiaID = "AphiaID")) %>% 
  select(
    mrgid, datasetID, datecollected, decimallongitude, decimallatitude,
    scientificname, AphiaID, Trophy, typeMX
  ) %>%
  group_by(Trophy) %>%
  summarize(`nr of observations` = n()) %>%
  ggplot(aes(Trophy, `nr of observations`)) + geom_col(aes(fill = Trophy))
```



## Bias in species per dataset

Not all datasets address the same species. For example, the Plankton Recorder data contain species that are retained on a 60 µm filter, which excludes a number of small species. A possible way to get a better grip of the effect of this bias is to only regards species that occur in all (or most) datasets. This has not been undertaken yet.  


# References 

*   Schneider LK, Flynn KJ, Herman PMJ, Troost TA, Stolte W. (2020a) Exploring the trophic spectrum: placing mixoplankton into marine protist communities of the Southern North Sea. Frontiers in Marine Science 7: 997 -  https://doi.org/10.3389/fmars.2020.586915
*   Schneider LK, Anestis K, Mansour J, Anschütz AA, Gypens N, Hansen PJ, John U, Klemm K, Lapeyra Martin J, Medic N, Not F, Stolte W. (2020b) A dataset on trophic modes of aquatic protists. Biodiversity Data Journal 8. https://doi.org/10.3897/BDJ.8.e56648.







