---
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")
source("workflow1_biomasse.r")

```
## Biomasse
```{r ,fig=TRUE}
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}

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

```

```{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 and vegetation types
```{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

### PCA all abiotic variables
```{r ,fig=TRUE}
releves_all = readRDS("releves_all_plus.rds")
releves_all$date600 = releves_all$last_snow_day.gdd300 + releves_all$gddspeed.gdd600
str(releves_all)

library(ade4)
vars = c("alti", "slope", "lf", "northing", "easting", "snowfree_time", "cumgdd0", "cumgdd60d", "albedo.gdd300",
"albedo.gdd600", "albedo.gdd900","temperature.gdd300", "frost.gdd300",
"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

### climate et type vegetation
```{r ,fig=TRUE}
ggplot(releves_all, aes(x = reorder(ref_typoveg, cumgdd0, na.rm=TRUE), y = cumgdd0)) +
  geom_boxplot() +
  labs(y="cum GDD 0degC", x="Type vegetation")
```
```{r ,fig=TRUE}
ggplot(releves_all, aes(x = reorder(ref_typoveg, cumgdd60d, na.rm=TRUE), y = cumgdd60d)) +
  geom_boxplot() +
  labs(y="cum GDD 60d", x="Type vegetation")
```

```{r ,fig=TRUE}
ggplot(releves_all, aes(x = reorder(ref_typoveg, temperature.gdd300, na.rm=TRUE), y = temperature.gdd300)) +
  geom_boxplot() +
  labs(y="temperature moy", x="Type vegetation")
```

```{r ,fig=TRUE}
  ggplot(releves_all, aes(x = reorder(ref_typoveg, last_snow_day.gdd300, na.rm=TRUE), y = last_snow_day.gdd300)) +
    geom_boxplot() +
    labs(y="LSD", x="Type vegetation")
```
```{r ,fig=TRUE}
ggplot(releves_all, aes(x = reorder(ref_typoveg, gddspeed.gdd600, na.rm=TRUE), y = gddspeed.gdd600)) +
    geom_boxplot() +
    labs(y="gddspeed 600", x="Type vegetation")
```

```{r ,fig=TRUE}
releves_all$date600 = releves_all$last_snow_day.gdd300 + releves_all$gddspeed.gdd600

ggplot(releves_all, aes(x = reorder(ref_typoveg, date600, na.rm=TRUE), y = date600)) +
    geom_boxplot() +
    labs(y="dateGDD600", x="Type vegetation")
```




### Modèle exploratoire complet
```{r ,fig=TRUE}
library(randomForest)
vars = c("alti", "slope", "lf", "northing", "easting", "snowfree_time", "cumgdd0", "cumgdd60d", "albedo.gdd300",
"albedo.gdd600", "albedo.gdd900","temperature.gdd300", "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", "ref_typoveg" )
dat = na.omit(releves_all[-which(releves_all$ref_typoveg %in% c("HUM", "LAND", "NA", "EBOU")),c(vars, "hmean")])
dat$lf = as.factor(dat$lf)
dim(dat)
mod = randomForest(hmean ~ ., data = dat, ntree = 100, mtry=16)

mod
varImpPlot(mod)
```

### Response graphs
```{r ,fig=TRUE}
library(plotmo)
# plotmo(mod, ylim = NA)
```

### Graphs univariés
```{r ,fig=TRUE}
ggplot(releves_all, aes(x = cumgdd0, y = hmean)) +
  geom_point() +
  geom_smooth(method="auto", se=TRUE, fullrange=FALSE, level=0.95) +
  labs(y="Hauteur", x="GDD 0degree")
```

```{r ,fig=TRUE}
ggplot(releves_all, aes(x = slope, y = hmean)) +
  geom_point() +
  geom_smooth(method="auto", se=TRUE, fullrange=FALSE, level=0.95) +
  labs(y="Hauteur", x="slope")
```

```{r ,fig=TRUE}
ggplot(releves_all, aes(x = gddspeed.gdd600, y = hmean)) +
  geom_point() +
  geom_smooth(method="auto", se=TRUE, fullrange=FALSE, level=0.95) +
  labs(y="Hauteur", x="GDD.600")
```


### Modèle par type

```{r ,fig=TRUE}
vars = c("alti", "slope", "lf", "northing", "easting", "snowfree_time", "cumgdd0", "cumgdd60d", "albedo.gdd300",
"albedo.gdd600", "albedo.gdd900","temperature.gdd300", "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" )
#vars = c("slope", "lf", "northing", "easting", "snowfree_time", "cumgdd60d", "ref_typoveg")

# dat = na.omit(releves_all[-which(releves_all$ref_typoveg %in% c("HUM", "LAND", "NA", "EBOU")),c(vars, "hmean")])
dat = na.omit(releves_all[which(releves_all$ref_typoveg =="PROD"),c(vars, "hmean")])
# dat$lf = as.factor(dat$lf)
dim(dat)
mod = randomForest(hmean ~., data = dat, ntree = 100, mtry=16)

mod
varImpPlot(mod)
```

```{r ,fig=TRUE}
vars = c("alti", "slope", "lf", "northing", "easting", "snowfree_time", "cumgdd0", "cumgdd60d", "albedo.gdd300",
"albedo.gdd600", "albedo.gdd900","temperature.gdd300", "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" )
#vars = c("slope", "lf", "northing", "easting", "snowfree_time", "cumgdd60d", "ref_typoveg")

# dat = na.omit(releves_all[-which(releves_all$ref_typoveg %in% c("HUM", "LAND", "NA", "EBOU")),c(vars, "hmean")])
dat = na.omit(releves_all[which(releves_all$ref_typoveg =="ALP"),c(vars, "hmean")])
# dat$lf = as.factor(dat$lf)
dim(dat)
mod = randomForest(hmean ~., data = dat, ntree = 100, mtry=16)

mod
varImpPlot(mod)
```

```{r ,fig=TRUE}
vars = c("alti", "slope", "lf", "northing", "easting", "snowfree_time", "cumgdd0", "cumgdd60d", "albedo.gdd300",
"albedo.gdd600", "albedo.gdd900","temperature.gdd300", "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" )
#vars = c("slope", "lf", "northing", "easting", "snowfree_time", "cumgdd60d", "ref_typoveg")

# dat = na.omit(releves_all[-which(releves_all$ref_typoveg %in% c("HUM", "LAND", "NA", "EBOU")),c(vars, "hmean")])
dat = na.omit(releves_all[which(releves_all$ref_typoveg =="SUB"),c(vars, "hmean")])
# dat$lf = as.factor(dat$lf)
dim(dat)
mod = randomForest(hmean ~., data = dat, ntree = 100, mtry=16)

mod
varImpPlot(mod)
```

```{r ,fig=TRUE}
vars = c("alti", "slope", "lf", "northing", "easting", "snowfree_time", "cumgdd0", "cumgdd60d", "albedo.gdd300",
"albedo.gdd600", "albedo.gdd900","temperature.gdd300", "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" )
#vars = c("slope", "lf", "northing", "easting", "snowfree_time", "cumgdd60d", "ref_typoveg")

# dat = na.omit(releves_all[-which(releves_all$ref_typoveg %in% c("HUM", "LAND", "NA", "EBOU")),c(vars, "hmean")])
dat = na.omit(releves_all[which(releves_all$ref_typoveg =="NIV"),c(vars, "hmean")])
# dat$lf = as.factor(dat$lf)
dim(dat)
mod = randomForest(hmean ~., data = dat, ntree = 100, mtry=16)

mod
varImpPlot(mod)
```

### Interannuel

```{r ,fig=TRUE}
vars = c("alti", "slope", "lf", "northing", "easting", "snowfree_time", "cumgdd0", "cumgdd60d", "albedo.gdd300",
"albedo.gdd600", "albedo.gdd900","temperature.gdd300", "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")
#vars = c("slope", "lf", "northing", "easting", "snowfree_time", "cumgdd60d", "ref_typoveg")

# dat = na.omit(releves_all[-which(releves_all$ref_typoveg %in% c("HUM", "LAND", "NA", "EBOU")),c(vars, "hmean")])
dat = na.omit(releves_all[-which(releves_all$ref_typoveg %in% c("HUM", "LAND", "NA", "EBOU")),c(vars, "hmean")])

# dat$lf = as.factor(dat$lf)
dim(dat)
dat
mod = randomForest(hmean ~., data = dat, ntree = 100, mtry=16)

mod
varImpPlot(mod)
```