---
title: "README"
author: "Isabelle Boulangeat"
date: "25/01/2021"
output:
  html_document:
    keep_md: yes
    variant: markdown_github
editor_options:
  chunk_output_type: console
always_allow_html: true
---

```{r setup, include=FALSE}
library(knitr)
# library(kableExtra)
knitr::opts_chunk$set(echo = TRUE)
# library(Jmisc)
library(tidyr)
sapply(list.files("R_fct"), function(x)source(paste0("R_fct/",x)))
data = readRDS("data.rds")

data$otable[which(is.na(data$otable$hauteur)),]

dat_h_veg = unique(merge(data$h, data$sites[,c("id_site", "ref_typoveg")], by.x = "ref_site", by.y ="id_site"))

length(unique(dat_h_veg$releve))
dim(dat_h_veg)

summary(dat_h_veg$hmean)

#==============================

library(dplyr)
# library(reshape)
library(tidyr)
library(ggplot2)

```
## Biomasse
```{r ,fig=TRUE}

library(dplyr)
library(tidyr)
library(ggplot2)
dat_long <- dat_h_veg %>%
 select(hmean, hmad, ref_typoveg) %>%
 drop_na() %>%
 group_by(ref_typoveg) %>%
 summarize(hauteur=mean(hmean), var_inter=mad(hmean), var_intra = mean(hmad)) %>%
 gather(stat, value, hauteur:var_intra)

ggplot(dat_long, aes(fill=stat, x=ref_typoveg, y=value)) +
  geom_bar(position = "dodge", stat = "identity") +
  facet_wrap(~stat, ncol = 1, scales = "free")

```

```{r ,fig=TRUE}

## boxplot par milieu ##

ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmean, na.rm=TRUE), y = hmean)) +
  geom_boxplot() +
  labs(y="Hauteur moyenne (est. biomasse)", x="Type vegetation")

```

```{r ,fig=TRUE}

## H series par milieu ##

pl <- ggplot(dat_h_veg, aes(x = date_releve, y = hmean)) +
  geom_line(aes(color = ref_site), show.legend = FALSE) +
  facet_wrap(~ref_typoveg)

pl + theme(legend.position = "none")


```

```{r ,fig=TRUE}
calage = read.csv("../PNE_calage_reformat.csv", sep= ";", dec = ",")
head(calage)
str(calage)
calage$releve = paste0(calage$site, calage$Releveur)

cal_mean = unique(calage %>% group_by(releve) %>% summarize(site = site, hmean = mean(hauteur, na.rm=TRUE)))

ggplot(cal_mean, aes(x = reorder(site, hmean), y = hmean)) +
  geom_boxplot() +
  geom_point() 

obs_error = cal_mean %>% group_by(site) %>% summarize(hmad = mad(hmean), hmad_st = mad(hmean)/mean(hmean), hm = mean(hmean))


```

```{r ,fig=TRUE}
ggplot(obs_error, aes(x = reorder(site, hm), y = hmad_st)) +
  geom_bar(stat = "identity")

```

```{r ,fig=TRUE}
### sd vs mad
ggplot(dat_h_veg, aes(x = hmad, y = hsd)) +
  geom_point() 

```

```{r ,fig=TRUE}
### dans le temps
ggplot(dat_h_veg, aes(x = date_releve, y = hmad)) +
  geom_line(aes(color = ref_site), show.legend = FALSE) +
  facet_wrap(~ref_typoveg)

## boxplot par milieu ##

```

```{r ,fig=TRUE}
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmad, na.rm=TRUE), y = hmad)) +
  geom_boxplot() +
  labs(y="Variation intra de hauteur", x="Type vegetation")

```

```{r ,fig=TRUE}
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmad/hmean, na.rm=TRUE), y = hmad/hmean)) +
  geom_boxplot() +
  labs(y="Variation intra de hauteur (stand. moy.)", x="Type vegetation")

```

## Sites topography
```{r ,echo=FALSE}
sites_all = readRDS("sites_all_plus.rds")
```

```{r ,fig=TRUE}
ggplot(sites_all@data, aes(x = reorder(ref_typoveg, alti, na.rm=TRUE), y = alti)) +
  geom_boxplot() +
  labs(y="alti", x="Type vegetation")
```
```{r ,fig=TRUE}
ggplot(sites_all@data, aes(x = reorder(ref_typoveg, slope, na.rm=TRUE), y = slope)) +
  geom_boxplot() +
  labs(y="slope", x="Type vegetation")
```
```{r ,fig=TRUE}
ggplot(sites_all@data, aes(x = reorder(ref_typoveg, northing, na.rm=TRUE), y = northing)) +
  geom_boxplot() +
  labs(y="northing", x="Type vegetation")
```
```{r ,fig=TRUE}
ggplot(sites_all@data, aes(x = reorder(ref_typoveg, easting, na.rm=TRUE), y = easting)) +
  geom_boxplot() +
  labs(y="easting", x="Type vegetation")
```

### Ensemble de variables abiotiques
```{r ,fig=TRUE}
releves_all = readRDS("releves_all_plus.rds")

library(ade4)
vars = c("alti", "slope", "lf", "northing", "easting", "snowfree_time", "cumgdd0", "cumgdd60d", "albedo.gdd300",
"albedo.gdd600", "albedo.gdd900", "frost_severe.gdd300", "frost_severe.gdd600",
"frost_severe.gdd900", "frost.gdd300", "frost.gdd600", "frost.gdd900",
"gddspeed.gdd300", "gddspeed.gdd600", "gddspeed.gdd900", "radiations.gdd300",
"radiations.gdd600", "radiations.gdd900", "rainfall.gdd300", "rainfall.gdd600",
"rainfall.gdd900", "snowdays.gdd300", "snowdays.gdd600", "snowdays.gdd900" )
dat = na.omit(releves_all[,vars])
dim(dat)
pca = dudi.pca(dat, scan=FALSE, nf = 4)
inertia.dudi(pca, row=FALSE, col =TRUE)
par(mfrow=c(1,2))
s.corcircle(pca$co)
s.corcircle(pca$co, xax=3, yax=4)
```

### Modélisation de la hauteur moyenne
```{r ,fig=TRUE}
dat = na.omit(releves_all[,c(vars, "hmean")])
dat$lf = as.factor(dat$lf)
dim(dat)
library(randomForest)
mod = randomForest(hmean ~ gddspeed.gdd300 + radiations.gdd300 + cumgdd60d + rainfall.gdd300 + rainfall.gdd900 + frost_severe.gdd300 + slope + northing + easting + lf, data = dat, ntree = 150, mtry=16)

mod
varImpPlot(mod)
```

```{r ,fig=TRUE}

mod = randomForest(hmean ~ gddspeed.gdd300 + radiations.gdd300 + cumgdd60d + rainfall.gdd900 + slope + northing + easting + lf, data = dat, ntree = 100, mtry=16)

mod

library(plotmo)
plotmo(mod, ylim = NA)

```