Commit 17ce3da7 authored by Boulangeat Isabelle's avatar Boulangeat Isabelle
Browse files

explo

parent c390ea65
......@@ -24,16 +24,7 @@ data = readRDS("data.rds")
```{r ,fig=TRUE}
dat_h_veg = merge(data$h, data$sites[,c("id_site", "ref_typoveg")], by.x = "ref_site", by.y ="id_site")
str(dat_h_veg)
summary(dat_h_veg$hmean)
#==============================
## H variation par milieu ##
#===============================
library(dplyr)
# library(reshape)
library(tidyr)
library(ggplot2)
dat_long <- dat_h_veg %>%
......
......@@ -37,7 +37,12 @@ import_from_db <- function(){
sql_biomass_h <- "SELECT * FROM sentinelle.biomasse_occurences"
data_biomass_h <- dbGetQuery(con, sql_biomass_h)
data_h = data_biomass_h %>% group_by(ref_releve) %>% summarize(hmean=mean(hauteur), hsd=sd(hauteur), hmad = mad(hauteur))
data_biomass_h$releve =
unlist(lapply(data_biomass_h$ref_releve, function(x){
paste(strsplit(x, "_")[[1]][-2], collapse="_")
}))
data_h = data_biomass_h %>% group_by(releve) %>% summarize(hmean=mean(hauteur), hsd=sd(hauteur), hmad = mad(hauteur))
# str(data_h)
# tail(data_h)
......@@ -45,11 +50,15 @@ import_from_db <- function(){
data_biomass_rel <- dbGetQuery(con, sql_biomass_rel)
# str(data_biomass_rel)
data_biomass_rel$releve =
unlist(lapply(data_biomass_rel$id_releve, function(x){
paste(strsplit(x, "_")[[1]][-2], collapse="_")
}))
data_h_all = merge(data_h, data_biomass_rel[,c("id_releve", "date_releve", "ref_observateur")], by.x="ref_releve", by.y = "id_releve")
data_h_all = merge(data_h, data_biomass_rel[,c("releve", "date_releve", "ref_observateur")], by="releve")
# str(data_h_all)
data_h_all$ref_site = unlist(lapply(data_h_all$ref_releve, function(x)strsplit(x, "_")[[1]][1]))
data_h_all$ref_site = unlist(lapply(data_h_all$releve, function(x)strsplit(x, "_")[[1]][1]))
dbGetQuery(con, "DROP TABLE chabli.biomasse_h")
......@@ -70,5 +79,5 @@ import_from_db <- function(){
saveRDS(data_h_all, file="data_h.rds")
saveRDS(data_biomass_site, file="data_sites.rds")
return(list(h = data_h_all, sites = data_biomass_site))
return(list(h = data_h_all, sites = data_biomass_site, otable = data_biomass_h))
}
No preview for this file type
No preview for this file type
library(tidync)
library(ncmeta)
library(tidyr)
library(raster)
library(sp)
#=============================================
#--- Data input
#=============================================
sites <- readRDS("data_sites.rds
sites <- sites[which(!is.na(sites$x)),]
coordinates(sites)= cbind(sites$y, sites$x)
proj4string(sites) = CRS("+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
NoP_raster = raster()
plot(NoP_raster)
points(sites)
spTransform()
extract()
# path_data_allslopes = "/Volumes/infogeo/Meteo_France/SAFRAN_montagne-CROCUS_2019/alp_allslopes"
path_data_allslopes = "/Volumes/ISA-RESEARCH/alp_allslopes"
safran_classes <- read.table("safran_classes.txt", h=T, sep="\t")
#=============================================
#--- Extract function
#=============================================
# list of variables:
#
extract_meteo <- function(path_data_allslopes, years, variables, NoP){
# library required
require(tidync)
require(ncmeta)
require(tidyr)
# check filepath dataset
if(!(length(list.files(path_data_allslopes)) > 0) ) stop("not able to connect to data path or empty path")
# temporal loop
for (y in years){
# read data and stock for a year y
pro_file = paste0(path_data_allslopes, "/reanalysis/pro/PRO_",as.character(y),"080106_",as.character(y+1),"080106.nc")
src <- tidync(pro_file)
data_y <- src %>% activate("D0,D3") %>% hyper_filter(Number_of_points = Number_of_points %in% NoP) %>% hyper_tibble(select_var = variables)
# TODO: compute values from data_y
}
# build output
# TODO!!
return(output)
}
......@@ -8,26 +8,34 @@ config.db.siddt()
data=import_from_db()
saveRDS(data, "data.rds")
data = readRDS("data.rds")
#==============================
## dataframe ##
#===============================
head(data$sites)
head(data$h)
head(data$otable)
dat_h_veg = merge(data$h, data$sites[,c("id_site", "ref_typoveg")], by.x = "ref_site", by.y ="id_site")
data$otable[which(is.na(data$otable$hauteur)),]
str(dat_h_veg)
summary(dat_h_veg$hmean)
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)
#==============================
## H variation par milieu ##
#===============================
library(dplyr)
# library(reshape)
library(tidyr)
library(ggplot2)
#==============================
## variation par milieu des mesures biomasse ##
#===============================
dat_long <- dat_h_veg %>%
select(hmean, hmad, ref_typoveg) %>%
drop_na() %>%
......@@ -39,16 +47,63 @@ ggplot(dat_long, aes(fill=stat, x=ref_typoveg, y=value)) +
geom_bar(position = "dodge", stat = "identity") +
facet_wrap(~stat, ncol = 1, scales = "free")
## 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")
## 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")
#==============================
## H boxplots ##
## effet obs / precision mesure ##
#===============================
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))
ggplot(obs_error, aes(x = reorder(site, hm), y = hmad_st)) +
geom_bar(stat = "identity")
#==============================
## H series ##
## heterogeneité ##
#===============================
pl <- ggplot(dat_h_veg, aes(x = date_releve, y = hmean)) +
### sd vs mad
ggplot(dat_h_veg, aes(x = hmad, y = hsd)) +
geom_point()
### 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)
pl + theme(legend.position = "none")
## boxplot par milieu ##
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")
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")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment