An error occurred while loading the file. Please try again.
-
Youcef Aouad authored4383599d
################################################################################################################################################################################################### LOAD and PROCEED CLIMATIC DATA################################################################################################################################################################################################source("./R/FUN.climate.R")## Data extracted from C. Piedallu data base with## Radition with cloud cover for each plots## 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.## Temperature and precipitation with temperature corrected by elevation based on## 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.######## LOAD PIDEALLU DATAdata.CLIM <- read.csv("./data/raw/DataFrance/climate_piedallu/placettesGK_avec_2011.csv", header=T,sep='\t',stringsAsFactors=FALSE,dec=",",na.strings="")data.CLIM$rumkg_500 <- as.numeric(gsub(",",".",data.CLIM$rumkg_500))###################################################################################### COMPUTE CLIMATIC VARIABLES !!!### NEED TO COMPUTE DDG5### NEED TO COMPUTE SOIL MAX WATER CONTENT### NEED TO COMPUTE PET### NEED TO COMPUTE WATER BUDGET## remove NA linesdata.CLIM.na <- is.na(data.CLIM[["tmoy6190_1_cor"]])data.CLIM2 <- data.CLIM[!data.CLIM.na,]### apply function sgddsgdd.vec <- apply(data.CLIM2[,82:93],MARGIN=1,FUN=fun.sgdd)MeanT.vec <- apply(data.CLIM2[,82:93],MARGIN=1,FUN=mean)plot(MeanT.vec,sgdd.vec)dim(data.CLIM)#################################################### Max Soil Water Content### compute based on## 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.### read data texture as in Piedallu et al. 2011data.texture <- read.table("./data/raw/DataFrance/climate_piedallu/texture.txt",header=T)### load ecological dataload(file="./data/process/ecologie_tot.Rdata")head(ecologie_tot)codeprof <- c(2.5,10,20,30,40,50,60,70,80,92.5) # code prof en cmcodecaillou <- c(2.5,10,20,30,40,50,60,70,80,90,97.5) ## code percentage en %#compute swhcswhc <- 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)swhc[is.na(ecologie_tot$prof2)] <- NA### add to ecological data based ecologie_tot$SWHC <- swhc######################################################################## COMPUTE PET WITH TURC FORMULA### 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 *10007172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
## 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()