---
title: "HMM plots for paper"
author: "Rachel Mawer"
date: "2023-10-05"
output: 
  html_document:
    toc: true
    toc_float: true
    theme: cerulean
---

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

#main packages #delete/add as needed
library(dplyr)
library(terra)
#library(momentuHMM)
library(ggplot2)
library(knitr)
library(ggpubr)
library(plotly)
library(pracma) #for standard error
library(reactable)
library(wesanderson)
library(amt)


theme_set(theme_bw())
theme_update()
```

## About

File to create plots relating to the HMM for the paper.

Plots will be made of:

* state distributions per species
* heatmaps of state occurence in river per species
* example track coloured by state
* histogram for % of steps per state for each fish
* step length turning angle distribution plots


```{r read in files}
barb_hmm <- readRDS("HMMs/barb_SI10.rds")
gr_hmm <- readRDS("HMMs/gr_SI10.rds")

dat_loc <-list.files("data/final ssf + hmm df",full.names = T) #location of datafiles

```

## Steps per state

```{r plots of steps per state}
#below df was previously made in the hmm+ssf explore stuff
step_percs <- read.csv("FINAL RESULTS/number steps per state per fish.csv")

st1 <- ggplot(step_percs,aes(x=perc_state_1))+
  geom_histogram(binwidth = 20)+
  facet_grid(species~.)+
  labs(x="Percentage of steps",y="Count",title="State 1")

st2 <- ggplot(step_percs,aes(x=perc_state_2))+
  geom_histogram(binwidth = 20)+
  facet_grid(species~.)+
  labs(x="Percentage of steps",y="",title="State 2")


ggarrange(plotlist=list(st1,st2),ncol=2,nrow = 1)

ggsave(filename="figures/hmm_steps_per_state.png",plot = st1,width=5,height=5)

```

Essentially, for barbel generally a 50-50 split, for grayling more 40-60. Some fish are only in one state.

Probably don't need to plot both states, since just inverses of each other

## State distributions


```{r plot state distributions,results='hide'}

cols <- wes_palette("Zissou1",n=5,type="discrete")[c(5,1)]

plot(barb_hmm,plotTracks = F,ask=F,col=cols)

ggsave(filename="figures/hmm_states_barbel.png",plot=last_plot(),width=6,height=6)
plot(gr_hmm,plotTracks = F,ask=F,col=cols)
ggsave(filename="figures/hmm_states_grayline.png",plot=last_plot(),width=6)



```



## Heat maps

Heat maps of state occurence in river per species.

```{r prep data for heatmaps}
for(i in 1:length(dat_loc)){
  dat <- read.csv(dat_loc[i])
  #remove false steps
  dat <- dat %>% filter(case_==TRUE)
  
  if(i==1){
    all_dat <- dat
  } else{
    all_dat <- rbind(all_dat,dat)
  }
}

#get river shapefile
library(rgdal,quietly=T)
river_shape <-readOGR('rivershp/new_river_shapefile.shp')
river_shape2 <- fortify(river_shape)
```

```{r make heatmaps}

#can facet grid by state and species

plt <- ggplot(all_dat,aes(x=x2_,y=y2_))+
  geom_polygon(data=river_shape2,aes(long,lat),col='black', fill=NA,linewidth=0.75)+
  stat_bin_2d(aes(fill=log(..density..)))+
  facet_grid(state~species)+
  scale_fill_continuous(type = "viridis")+
  scale_x_continuous(breaks=c(591600,591800))+
  scale_y_continuous(breaks=c(5296850,5296950))+
  geom_point(x=591683.2,y=5296867.5,col="red")+
  coord_fixed()+
  labs(x="x",y="y")+
  coord_cartesian(ylim=c(5296780,5297010),xlim=c(591500,591870),expand=F)#+
  #scale_fill_gradientn(colours=wes_palette("Zissou1",type="continuous"))

plt

ggsave(filename="figures/hmm_heatmaps_size_update.png",plot=plt,width=5,height=3.5)


```

Identifiable hotspots for state 2 by fish pass entrance (red spot), though for grayling its more in the area below - maybe resting?

## One fish track

using the all_dat from above, filter to 1 fish and 1 approach, then ggplot

Fish id 46852 track 1 looks nice.

```{r filter data and plot}

unique(all_dat$fish_id)
one_fish <- all_dat %>%filter(fish_id==46852) %>% arrange(t1_)
unique(one_fish$approach)

one_track <- one_fish %>% filter(approach==1)

