Commit 8ab10cda authored by Boulangeat Isabelle's avatar Boulangeat Isabelle
Browse files

analyse sites topo/ calc var clim

parent 0c650664
No preview for this file type
......@@ -37,7 +37,7 @@ library(tidyr)
library(ggplot2)
```
## Biomasse
```{r ,fig=TRUE}
library(dplyr)
......@@ -131,4 +131,31 @@ ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmad/hmean, na.rm=TRUE), y = hmad
geom_boxplot() +
labs(y="Variation intra de hauteur (stand. moy.)", x="Type vegetation")
```
\ No newline at end of file
```
## 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")
```
This diff is collapsed.
......@@ -12,7 +12,7 @@ always_allow_html: true
---
## Biomasse
```r
library(dplyr)
......@@ -189,3 +189,44 @@ ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmad/hmean, na.rm=TRUE), y = hmad
```
![](README_files/figure-html/unnamed-chunk-9-1.png)<!-- -->
## Sites topography
```
## Warning in readRDS("sites_all_plus.rds"): strings not representable in native
## encoding will be translated to UTF-8
```
```r
ggplot(sites_all@data, aes(x = reorder(ref_typoveg, alti, na.rm=TRUE), y = alti)) +
geom_boxplot() +
labs(y="alti", x="Type vegetation")
```
![](README_files/figure-html/unnamed-chunk-11-1.png)<!-- -->
```r
ggplot(sites_all@data, aes(x = reorder(ref_typoveg, slope, na.rm=TRUE), y = slope)) +
geom_boxplot() +
labs(y="slope", x="Type vegetation")
```
![](README_files/figure-html/unnamed-chunk-12-1.png)<!-- -->
```r
ggplot(sites_all@data, aes(x = reorder(ref_typoveg, northing, na.rm=TRUE), y = northing)) +
geom_boxplot() +
labs(y="northing", x="Type vegetation")
```
![](README_files/figure-html/unnamed-chunk-13-1.png)<!-- -->
```r
ggplot(sites_all@data, aes(x = reorder(ref_typoveg, easting, na.rm=TRUE), y = easting)) +
geom_boxplot() +
labs(y="easting", x="Type vegetation")
```
![](README_files/figure-html/unnamed-chunk-14-1.png)<!-- -->
......@@ -42,7 +42,11 @@ import_from_db <- function(){
paste(strsplit(x, "_")[[1]][-2], collapse="_")
}))
data_h = data_biomass_h %>% group_by(releve) %>% summarize(hmean=mean(hauteur), hsd=sd(hauteur), hmad = mad(hauteur))
### Attention! 3605 NA entries in hauteur
na.entries = data_biomass_h[which(is.na(data_biomass_h$hauteur)),]
unique(na.entries$ref_releve)
write.table(na.entries, file = "na_entries.txt", sep="\t", row.names=F, quote =F)
data_h = data_biomass_h %>% group_by(releve) %>% summarize(hmean=mean(hauteur), hsd=sd(hauteur), hmad = mad(hauteur), na.rm = TRUE)
# str(data_h)
# tail(data_h)
......
......@@ -10,12 +10,12 @@ extract_ncdat <- function(path, variables, NoPs){
}
calc_meteo_variables <- function(path_data_allslopes, year, days, NoPs, periods, tbase = 5.5){
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", "DSN_T_ISBA", "RAINF_ISBA", "TALB_ISBA")
# periods = c(30,60)
variables = c("RN_ISBA", "TS_ISBA", "WSN_T_ISBA", "RAINF_ISBA", "TALB_ISBA")
# periods = c(300,600,900)
# library required
require(tidync)
......@@ -34,92 +34,93 @@ 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)
output = array(data = NA, dim = list(length(days),8, length(periods)), dimnames = list(1:length(days), c("snowfree", "cumgdd", "frost", "gddspeed", "rainfall", "radiations", "albedo", "snowstock"), periods))
# releves loop
for (rel in 1:length(days)) {
# rel = 1
nop = which(NoPs[rel]%in%unique(NoPs))
# 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))
for(nop in 1:length(unique(NoPs))){
snow = datAll[nop,1:460 , which(variables == "WSN_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
gdd = (temp[LSD[nop]:460] - tbase)
gdd[gdd<0]=0
gdd[snow[LSD[nop]:460]>0] = 0
cumgdd_series = cumsum(gdd)
# number of days to reach treshold cumgdd
ndays = sapply(periods, function(x) which(cumgdd_series>x)[1])
output[nop,"gddspeed",] = ndays
# number of snow days after "winter" until threshold gdd
output[nop,"snowdays",] = sapply(ndays, function(x){
if(is.na(x)) out = NA
else{out = sum(snow[LSD[nop]:(LSD[nop]+x)]>0)}
out
} )
# number of frost days after "winter" until threshold gdd
output[nop,"frost",] = sapply(ndays, function(x){
if(is.na(x)) out = NA
else{out = sum(snow[LSD[nop]:(LSD[nop]+x)]<5 & temp[LSD[nop]:(LSD[nop]+x)]<0)}
out
} )
# number of severe frost days after "winter" until threshold gdd
output[nop,"frost_severe",] = sapply(ndays, function(x){
if(is.na(x)) out = NA
else{out = sum(snow[LSD[nop]:(LSD[nop]+x)]<5 & temp[LSD[nop]:(LSD[nop]+x)]<(-5))}
out
} )
# radiations sum until treshold gdd
output[nop,"radiations",] = sapply(ndays, function(x){
if(is.na(x)) out = NA
else{out = sum(datAll[nop,LSD[nop]:(LSD[nop]+x), which(variables == "RN_ISBA")])}
out
} )
# radiations sum until treshold gdd
output[nop,"albedo",] = sapply(ndays, function(x){
if(is.na(x)) out = NA
else{out = sum(datAll[nop,LSD[nop]:(LSD[nop]+x), which(variables == "TALB_ISBA")])}
out
} )
# precipitations sum until treshold gdd
output[nop,"rainfall",] = sapply(ndays, function(x){
if(is.na(x)) out = NA
else{out = sum(datAll[nop,LSD[nop]:(LSD[nop]+x), which(variables == "RAINF_ISBA")])}
out
} )
#
}
# sum GDD and time until measured date (for each rel)
# sum GDD 60 days before measure
time_after_nosnow = cumGDD0 = cumGDD60d = rep(NA, length(days))
for(rel in 1:length(days)){
d = which(days[rel]==datesAll)
nop = which(NoPs[rel]%in%unique(NoPs))
time_after_nosnow[rel] = d - LSD[nop]
for (p in 1:length(periods)){
# p = 1
temp = datAll[nop,(d-periods[p]):d, which(variables == "TS_ISBA")]-273.15
snow = datAll[nop,(d-periods[p]):d, which(variables == "DSN_T_ISBA")]
snowfree = snow==0
gdd = (temp - tbase)
gdd[gdd<0]=0
#### snowfree
output[rel, "snowfree", p] = sum(snowfree)
#### cumgdd (snowfree)
output[rel, "cumgdd", p] = sum(gdd[snowfree])
#### snowfreefrost
output[rel, "frost", p] = sum(which(temp<0 & snowfree))
#### gddspeed
output[rel, "gddspeed", p] = mean(gdd[gdd>0])
#### rainfall
output[rel, "rainfall", p] = sum(datAll[nop,(d-periods[p]):d, which(variables == "RAINF_ISBA")])
#### radiations
output[rel, "radiations", p] = sum(datAll[nop,(d-periods[p]):d, which(variables == "RN_ISBA")])
#### albedo
output[rel, "albedo", p] = sum(datAll[nop,(d-periods[p]):d, which(variables == "TABL_ISBA")])
#### snowstock
output[rel, "snowstock", p] = max(datAll[nop,(d-periods[p]):d, which(variables == "WSN_T_ISBA")])
}
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)
}
return(output)
return(list(var_gdd_treshold = output, last_snow = LSD, snowfree_time = time_after_nosnow, cumGDD0 = cumGDD0, cumGDD60d = cumGDD60d))
}
#=============================================
#--- build NoP raster functions
#=============================================
NoP_from_values <- function(x, safran_classes_complet, alti_field="ZS", slope_field="slope", aspect_field="aspect", massif_field="massif_num"){
num_point = safran_classes_complet[which(((safran_classes_complet[,alti_field] == x[alti_field] & safran_classes_complet[,slope_field] == x[slope_field]) & safran_classes_complet[,aspect_field]==x[aspect_field]) & safran_classes_complet[,massif_field]==x[massif_field]), "Number_of_points"]
if(length(num_point)==0) num_point=9999
return(num_point)
### function annexe
cumgdd <- function(temp, snow, tbase){
snowfree = snow==0
gdd = (temp - tbase)
gdd[gdd<0]=0
return(sum(gdd[snowfree]))
}
NoP_from_rasters <- function(r_alti, r_slope, r_aspect, r_massif, safran_classes_complet, parallel=FALSE, nclust=2){
if(resolution(r_alti)!=resolution(r_slope) || resolution(r_alti)!=resolution(r_aspect) || resolution(r_alti)!=resolution(r_massif)) stop("pb resolution")
if(projection(r_alti)!=projection(r_slope) || projection(r_alti)!=projection(r_aspect) || projection(r_alti)!=projection(r_aspect)) stop("pb proj")
crocus_alti <- cut(r_alti, breaks = seq(0, 5200, by = 300) , include.lowest = TRUE)
crocus_aspect <- cut(r_aspect, breaks = c(-1,seq(0, 360, by = 45)) , include.lowest = TRUE)
crocus_slope <- cut(r_slope, breaks = c(0, 10, 30, 90) , include.lowest = TRUE)
crocus_stack <- stack(crocus_alti, crocus_slope, crocus_aspect, crocus_massif_raster)
names(crocus_stack) <- c("altitude", "slope", "aspect", "massif")
corresp <- function(x, na.rm){
if(sum(is.na(x))==0){
num_point = safran_classes_complet[which(((sf4$cat_raster_alti == x[1] & safran_classes_complet$cat_raster_aspect == x[3]) & safran_classes_complet$cat_raster_slope==x[2]) & safran_classes_complet$massif_num==x[4]), "Number_of_points"]
}else num_point=NA
if(length(num_point)==0) num_point=9999
return(num_point)
}
if(parallel==TRUE){
library(doParallel)
beginCluster(n=nclust, type="SOCK")
cl = getCluster()
clusterExport(cl, "safran_classes_complet", envir=environment())
NoP = clusterR(crocus_stack, stackApply, args = list(indices=rep(1, 4), fun = corresp, na.rm=TRUE))
endCluster()
}else{
NoP <- stackApply(crocus_stack, indices=rep(1, 4),fun = corresp, na.rm=TRUE)
}
return(NoP)
}
This source diff could not be displayed because it is too large. You can view the blob instead.
No preview for this file type
......@@ -24,6 +24,7 @@ length(unique(dat_h_veg$releve))
dim(dat_h_veg)
summary(dat_h_veg$hmean)
dat_h_veg[is.na(dat_h_veg$hmean),]
#==============================
......@@ -76,7 +77,7 @@ cal_mean = unique(calage %>% group_by(releve) %>% summarize(site = site, hmean =
ggplot(cal_mean, aes(x = reorder(site, hmean), y = hmean)) +
geom_boxplot() +
geom_point()
geom_point()
obs_error = cal_mean %>% group_by(site) %>% summarize(hmad = mad(hmean), hmad_st = mad(hmean)/mean(hmean), hm = mean(hmean))
......@@ -89,7 +90,7 @@ ggplot(obs_error, aes(x = reorder(site, hm), y = hmad_st)) +
### sd vs mad
ggplot(dat_h_veg, aes(x = hmad, y = hsd)) +
geom_point()
geom_point()
### dans le temps
ggplot(dat_h_veg, aes(x = date_releve, y = hmad)) +
......@@ -106,4 +107,3 @@ ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmad, na.rm=TRUE), y = hmad)) +
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")
......@@ -21,9 +21,9 @@ coordinates(sites2)= c("x","y")
proj4string(sites2) <- CRS("+init=epsg:2154") # guess lambert 93
##
NoP_raster = raster("_DATA/safran_raster_NoP_stack.img")
NoP_raster = raster("_DATA/safran_raster_NoP.img")
#
# plot(NoP_raster)
projection(NoP_raster)==projection(sites1)
sites1_p = spTransform(sites1, projection(NoP_raster))
sites1_p$NoP = extract(NoP_raster, sites1_p)
......@@ -39,15 +39,21 @@ dim(sites_all)
length(unique(sites_all$NoP))
summary(sites_all$NoP)
saveRDS(sites_all, file = "sites_all.rds")
# plot(sites_all)
# points(sites_all[which(sites_all$NoP>9999),], col = "blue")
# sites_all[which(sites_all$NoP>9999),]@data
sf = read.csv("_DATA/safran_classes_complet.csv", sep=";", dec=",")
sites_all2 = merge(sites_all, sf[,c("Number_of_points", "NoP_stack")], by.x="NoP", by.y="NoP_stack")
tail(sites_all2)
saveRDS(sites_all2, file = "sites_all.rds")
sites_all = sites_all2
# sf = read.csv("_DATA/safran_classes_complet.csv", sep=";", dec=",")
# sites_all2 = merge(sites_all, sf[,c("Number_of_points", "NoP_stack")], by.x="NoP", by.y="NoP_stack")
# tail(sites_all2)
# sum(is.na(sites_all2$Number_of_points))
# sites_all2@data[is.na(sites_all2$Number_of_points),c("id_site", "NoP")]
#
# sf[which(sf$cat_raster_alti>8),"massif_num"]
# saveRDS(sites_all2, file = "sites_all.rds")
# sites_all = sites_all2
#=============================================
#--- Extract meteo data
# voir metadonnées https://en.aeris-data.fr/catalogue-en/ find "safran" S2M
......@@ -66,7 +72,7 @@ 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/infogeo/Meteo_France/SAFRAN_montagne-Crocus_2020/alp_allslopes"
path_data_allslopes = "/Volumes/ISA-RESEARCH/alp_allslopes"
......@@ -74,31 +80,114 @@ path_data_allslopes = "/Volumes/ISA-RESEARCH/alp_allslopes"
# ncatt_get( nc, "TALB_ISBA")
library(dplyr)
# test = releves %>% select(releve, NoP, year, date_releve) %>% filter(year == "2017")
# test = releves %>% select(releve, Number_of_points, 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))
# dim(res)
# dimnames(res)[[1]] = test$releve
releves2 = releves[-which(releves$year>2018),]
dim(releves)
dim(releves2)
releves2 = releves[-which(releves$year>2017),]
list_outputs <- lapply(unique(releves2$year), function(y){
print(y)
# y=2012
select = releves2 %>% select(releve, Number_of_points, year, date_releve) %>% filter(year == y)
res = calc_meteo_variables(path_data_allslopes, as.numeric(select$year[1]), select$date_releve, select$Number_of_points, c(30,60))
dimnames(res)[[1]] = select$releve
# y=2017
selected = releves2 %>% dplyr::select(releve, NoP, year, date_releve) %>% filter(year == y)
res = calc_meteo_variables(path_data_allslopes, as.numeric(selected$year[1]), selected$date_releve, selected$NoP, c(300, 600,900), tbase=0)
res
})
final_array = abind(list_outputs, along = 1)
names_releves_list = lapply(unique(releves2$year), function(y){
selected = releves2 %>% dplyr::select(releve, NoP, year, date_releve) %>% filter(year == y)
selected$releve
})
NoP_year_list = lapply(unique(releves2$year), function(y){
selected = releves2 %>% dplyr::select(releve, NoP, year, date_releve) %>% filter(year == y)
paste0(selected$NoP, "_", selected$year)
})
# res = list(var_gdd_treshold = output, last_snow = LSD, snowfree_time = time_after_nosnow, cumGDD0 = cumGDD0, cumGDD60d = cumGDD60d)
# 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))
length(list_outputs)
length(list_outputs[[1]])
dim(list_outputs[[1]][[1]])
length(list_outputs[[1]][[2]])
list_opt = lapply(list_outputs, function(x)x[[1]])
length(list_opt)
dim(list_opt[[1]])
last_snow_day = unlist(lapply(list_outputs, function(x)x[[2]]))
length(last_snow_day)
df_rel = data.frame(snowfree_time = unlist(lapply(list_outputs, function(x)x[[3]])), cumgdd0 = unlist(lapply(list_outputs, function(x)x[[4]])), cumgdd60d = unlist(lapply(list_outputs, function(x)x[[5]])))
rownames(df_rel) = unlist(names_releves_list)
final_array = abind(list_opt, along = 1)
dim(final_array)
head(final_array)
# build output (array to df)
df <- as.data.frame.table(final_array)
names(df) = c("releve", "meteo_variable", "period", "value")
releves_all = merge(releves, df, all=FALSE)
df_list = lapply(1:3, function(x){
df <- as.data.frame.table(cbind(final_array[,,x], last_snow_day))
df$gdd = c("gdd300", "gdd600", "gdd900")[x]
colnames(df) = c("NoP_year", "meteo_variable", "value", "gdd_threshold")
df
})
df_flat = do.call(rbind.data.frame, df_list)
head(df_flat)
tail(df_flat)
df_gdd = df_flat %>% unite("meteo", meteo_variable, gdd_threshold, sep = ".") %>% spread(meteo, value)
head(df_gdd)
## combine data.frame(x,y)
head(df_rel)
df_rel$NoP_year = unlist(NoP_year_list)
df_rel$releve = rownames(df_rel)
df_final = merge(df_rel, df_gdd, by = "NoP_year")
head(df_final)
releves_all = merge(releves, df_final, all=FALSE)
head(releves_all[order(releves_all$releve),])
saveRDS("releves_all", file = "releves_all.rds")
########
mnt25m <- raster("/Volumes/ISA-RESEARCH/_DATA/zaa_mnt/mntAlpes_25m.tif")
aspect_alps <- raster("/Volumes/ISA-RESEARCH/_DATA/zaa_CHABLI_DB/_spatialManips/aspect_trigo_alps.tif")
slope_alps <- raster("/Volumes/ISA-RESEARCH/_DATA/zaa_CHABLI_DB/_spatialManips/slopeAlps.tif")
lf = raster("/Volumes/ISA-RESEARCH/_DATA/zaa_landform/landform_tpi_saga_juill2020.tif")
projection(mnt25m)==projection(lf)
projection(sites_all)
plot(mnt25m)
points(sites_all)
sites_all$alti = extract(mnt25m, sites_all)
sites_all$slope = extract(slope_alps, sites_all)
sites_all$lf = extract(lf, sites_all)
library(shadow)
sites_all$northing = cos(deg2rad(extract(aspect_alps, sites_all)))
sites_all$easting = sin(deg2rad(extract(aspect_alps, sites_all)))
sites_all$aspect = extract(aspect_alps, sites_all)
summary(sites_all)
head(sites_all)
reltot = merge(sites_all@data[,c("id_site", "alti", "slope", "lf", "northing", "easting")], releves_all, by.x = "id_site", by.y="ref_site")
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")
releves_all = readRDS("releves_all.rds")
head(releves_all)
length(unique(releves_all$releve))
length(unique(releves_all$ref_site))
library(ggplot2)
library(randomForest)
sub = releves_all[releves_all$ref_typoveg=="SUB" & !(is.na(releves_all$hmean)),]
que = releves_all[releves_all$ref_typoveg=="QUE" & !(is.na(releves_all$hmean)),]
prod = releves_all[releves_all$ref_typoveg=="PROD" & !(is.na(releves_all$hmean)),]
alp = releves_all[releves_all$ref_typoveg=="ALP" & !(is.na(releves_all$hmean)),]
mod = randomForest(hmean ~ cumgdd.90 + frost.90 + rainfall.90, data = alp)
mod
ggplot(hmean~)
plot(hmean~cumgdd.90, data = releves_all)
dat = releves_all %>% filter(meteo_variable == "cumgdd")
ggplot(dat, aes(x=value, y=hmean)) +
library(ggplot2)
dat = ggplot(dat, aes(x=value, y=hmean)) +
geom_point() +
facet_wrap(~period, ncol = 2, scales = "free")
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