Skip to content
GitLab
Select Git revision
  • master default
  • BRIEFCOM
  • ctfs
  • test.dplyr
  • package
Find file BlameHistoryPermalink
CLIMATE.FRANCE.R 10.37 KiB
1 
################################################
2 
################################################
3 
################################################
4 
################################################
5 
### LOAD and PROCEED CLIMATIC DATA
6 
################################################
7 
################################################
8 
################################################
9 
################################################
10 
source("./R/FUN.climate.R")
11 
## Data extracted from C. Piedallu data base with
12 
## Radition with cloud cover for each plots
13 
## Piedallu, C, and J Gegout. 2008. “Efficient Assessment of Topographic Solar Radiation to Improve Plant Distribution Models.” Agricultural and Forest Meteorology 148 (11) (October): 1696–1706. doi:10.1016/j.agrformet.2008.06.001.
14 
## Temperature and precipitation with temperature corrected by elevation based on
15 
## Piedallu, C., J. C. Gégout, V. Perez, F. Lebourgeois, and R. Field. 2012. “Soil Water Balance Performs Better Than Climatic Water Variables in Tree Species Distribution Modelling.” Global Ecology and Biogeography. http://onlinelibrary.wiley.com/doi/10.1111/geb.12012/full.
16
17 
######
18 
## LOAD PIDEALLU DATA
19 
data.CLIM <- read.csv("./data/raw/DataFrance/climate_piedallu/placettesGK_avec_2011.csv",
20 
                      header=T,sep='\t',stringsAsFactors=FALSE,dec=",",na.strings="")
