MAJ analyse CBNA

.gitignore

@@ -2,3 +2,7 @@
 *.xlsx
 *.csv
 CEPAZ_Zone_Etude_et_ZP/*
+Extract_raster_zp(20_02_2020)/*
+*.shp
+*.tif
+*.rds

diversityEcosystem.R

+# load data
+library(openxlsx)
+
+lf.tab = read.xlsx("LANDFORM/Export_ZP_landform.xlsx")
+head(lf.tab)
+colnames(lf.tab)[2] = "prop_NA"
+
+library(rgdal)
+zp = readOGR("Extract_raster_zp(20_02_2020)", "CEPAZ_ZP")
+head(zp@data)
+
+library(raster)
+expoeo = raster("Extract_raster_zp(20_02_2020)/expoeo.tif")
+projection(expoeo)
+expons = raster("Extract_raster_zp(20_02_2020)/expons.tif")
+projection(expons)
+tpi = raster("Extract_raster_zp(20_02_2020)/extrtpi.tif")
+projection(tpi)
+pente = raster("Extract_raster_zp(20_02_2020)/pentezp2.tif")
+projection(pente)
+
+# ray = raster("Extract_raster_zp(20_02_2020)/rayextr.tif")
+# projection(ray) == projection(pente)
+# ray2 = resample(ray, pente, method = "ngb")
+# saveRDS(ray2, "ray2.rds")
+ray2 = readRDS("ray2.rds")
+ray2
+
+lf = raster("LANDFORM/landformsTPI_saga.tif")
+projection(lf)
+lf_proj = projectRaster(lf, pente, method = "ngb")
+lf.clip = crop(lf_proj, extent(pente))
+# saveRDS(lf.clip, file = "lf_clip.rds")
+# lf.clip = readRDS("lf_clip.rds")
+lf.clip
+projection(zp) == projection(expoeo)
+
+library(sp)
+zp_proj = spTransform(zp, crs(projection(pente)))
+
+#=========================
+
+zp_x = zp_proj[1,]
+expoeo.vec = extract(expoeo, zp_x)[[1]]
+expons.vec = extract(expons, zp_x)[[1]]
+
+projection(lf.clip) == projection(zp_x)
+lf.vec = extract(lf.clip, zp_x)[[1]]
+
+
+tab = data.frame( expoeo = expoeo.vec, 
+                  expons = expons.vec,
+                  lf = ,
+                  ray = extract(ray2, zp_x)[[1]],
+                  pente = extract(pente, zp_x)[[1]],
+                  tpi = extract(tpi, zp_x)[[1]] )
+tab$lf = as.factor(tab$lf)
+
+library(FD)
+Dij = gowdis(tab)
+
+res_x = mean(as.vector(Dij))
+                  
+

exploration.Rmd

 ---
-title: "EXPLORATION HABITATS CEPAZ"
+title: "ANALYSE CBNA CEPAZ"
 author: "Isabelle Boulangeat"
-date: "22/01/2020"
-output: 
+date: "27/07/2020"
+output:
   md_document:
      variant: markdown_github
-editor_options: 
+editor_options:
   chunk_output_type: console
 ---
 
@@ -15,131 +15,172 @@ knitr::opts_chunk$set(echo = TRUE, warning=FALSE)
 
 ## Reading data
 
+Tous les relevés phyto du CBNA après 2013 (incluant les zones CEPAZ et les autres) sont chargés.
+
+La table des Zones Pastorales avec infos conplémentaires est chargée.
+
+La liste de taxons du CBNA relevés après 2013 avec le statut des espèces, et la liste des espèces de la zp avec leur statuts et autres stats.
+
 ```{r load_data}
-library("openxlsx")
-datphyto = read.xlsx("releves_phyto_CEPAZ_2013+.xlsx")
-head(datphyto,3)  
+library(openxlsx)
+# datphyto = read.xlsx("releves_phyto_CEPAZ_2013+.xlsx")
+# head(datphyto,3)
 
 library(rgdal)
-polyZone = readOGR("CEPAZ_Zone_Etude_et_ZP", "CEPAZ_Zones_Pastorales")
-head(polyZone@data,3)
-nrow(polyZone@data)
-
+datphyto = readOGR(dsn = "releves_phyto_cbna_sup2013_v3.shp", layer = "releves_phyto_cbna_sup2013_v3")
+# nb 355 035 releves-especes
+head(datphyto@data,3)
+colnames(datphyto@data)
+head(datphyto@data$CODE,100) # code ZP ou NA
+
+insideZP = datphyto@data[which(!is.na(datphyto$CODE)),]
+dim(insideZP)
+
+# polyZone = readOGR("CEPAZ_Zone_Etude_et_ZP", "CEPAZ_Zones_Pastorales")
+# head(polyZone@data,3)
+# nrow(polyZone@data[which(polyZone$INSEEDEP %in% c("01", "74", "73", "38", "26", "05", "04")),])
+#
 tabZP = read.csv2("ZP_table.csv", encoding = "UTF-8")
 head(tabZP,3)
 colnames(tabZP)[1] = "CODE"
 nrow(tabZP)
 
-status = read.xlsx("stat_taxons_zp_statuts.xlsx")
-head(status,3)
-nrow(status)
+statuts_all = read.xlsx("taxon_cbna_2013+.xlsx")
+head(statuts_all,3)
+nrow(statuts_all)
+statuts_all$UICN = statuts_all$UICN_PACA
+statuts_all$UICN[which(is.na(statuts_all$UICN_PACA))] = statuts_all$UICN_RA[which(is.na(statuts_all$UICN_PACA))]
+
+sp_zp = read.xlsx("stat_taxons_zp_statuts.xlsx")
+head(sp_zp)
 ```
 
-## Taille des ZP repr??sent??es
+## Taille des ZP représentées
 ```{r ZP_size, fig=TRUE}
 library(dplyr)
 library(ggplot2)
 
-# taille des ZP repr??sent??es
-tabZP_phyto = merge(tabZP, datphyto, by = "CODE", all=FALSE)
-tabZP_phyto$dataset = "phyto" 
-tabZP_points = merge(tabZP, status, by = "CODE", all=FALSE)
-tabZP_points$dataset = "points" 
+# taille des ZP représentées
+tabZP_phyto = merge(tabZP, insideZP, by = "CODE", all=FALSE)
+tabZP_phyto$dataset = "phyto"
 
 tabZP$dataset = "all"
 selectCol = c("CODE", "dataset", "SURFACE")
-rbind(tabZP_phyto[, selectCol], tabZP_points[, selectCol], tabZP[, selectCol]) %>%
+rbind(tabZP_phyto[, selectCol], tabZP[, selectCol]) %>%
   ggplot(aes(SURFACE, fill=dataset, colour = dataset, stat(density))) +
   geom_histogram(alpha=0.2, position = "identity")
 ```
 
-## Simple stats
-```{r simple_statistiques}
-# nb ZP o?? il y a (au moins) un releve phyto 
-length(unique(datphyto$CODE))
+On a une sous-représentation des ZP de petite taille.
+
+## Taille des ZP par departement
+```{r ZP_size_dep, fig=TRUE}
+selectCol = c("CODE", "dep", "SURFACE")
+tabZP_phyto[, selectCol] %>%
+  ggplot(aes(SURFACE, fill=dep, colour = dep, stat(density))) +
+  geom_histogram(alpha=0.2, position = "identity")
+```
+
+## Taille des ZP par milieu
+```{r ZP_size_milieu, fig=TRUE}
+selectCol = c("CODE", "MILIEU", "SURFACE")
+tabZP_phyto[, selectCol] %>%
+  ggplot(aes(SURFACE, fill=MILIEU, colour = MILIEU, stat(density))) +
+  geom_histogram(alpha=0.2, position = "identity")
+```
+
+
+
+## Nb de relevés par ZP
+```{r nbrel}
+# nb ZP où il y a (au moins) un releve phyto
+length(unique(insideZP$CODE))
+length(unique(tabZP$CODE))
 
 # nombre de releve phyto par ZP
-zp_releve = datphyto %>% group_by(CODE) %>% summarise(n = length(unique(numchrono)))
+zp_releve = insideZP %>% group_by(CODE) %>% summarise(n = length(unique(numchrono)))
 summary(zp_releve$n)
 
-# nb de ZP o?? il y a (au moins) une observation d'espece
-length(unique(status$CODE))
-
-# ann??es 1ere et derniere observations
-summary(status$`_min`)
-summary(status$`_max`)
 
-# couper ?? 2013 (pour avoir la m??me donn??e que les relev??s phyto)
-status2013 = status[which(status$`_max`>=2013),]
-length(unique(status2013$CODE))
 ```
 
 
-# Diversit?? alpha des ZP
+# Diversité alpha des ZP
 ```{r diversite_alpha, fig=TRUE}
-# nombre d'especes par ZP, par observations points
-zp_espece = status2013 %>% group_by(CODE) %>% summarise(sp_richness = length(unique(numtaxon)))
-summary(zp_espece$sp_richness)
-
-# nombre d'especes par ZP, par releves phyto 
-zp_espece_phyto = datphyto %>% group_by(CODE) %>% summarise(sp_richness_phyto = length(unique(numtaxon)))
-summary(zp_espece_phyto$sp_richness_phyto)
-
-# comparaison des 2 jeux de donn??es
-zp_espece.all = merge(zp_espece, zp_espece_phyto, all.x=TRUE, all.y = TRUE, by = "CODE")
-head(zp_espece.all)
-nrow(zp_espece.all)==length(unique(status2013$CODE))
-
-zp_espece.all %>% 
-  ggplot(aes(sp_richness, sp_richness_phyto)) +
-  geom_point(shape=19) +
-  geom_smooth() +
-  geom_abline(intercept = 0, slope=1)
-```
 
-Le jeu de donn??es phytosociologique sous-estime la diversit?? par rapport au jeu de donn??es d'observations individuelles, surtout pour les niveaux de richesse > 100 esp??ces.
+# nombre d'especes par ZP, pour releves phyto
+zp_espece_phyto = insideZP %>% group_by(numchrono) %>% summarise(sp_richness_phyto = length(unique(numtaxon)))
 
-## Analyse par d??partement et milieu
-```{r diversity_groups, fig=TRUE}
-zp_espece$dataset = "points"
-zp_espece_phyto$dataset = "phyto"
-names(zp_espece_phyto)[2] = "sp_richness"
-zp_espece.all2 = rbind(zp_espece, zp_espece_phyto)
-dat_div = merge(tabZP[,-which(colnames(tabZP)=="dataset")], zp_espece.all2, by="CODE", all.x=TRUE)
-head(dat_div)
+div_alpha_zp = unique(merge(insideZP[, c("CODE", "numchrono")], zp_espece_phyto, by = "numchrono"))
+head(div_alpha_zp)
+summary(div_alpha_zp)
 
-dat_div %>% 
-  ggplot(aes(MILIEU, sp_richness, fill = dataset)) +
-  geom_boxplot(notch=TRUE) + 
-  stat_boxplot(na.rm=TRUE) +
-  ylim(0,100) 
+stat_alpha = div_alpha_zp %>% group_by(CODE) %>% summarise(alpha_mean = mean(sp_richness_phyto), nb_rel = length(unique(numchrono)))
+summary(stat_alpha)
+
+# richesse des releves phyto hors ZP (ou total)
+div_alpha_tot = datphyto@data %>% group_by(numchrono) %>% summarise(sp_richness_phyto = length(unique(numtaxon)))
+summary(div_alpha_tot)
+
+# detail par milieu et departement (figure)
+stat_alpha_details = unique(merge(stat_alpha, tabZP, by = "CODE"))
+head(stat_alpha_details)
 
-dat_div %>% group_by(MILIEU) %>%
-  ggplot(aes(MILIEU, sp_richness, fill = dataset)) +
-  geom_boxplot(notch=TRUE) + 
+stat_alpha_details %>% group_by(MILIEU) %>%
+  ggplot(aes(MILIEU, alpha_mean)) +
+  geom_boxplot(notch=TRUE) +
   stat_boxplot(na.rm=TRUE) +
   ylim(0,100) +
   facet_wrap(~INSEEDEP)
 
+
 ```
 
-# Statut des esp??ces des ZP
-Nombre et proportion d'esp??ces liste rouge parmi les esp??ces relev??es
-Rappel:  Eteinte (EX), Eteinte ?? l?????tat sauvage (EW), En danger critique (CR), En danger (EN), Vuln??rable (VU), Quasi menac??e (NT), Pr??occupation mineure (LC), Donn??es insuffisantes (DD), Non ??valu??e (NE).
 
-```{r status_especes}
-sp_phyto = unique(datphyto$numtaxon)
-length(sp_phyto)
-sp_points = unique(status2013$numtaxon)
-length(sp_points)
-length(unique(status$numtaxon))
+# Statut des espèces des ZP
+Nombre et proportion d'espèces liste rouge parmi les espèces relevées
+Rappel:  Eteinte (EX), Eteinte à l'état sauvage (EW), En danger critique (CR), En danger (EN), Vulnérable (VU), Quasi menacée (NT), Préoccupation mineure (LC), Données insuffisantes (DD), Non évaluée (NE).
+
+Evalué avec les données "points" comprenant les relevés phyto et des relevés supplémentaires d'espèces cible.
+
+```{r status_especes, fig=TRUE}
+sp_phyto = unique(datphyto@data$numtaxon)
+length(sp_phyto) # 3789 especes
+sp_phyto_zp = unique(insideZP$numtaxon)
+length(sp_phyto_zp) # 1587 especes
+
+sp_points = unique(statuts_all$numtaxon) # après 2013
+length(sp_points) # 5606
+sp_zp_2013 = sp_zp[which(sp_zp$`_max`>=2013),]
+sp_points_zp = unique(sp_zp_2013$numtaxon)
+length(sp_points_zp) # 2146
+
+
+statuts_all$UICN =as.factor(statuts_all$UICN)
+levels(statuts_all$UICN) = c("CR","CR","DD","EN","EW","LC", "NE"  ,"NT",  "RE",  "VU")
+aa = table(unique(statuts_all[, c("numtaxon", "UICN")])$UICN) / length(sp_points)
+sum(is.na(unique(statuts_all[, c("numtaxon", "UICN")])$UICN))/ length(sp_points)
+
+
+sp_zp_2013$UICN =as.factor(sp_zp_2013$lr_fr)
+
+zp = table(unique(sp_zp_2013[, c("numtaxon", "UICN")])$UICN) / length(sp_points_zp)
+sum(is.na(unique(sp_zp_2013[, c("numtaxon", "UICN")])$UICN))/ length(sp_points_zp)
+zp_vec = unclass(zp)
+
+library(reshape2)
+rbind(unclass(aa),c(CR=0, zp_vec[1:2], EW=0, zp_vec[3], NE=0, zp_vec[4], RE=0, zp_vec[5])) %>% melt() %>%
+  ggplot(aes(x = Var2, y = value, fill = factor(Var1))) +
+  geom_bar(stat = "identity", position = position_dodge2()) +
+  scale_x_discrete(name="statut UICN") +
+  scale_y_continuous(name = "proportion de taxons") +
+  labs(fill = "emprise") +
+  scale_fill_hue(labels = c("CBNA", "ZP"))
 
-table(unique(status2013[, c("numtaxon", "lr_fr")])$lr_fr)
 
-status_phyto= status[which(status$numtaxon %in% datphyto$numtaxon),]
-length(unique(status_phyto$numtaxon)) == length(unique(datphyto$numtaxon))
-# attention 3 especes n'ont pas de statut attribu??
-table(unique(status_phyto[, c("numtaxon", "lr_fr")])$lr_fr)
 ```
-On voit bien que les observations ponctuelles r??pertorient plus d'esp??ces vuln??rables.
+On retrouve 38% de la flore de la zone d'agrément du CBNA échantillonnée dans les ZP.
+
+Attention pour les statuts il manque 36% de données dans la liste complète et 22% dans les espèces des ZP.
 
+Attention, toute cette évaluation exclue les départements hors zone CBNA (au Sud).