Commit 3a946e7e authored by Boulangeat Isabelle's avatar Boulangeat Isabelle
Browse files

prepa reu 15-10

parent e2f8c33c
extract_ncdat <- function(path, variables, NoPs){
# path = paste0(path_data_allslopes, "/reanalysis/pro/PRO_",as.character(year),"080106_",as.character(year+1),"080106.nc")
# path = "/Volumes/infogeo/Meteo_France/SAFRAN_montagne-Crocus_2020/alp_allslopes/reanalysis/meteo/FORCING_2017080106_2018080106.nc"
nc.read = nc_open(path)
# nc.read
# attributes(nc.read)$name
# attributes(nc.read$var)$name
# attributes(nc.read$dim)$name
# ncvar_get(nc.read, "Tair")->aa
time<-ncvar_get(nc.read,"time")
tunits<-ncatt_get(nc.read,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
dates<-as.Date(time,origin=unlist(tustr)[3])
nc.data = lapply(variables, function(x) ncvar_get(nc.read,x)[unique(NoPs),])
nc.data = abind(lapply(variables, function(x) ncvar_get(nc.read,x)[unique(NoPs),]), along = 3)
nc_close(nc.read)
return(list(data = abind(nc.data, along = 3), dates = dates))
return(list(data = nc.data, dates = dates))
}
calc_meteo_variables <- function(path_data_allslopes, year, days, NoPs, periods, tbase = 0){
# year = as.numeric(test$year[1])
# days = test$date_releve
# NoPs = test$NoP
variables = c("RN_ISBA", "TS_ISBA", "WSN_T_ISBA", "RAINF_ISBA", "TALB_ISBA")
# periods = c(300,600,900)
# library required
require(tidync)
require(ncmeta)
......@@ -25,6 +30,13 @@ calc_meteo_variables <- function(path_data_allslopes, year, days, NoPs, periods,
require(ncdf4.helpers)
require(abind)
# year = as.numeric(test$year[1])
# days = test$date_releve
# NoPs = test$NoP
variables = c("RN_ISBA", "TS_ISBA", "DSN_T_ISBA", "RAINF_ISBA", "TALB_ISBA")
variables_meteo = c("Tair", "HUMREL")
# periods = c(300,600,900)
# check filepath dataset
if(!(length(list.files(path_data_allslopes)) > 0) ) stop("not able to connect to data path or empty path")
......@@ -34,24 +46,32 @@ calc_meteo_variables <- function(path_data_allslopes, year, days, NoPs, periods,
datAll = abind(dat2$data, dat1$data, along = 2)
datesAll = c(dat2$dates, dat1$dates)
dat1_meteo = extract_ncdat(paste0(path_data_allslopes, "/reanalysis/meteo/FORCING_",as.character(year),"080106_",as.character(year+1),"080106.nc"), variables_meteo, NoPs)
dat2_meteo = extract_ncdat(paste0(path_data_allslopes, "/reanalysis/meteo/FORCING_",as.character(year-1),"080106_",as.character(year),"080106.nc"), variables_meteo, NoPs)
datAll_meteo = abind(dat2_meteo$data, dat1_meteo$data, along = 2)
datesAll_meteo = c(dat2$dates, dat1$dates)
# calc LSD = last snow day (10j no snow after longest snow period)
# other variables for each nop (sitexyear)
LSD = rep(NA, length(unique(NoPs)))
output = array(data = NA, dim = list(length(unique(NoPs)),8, length(periods)), dimnames = list(paste0(unique(NoPs), "_", year), c("gddspeed", "snowdays", "frost", "frost_severe", "radiations", "albedo", "rainfall", "hydro"), periods))
output = array(data = NA, dim = list(length(unique(NoPs)),8, length(periods)), dimnames = list(paste0(unique(NoPs), "_", year), c("gddspeed", "snowdays", "frost", "frost_severe", "radiations", "albedo", "rainfall", "temperature"), periods))
for(nop in 1:length(unique(NoPs))){
snow = datAll[nop,1:460 , which(variables == "WSN_T_ISBA")]
snow = datAll[nop,1:460 , which(variables == "DSN_T_ISBA")]
slide10_nosnow = which(mapply(function(X,Y) sum(snow[X:Y]) ,1:(length(snow)-9) , 10:length(snow))==0)
LSD[nop] = slide10_nosnow[which.max(diff(slide10_nosnow, lag=1))+1]
temp = datAll[nop,1:460, which(variables == "TS_ISBA")]-273.15
temp = datAll_meteo[nop,1:460, which(variables_meteo == "Tair")]-273.15
gdd = (temp[LSD[nop]:460] - tbase)
gdd[gdd<0]=0
gdd[snow[LSD[nop]:460]>0] = 0
cumgdd_series = cumsum(gdd)
output[nop,"temperature",] = mean(temp[LSD[nop]:460])
# number of days to reach treshold cumgdd
ndays = sapply(periods, function(x) which(cumgdd_series>x)[1])
......@@ -110,8 +130,8 @@ calc_meteo_variables <- function(path_data_allslopes, year, days, NoPs, periods,
nop = which(NoPs[rel]%in%unique(NoPs))
time_after_nosnow[rel] = d - LSD[nop]
cumGDD0[rel] = cumgdd(temp = datAll[nop,LSD[nop]:d, which(variables == "TS_ISBA")]-273.15, snow = datAll[nop,LSD[nop]:d, which(variables == "WSN_T_ISBA")], tbase=tbase)
cumGDD60d[rel] = cumgdd(temp = datAll[nop,(d-60):d, which(variables == "TS_ISBA")]-273.15, snow = datAll[nop,(d-60):d, which(variables == "WSN_T_ISBA")], tbase=tbase)
cumGDD0[rel] = cumgdd(temp = datAll_meteo[nop,LSD[nop]:d, which(variables_meteo == "Tair")]-273.15, snow = datAll[nop,LSD[nop]:d, which(variables == "DSN_T_ISBA")], tbase=tbase)
cumGDD60d[rel] = cumgdd(temp = datAll_meteo[nop,(d-60):d, which(variables_meteo == "Tair")]-273.15, snow = datAll[nop,(d-60):d, which(variables == "DSN_T_ISBA")], tbase=tbase)
}
return(list(var_gdd_treshold = output, last_snow = LSD, snowfree_time = time_after_nosnow, cumGDD0 = cumGDD0, cumGDD60d = cumGDD60d))
......
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
#==============================
## Import data from database ##
#===============================
require(RPostgreSQL)
config.db.siddt()
data=import_from_db()
saveRDS(data, "data.rds")
data = readRDS("data.rds")
\ No newline at end of file
sapply(list.files("R_fct"), function(x)source(paste0("R_fct/",x)))
#==============================
## Import data from database ##
#===============================
require(RPostgreSQL)
config.db.siddt()
data=import_from_db()
saveRDS(data, "data.rds")
data = readRDS("data.rds")
#==============================
## dataframe ##
#===============================
......@@ -43,27 +34,27 @@ dat_long <- dat_h_veg %>%
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")
## 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")
#
# 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")
#
#
#==============================
## effet obs / precision mesure ##
#===============================
......@@ -74,36 +65,36 @@ 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()
#
# 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")
#
# ggplot(obs_error, aes(x = reorder(site, hm), y = hmad_st)) +
# geom_bar(stat = "identity")
#==============================
## heterogeneité ##
#===============================
### 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)
## 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")
#
# ### 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)
#
# ## 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")
......@@ -74,14 +74,15 @@ releves$year= unlist(lapply(releves$date_releve, substr, start = 1,stop=4))
# summary(releves$day_number)
head(releves)
path_data_allslopes = "/Volumes/infogeo/Meteo_France/SAFRAN_montagne-Crocus_2020/alp_allslopes"
path_data_allslopes = "/Volumes/ISA-RESEARCH/alp_allslopes"
#path_data_allslopes = "/Volumes/ISA-RESEARCH/alp_allslopes"
# nc = nc_open(paste0(path_data_allslopes, "/reanalysis/pro/PRO_",as.character(year),"080106_",as.character(year+1),"080106.nc"))
# ncatt_get( nc, "TALB_ISBA")
library(dplyr)
# test = releves %>% select(releve, Number_of_points, year, date_releve) %>% filter(year == "2017")
# test = releves %>% select(releve, NoP, year, date_releve) %>% filter(year == "2017")
# test = test[1:6,]
#
# res = calc_meteo_variables(path_data_allslopes, as.numeric(test$year[1]), test$date_releve, test$NoP, c(30,60))
......@@ -181,14 +182,14 @@ head(reltot)
saveRDS(reltot, file = "releves_all_plus.rds")
saveRDS(sites_all, file = "sites_all_plus.rds")
###### analyse sites - type - topo
library(ggplot2)
ggplot(sites_all@data, aes(x = reorder(ref_typoveg, alti, na.rm=TRUE), y = alti)) +
geom_boxplot() +
labs(y="alti", x="Type vegetation")
ggplot(sites_all@data, aes(x = reorder(ref_typoveg, slope, na.rm=TRUE), y = slope)) +
geom_boxplot() +
labs(y="slope", x="Type vegetation")
# ###### analyse sites - type - topo
# library(ggplot2)
#
# ggplot(sites_all@data, aes(x = reorder(ref_typoveg, alti, na.rm=TRUE), y = alti)) +
# geom_boxplot() +
# labs(y="alti", x="Type vegetation")
#
#
# ggplot(sites_all@data, aes(x = reorder(ref_typoveg, slope, na.rm=TRUE), y = slope)) +
# geom_boxplot() +
# labs(y="slope", x="Type vegetation")
......@@ -6,6 +6,27 @@ length(unique(releves_all$id_site))
head(releves_all[which(is.na(releves_all$hmean)),])
library(ggplot2)
ggplot(releves_all, aes(x = reorder(ref_typoveg, cumgdd0, na.rm=TRUE), y = cumgdd0)) +
geom_boxplot() +
labs(y="gdd", x="Type vegetation")
ggplot(releves_all, aes(x = reorder(ref_typoveg, cumgdd0, na.rm=TRUE), y = last_snow_day.gdd300)) +
geom_boxplot() +
labs(y="LSD", x="Type vegetation")
ggplot(releves_all, aes(x = reorder(ref_typoveg, cumgdd0, na.rm=TRUE), y = gddspeed.gdd600)) +
geom_boxplot() +
labs(y="gddspeed 600", x="Type vegetation")
releves_all$date600 = releves_all$last_snow_day.gdd300 + releves_all$gddspeed.gdd600
str(releves_all)
ggplot(releves_all, aes(x = reorder(ref_typoveg, cumgdd0, na.rm=TRUE), y = date600)) +
geom_boxplot() +
labs(y="dateGDD600", x="Type vegetation")
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",
......
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