21
22 
data.CLIM$rumkg_500 <- as.numeric(gsub(",",".",data.CLIM$rumkg_500))
23
24
25 
#########################################
26 
#########################################
27 
#### COMPUTE CLIMATIC VARIABLES !!!
28
29 
### NEED TO COMPUTE DDG5
30 
### NEED TO COMPUTE SOIL MAX WATER CONTENT
31 
### NEED TO COMPUTE PET
32 
### NEED TO COMPUTE WATER BUDGET
33
34 
## remove NA lines
35 
data.CLIM.na <- is.na(data.CLIM[["tmoy6190_1_cor"]])
36 
data.CLIM2 <-  data.CLIM[!data.CLIM.na,]
37
38 
### apply function sgdd
39 
sgdd.vec <- apply(data.CLIM2[,82:93],MARGIN=1,FUN=fun.sgdd)
40 
MeanT.vec <- apply(data.CLIM2[,82:93],MARGIN=1,FUN=mean)
41 
plot(MeanT.vec,sgdd.vec)
42 
dim(data.CLIM)
43
44 
#################################################
45 
### Max Soil Water Content
46 
### compute based on
47 
## Piedallu, C., J. C. Gégout, A. Bruand, and I. Seynave. 2011. “Mapping Soil Water Holding Capacity over Large Areas to Predict Potential Production of Forest Stands.” Geoderma 160 (3): 355–366.
48
49 
### read data texture as in Piedallu et al. 2011
50 
data.texture <- read.table("./data/raw/DataFrance/climate_piedallu/texture.txt",header=T)
51 
### load ecological data
52 
load(file="./data/process/ecologie_tot.Rdata")
53 
head(ecologie_tot)
54
55 
codeprof <- c(2.5,10,20,30,40,50,60,70,80,92.5) # code prof en cm
56 
codecaillou <- c(2.5,10,20,30,40,50,60,70,80,90,97.5) ## code percentage en %
57
58 
#compute swhc
59 
swhc <- fun.swhc(affroc=ecologie_tot$affroc,cailloux=ecologie_tot$cailloux,text2=ecologie_tot$text2,text1=ecologie_tot$text1,prof2=ecologie_tot$prof2,prof1=ecologie_tot$prof1,codeprof,codecaillou,data.texture)
60 
swhc[is.na(ecologie_tot$prof2)] <-  NA
61
62 
### add to ecological data based
63 
ecologie_tot$SWHC <- swhc
64
65
66 
##################################
67 
##################################
68 
#### COMPUTE PET WITH TURC FORMULA
69 
### unit to convert radneb61_ from J/cm2/month in MJ/m2/day /100/30 then to convert in KJ/m2/day as in Nick formula *1000
70
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
## RAD and Temperature of the 12 months radneb61_1to12 <- 66:77 tmoy6190_1_cor.1to12 <- 82:93 ### apply in parallel library(doParallel) registerDoParallel(cores=6) ## affect automaticaly half of the core detected to the foreach getDoParWorkers() ## here 8 core so 4 core if want to use more registerDoParallel(cores=6) PET.matrix <- foreach(i=1:length(data.CLIM$idp), .combine=rbind) %dopar% { fun.PET(i,rad=data.CLIM[,radneb61_1to12],temp=data.CLIM[,tmoy6190_1_cor.1to12]) } PET.matrix <- PET.matrix[!data.CLIM.na,] PET.matrix[PET.matrix<0] <- 0 ## affect zero if negative PET ## plot to check PET computed by me vs PET of Christian ## par(mfrow=c(3,4)) ## for (i in 1:12) {plot(PET.matrix[,i],data.CLIM2[,i+48]) ; lines(0:100,0:100)} colnames(PET.matrix) <- paste("PET.cor.",1:12,sep="") ############## ### MERGE CLIMATE and PET data.CLIM2 <- cbind(data.CLIM2,PET.matrix,sgdd=sgdd.vec) ### MERGE WITH ECOLOGICAL DATA ecologie.clim <- merge(ecologie_tot,data.CLIM2,by="idp",all.x=T) dim(ecologie_tot) dim(ecologie.clim) names(ecologie.clim) ecologie.clim2 <- ecologie.clim[!is.na(ecologie.clim$PET.cor.1) &!is.na(ecologie.clim$SWHC),] ########################################################################################### ###### COMPUTE WATER BUDGET ## ### test function ## fun.WaterBudget(i=16,prcp.m=(ecologie.clim2[,paste("prec6190_",1:12,sep="")]), ## PET.m=(ecologie.clim2[,paste("PET.cor.",1:12,sep="")]), ## Ta.m=(ecologie.clim2[,c("tmoy6190_1_cor",paste("tmoy6190_1_cor.",1:11,sep=""))]), ## SWHC.v=(ecologie.clim2[ ,"SWHC"]),n=2) ### apply function in parallel library(doParallel) registerDoParallel(cores=6) ## affect automaticaly half of the core detected to the foreach getDoParWorkers() ## here 8 core so 4 core if want to use more registerDoParallel(cores=6) WBWS.matrix <- foreach(i=1:length(ecologie.clim2$idp), .combine=rbind) %dopar% { fun.WaterBudget(i,prcp.m=(ecologie.clim2[,paste("prec6190_",1:12,sep="")]), PET.m=(ecologie.clim2[,paste("PET.cor.",1:12,sep="")]), Ta.m=(ecologie.clim2[,c("tmoy6190_1_cor",paste("tmoy6190_1_cor.",1:11,sep=""))]), SWHC.v=(ecologie.clim2[ ,"SWHC"]),n=2) } ### MERGE all and saved WBWS.matrix2 <- cbind("idp"=ecologie.clim2$idp,WBWS.matrix) ecologie.clim.data <- merge(ecologie.clim,WBWS.matrix2,by="idp",all.x=T) ##### # change tplant not good format load("./data/process/placette_tot.Rdata") ecologie.clim.data$tplant <- placette_tot$tplant ### Compute mean T and annual sum of precip
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
ecologie.clim.data$MAT <- apply(ecologie.clim.data[,114:125],MARGIN=1,FUN=mean) ecologie.clim.data$SAP <- apply(ecologie.clim.data[,69:80],MARGIN=1,FUN=sum) saveRDS(ecologie.clim.data,file="./data/process/ecologie.clim.data.rds") ######################################################### ######################################################### ###### FIGURES OF CLIMATIC DATA ############################# ecologie.clim.data <- readRDS("./data/process/ecologie.clim.data.rds") names(ecologie.clim.data) ###check climatic data pdf("./figs/sgdd.tmin.map.pdf") par(mfrow=c(2,2)) plot(ecologie.clim.data$xl93,ecologie.clim.data$yl93, col=rev(heat.colors(10))[cut(ecologie.clim.data$sgdd,quantile(ecologie.clim.data$sgdd,probs=(0:10)/10,na.rm=T),labels=F)], cex=0.2,pty="s",xlab=NA,ylab=NA,main="SGDD") plot(ecologie.clim.data$xl93,ecologie.clim.data$yl93, col=rev(heat.colors(10))[cut(ecologie.clim.data$tmin6190_min_cor,quantile(ecologie.clim.data$tmin6190_min_cor,probs=(0:10)/10,na.rm=T),labels=F)], cex=0.2,pty="s",xlab=NA,ylab=NA,main="Tmin") plot( ecologie.clim.data$tmin6190_min_cor,ecologie.clim.data$sgdd,cex=0.2,xlab="Tmin",ylab="SGDD") dev.off() pdf("./figs/Water.map.pdf") par(mfrow=c(2,2),mar=c(1,1,1,1)) plot(ecologie.clim.data$xl93,ecologie.clim.data$yl93, col=rev(colorRampPalette(c("blue", "red"))( 10 ))[cut(ecologie.clim.data$WB.y,quantile(ecologie.clim.data$WB.y,probs=(0:10)/10,na.rm=T),labels=F)], cex=0.2,pty="s",xlab=NA,ylab=NA) plot(ecologie.clim.data$xl93,ecologie.clim.data$yl93, col=rev(colorRampPalette(c("blue", "red"))( 10 ))[cut(ecologie.clim.data$WB.s,quantile(ecologie.clim.data$WB.s,probs=(0:10)/10,na.rm=T),labels=F)], cex=0.2,pty="s",xlab=NA,ylab=NA) plot(ecologie.clim.data$xl93,ecologie.clim.data$yl93, col=rev(colorRampPalette(c("blue", "red"))( 10 ))[cut(ecologie.clim.data$WS.y,quantile(ecologie.clim.data$WS.y,probs=(0:10)/10,na.rm=T),labels=F)], cex=0.2,pty="s",xlab=NA,ylab=NA) plot(ecologie.clim.data$xl93,ecologie.clim.data$yl93, col=rev(colorRampPalette(c("blue", "red"))( 10 ))[cut(ecologie.clim.data$WS.s,quantile(ecologie.clim.data$WS.s,probs=(0:10)/10,na.rm=T),labels=F)], cex=0.2,pty="s",xlab=NA,ylab=NA) dev.off() pdf("./figs/Water.var,cor.pdf") par(mfrow=c(2,2)) plot(ecologie.clim.data[['WB.y']],ecologie.clim.data[['WS.y']],xlab="WB year",ylab="ratio AET/D year",cex=0.2) plot(ecologie.clim.data[['WB.s']],ecologie.clim.data[['WS.s']],xlab="WB growing season",ylab="ratio AET/D growing season",cex=0.2) plot(ecologie.clim.data[['WB.s']],ecologie.clim.data[['WB.y']],xlab="WB growing season",ylab="WB growing year",cex=0.2) plot(ecologie.clim.data[['WS.s']],ecologie.clim.data[['WS.y']],xlab="ratio AET/D growing season",ylab="ratio AET/D year",cex=0.2) dev.off() pdf("./figs/sgdd.water.cor.pdf") plot(ecologie.clim.data[['sgdd']],ecologie.clim.data[['WS.y']],xlab="SGDD",ylab="ratio AET/D year",cex=0.2) dev.off() ### DO PCA OF CLIMATIC DATA data.prc <- prcomp(ecologie.clim.data[apply(is.na(ecologie.clim.data[,c('sgdd','tmin6190_min_cor','WB.s','WS.s')]),MARGIN=1,FUN=sum)<1,c('sgdd','tmin6190_min_cor','WB.s','WS.s')],na.action=na.omit,scale=T) summary(data.prc) scores <-data.prc$x pdf("./figs/climate.pca.pdf") biplot(data.prc,pch=1,xlabs=rep('.',length(data.prc$x[,1])),ylabs=c('SGDD','Tmin','WB','WS')) text(150,180,labels=paste("Var PC1 ", ((summary(data.prc)))$importance[2,1],sep="")) text(150,165,labels=paste("Var PC2 ", ((summary(data.prc)))$importance[2,2],sep="")) dev.off()
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243
##### FIGS on elevation correction of temperature done by Christian ### map of diff between mean altitude of the 1x1km cell and actual elevation of the plots pdf("./figs/map.error.alti.pdf") plot(data.CLIM$xl2,data.CLIM$yl2,col="grey",cex=0.3,xlab="X",ylab="Y") data.temp <- data.CLIM[data.CLIM$diff..altIFN...mnt50 < (-100) | (data.CLIM$diff..altIFN...mnt50>100),] points(data.temp$xl2,data.temp$yl2,col=c("red","coral3","coral","grey","coral","coral3","red")[cut(data.temp$diff..altIFN...mnt50,breaks=c(-500,-300,-200,-100,100,200,300,500),labels=F)],cex=0.5,pch=1) dev.off() ##### PLOT HIST OF ALL VARIABLES pdf('./figs/climatvar.hist.pdf') par(mfrow=c(3,4)) for (i in 37:48) hist(data.CLIM[,i],xlab='Monthly Precip (mm)',main=names(data.CLIM)[i]) par(mfrow=c(3,4)) for (i in 49:60) hist(data.CLIM[,i],xlab='PET (mm)',main=names(data.CLIM)[i]) par(mfrow=c(3,4)) for (i in 66:77) hist(data.CLIM[,i],xlab='Radiation (?)',main=names(data.CLIM)[i]) par(mfrow=c(3,4)) for (i in 82:93) hist(data.CLIM[,i],xlab='temperature (C)',main=names(data.CLIM)[i]) dev.off() pdf('./figs/xyplot.Tunco.Tcor.pdf') par(mfrow=c(3,4)) for (i in 1:12) plot(data.CLIM[,(23:34)[i]],data.CLIM[,(82:94)[i]],xlab='temperature (C) uncorrected',ylab='temperature (C) corrected', main=names(data.CLIM)[(23:34)[i]]) dev.off()