plt <- 
ggplot(one_track,aes(x=x1_,y=y1_))+
  #geom_polygon(data=river_shape2,aes(long,lat),col='black', fill=NA)+
  geom_path(linewidth=0.75,aes(group=1,col=as.factor(state)))+
  geom_point(aes(col=as.factor(state)))+
  scale_color_manual(values=wes_palette("Zissou1",n=5,type="discrete")[c(5,1)])+
  labs(x="x",y="y",col="State")

plt

ggsave(filename="figures/hmm_example_track.png",plot=plt,width=5,height=3.5)


```

## Step length and turning angle distributions per state

Pooling all fish together, step length and turning angle distributions are plotted per state

```{r calculate distribution parameters,warning=F}
#read in all data for barbel and grayling, split into state

barb <- read.csv("data/steps with states/barbel_all.csv")
gray <- read.csv("data/steps with states/grayling_all.csv")


#rbind and do as one

all_d <- rbind(barb,gray)

sl_distrs <- all_d %>% group_by(species,State) %>% summarise(sl_shape = fit_distr(sl_,"gamma")[[2]][[1]],sl_scale=fit_distr(sl_,"gamma")[[2]][[2]], ta_kappa =fit_distr(ta_,"vonmises")$params$kappa,ta_mu =fit_distr(ta_,"vonmises")$params$mu )

reactable(sl_distrs)

```

Looking at the actual values for distributions, SL and TA distributions suggest that more directed movement in state 1 and with bigger step lengths, as expected for transit state.

Density plots per species/state:

Step lengths (with x axis limited to max 10):

```{r density plots for step lengths}

#generate 1000 estimates # of means? or 1000 steps

for(i in 1:nrow(sl_distrs)){
  row_i <- sl_distrs[i,]

  r_step <- rgamma(1000000,shape = row_i$sl_shape, 
                     scale = row_i$sl_scale) %>% as.data.frame()
  
  df <- data.frame(step_length = r_step,species=row_i$species,state=row_i$State)
  names(df)[1] <- "step_length"
  
  if(i==1){
    steps <- df
  } else{
    steps <- rbind(steps,df)
  }

    
    
}


#steps %>% group_by(mean_sl) %>% summarise(mean=mean(step_length))


s_plot <- ggplot(steps,aes(x=step_length,col=as.factor(state)))+
  geom_density()+
  scale_color_manual(values=wes_palette("Darjeeling1",n=5,type="discrete"),name="State")+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_grid(.~species)+
  labs(x="Step length (m)",y="Density")+
  coord_cartesian(xlim=c(0,10)) #limit x coords since go up to 200 :<

ggsave(filename="figures/step_length_dist.png",s_plot,width=6,height=2.5)


```


Turning angles:


```{r density plots for ta}


#generate 1000 estimates # of means? or 1000 steps

for(i in 1:nrow(sl_distrs)){
  row_i <- sl_distrs[i,]

  r_ang <-  circular::rvonmises(10000,mu = row_i$ta_mu, 
                     kappa = row_i$ta_kappa) %>% as.data.frame()
  
  df <- data.frame(ta = r_ang,species=row_i$species,state=row_i$State)
  names(df)[1] <- "turning_angle"
  
  if(i==1){
    ta <- df
  } else{
    ta <- rbind(ta,df)
  }

    
    
}

radtodeg <- function(rad) {((rad * 180) / (pi)) %% 360} #function since radians not limited to 0 and 2pi

#ta$deg <- ta$turning_angle *180/pi
ta$deg <- radtodeg(ta$turning_angle)

ggplot(ta,aes(x=turning_angle,col=as.factor(state)))+
  geom_density()+
  scale_color_manual(values=wes_palette("Darjeeling1",n=5,type="discrete"))+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_grid(species~.)#+


```

To make a nicer plot and centered on 0 (e.g. directed movement), for values >pi, I will subtract 2pi

```{r ta plot 2}

ta$ta_2 <- case_when(ta$turning_angle>pi ~ta$turning_angle-2*pi,
                     ta$turning_angle<=pi ~ ta$turning_angle)

ta_plot <- ggplot(ta,aes(x=ta_2,col=as.factor(state)))+
  geom_density()+
  scale_color_manual(values=wes_palette("Darjeeling1",n=5,type="discrete"),name="State")+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_grid(species~.)+
  labs(y="Density",x="Radians")

ggsave(filename="figures/ta_plot.png",ta_plot,width=6,height=4)
```

Can see that state 1 will have more directed movement.