WAP-Heatwaves_KrillRecruitment_Chlorophylla_PART-II.Rmd

---
title: Detection of extreme heatwaves and low sea ice cover in the Antarctic Peninsula  for evaluation of effects over Krill, Part II - process and analyse downloaded data
author: "Lucas Krüger, Maurício Mardones, Lorena Rebolledo"
date: "`r Sys.Date()`"
output: html_document
---

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

#---------- Download NOAA erddap optimum interpolation SST----

#this part of the script is a modification of the scripts in the heatwaveR vignette
#available at https://robwschlegel.github.io/heatwaveR/articles/OISST_preparation.html


```{r, echo =FALSE}
#---------- Download NOAA erddap optimum interpolation SST----

# The packages we will use
library(dplyr) # A staple for modern data management in R
library(lubridate) # Useful functions for dealing with dates
library(ggplot2) # The preferred library for data visualisation
library(patchwork)
library(tidync) # For easily dealing with NetCDF data
library(rerddap) # For easily downloading subsets of data
library(doParallel) # For parallel processing
library(raster)
library(terra)
library(sf)

library(mapdata)
library(mapproj)

library(reshape2)

```


###----------- calculating marine heatwaves-------------

### load and process krill data

```{r}

  
krill.strata<-readRDS("KrillData.Rdata") 


krill.strata$timestamp<-as.POSIXct(strptime(krill.strata$datetime_set_start, format="%Y-%m-%d %H:%M:%S", tz="GMT"))

krill.strata$Week<-week(krill.strata$timestamp)
krill.strata$Year<-year(krill.strata$timestamp)
krill.strata$Month<-month(krill.strata$timestamp)
krill.strata$Quarter<-quarter(krill.strata$timestamp)

krill.strata$Season[krill.strata$Quarter=="1"]<-"JFM"
krill.strata$Season[krill.strata$Quarter=="2"]<-"AMJ"
krill.strata$Season[krill.strata$Quarter=="3"]<-"JAS"
krill.strata$Season[krill.strata$Quarter=="4"]<-"OND"

krill.strata$sector<-krill.strata$ID


krill.strata$stage<-as.factor(krill.strata$maturity_stage)

krist<-na.omit(data.frame(krill.strata[13],krill.strata[,16:23]))

#correct strata name to match the climatology data

krist$sector[krist$sector=="Extra"]<-"EX"
krist$sector[krist$sector=="SSIW"]<-"SS"

krist$Talla<-ifelse(krist$length_total_cm<4,"Recruit","Adult") # first step to calculate the recruitment index


ad<-(plyr::ddply(krist, c("Year","sector"), summarise,
                        Ntot=length(na.omit(length_total_cm))))  # number of adults

rec<-(plyr::ddply(subset(krist,Talla=="Recruit"), c("Year","sector"), summarise,
                        Nrec=length(na.omit(length_total_cm))))  # number of recruiters


recind<-merge(ad,rec,all.x=TRUE)  # merge to calculate recruitment index

recind$Nrec[is.na(recind$Nrec)]<-0


recind$recind<-recind$Nrec/recind$Ntot  # recruitment index

summary(recind$Year)


ggplot(recind,aes(Year,recind))+
  geom_smooth(se=F,method="gam")+
  geom_point()+
  #facet_wrap(sector~.)+
  theme_bw()+ylab("Proportion of recruiting individuals")+xlim(2000,2015)

```



```{r}

chlall<-base::readRDS("CHL_summer.Rds")

climS<-base::readRDS("climatologyICE.Rds")
climT<-base::readRDS("climatologyTEMP.Rds")

summary((year(climS$t)))

climS$Year<-year(climS$t)
climS$Month<-month(climS$t)
climS$Quarter<-quarter(climS$t)

climT$Year<-year(climT$t)
climT$Month<-month(climT$t)
climT$Quarter<-quarter(climT$t)



climWS<-plyr::ddply(subset(climS,Quarter=="2"|Quarter=="3"), c("Year","sector"), summarise,
                    SICPeak=max(na.omit(sicMean)),SICmean=mean(na.omit(sicMean)),Sthresh=mean(na.omit(thresh)))

climST<-plyr::ddply(subset(climT,Quarter=="1"|Quarter=="4"), c("Year","sector"), summarise,
                    
                    TempMax=max(na.omit(tMax)),TempMean=mean(na.omit(tMax)),Tthresh=mean(thresh))

climWS<-merge(climWS,climST)

head(climWS)

evesum<-base::readRDS("events_summary.Rds")

eveclim<-merge(evesum,climWS,all=T)

mean(chlall$chlm)

ecc<-merge(eveclim,chlall,all=T)

chlk<-merge(ecc,recind)


chlk<-merge(ecc,recind)

head(chlk)
summary(chlk$SICmean)
summary(chlk$SICmean[chlk$Year<2005])
summary(chlk$SICmean[chlk$Year>2004])



chlk$sicc[chlk$SICmean<0.47 & chlk$Year>=2005]<-"Low"
chlk$sicc[chlk$SICmean<0.15 & chlk$Year<2005]<-"Low"

chlk$sicc[is.na(chlk$sicc)]<-"High"

summary(chlk$TempMean[chlk$Year<2005])
summary(chlk$TempMean[chlk$Year>2004])


chlk$tempc[chlk$TempMean<0.877 & chlk$Year>=2005]<-"Low"
chlk$tempc[chlk$TempMean<0.541 & chlk$Year<2005]<-"Low"

chlk$tempc[is.na(chlk$tempc)]<-"High"


ggplot()+
  geom_smooth(data=chlk,aes(x=Year,y=recind),
              method="gam",se=F,colour="grey30")+xlim(2000,2015)+
  geom_point(data=chlk,aes(x=Year,y=recind,colour=sicc,shape=sicc),size=2)+
  scale_colour_manual(values=c("blue2","red2"))+theme_bw()
  
  




chlk$TN[is.na(chlk$TN)]<-0

chlk$mhws<-chlk$TN
chlk$MHW<-ifelse(chlk$mhws>0,"with events","without events")

summary(chlk$chlm)

chlk$chlc<-ifelse(chlk$chlm<median(na.omit(chlk$chlm)),"Low","High")
chlk$chlc[is.na(chlk$chlc)]<-"Low"

median(chlk$recind)

kdf<-(data.frame(CHL=chlk$chlc,SIC=chlk$sicc,SST=chlk$tempc,MHW=chlk$MHW,
                KRI=ifelse(chlk$recind<0.16,"Low KRI","High KRI")))


library(FactoMineR)                



mc1<-MCA(kdf, ncp = 2, ind.sup = NULL, quanti.sup =NULL,
         quali.sup = NULL, excl=NULL, graph = TRUE,
         level.ventil = 0.5, axes = c(1,2), row.w = NULL,
         method="Burt", na.method="NA", tab.disj=NULL)

plot(mc1,invisible=c("ind","quali.sup"),hab="quali")


dimdesc(mc1,axes=1:2)



```