Nothing
#' Combine multiple phenology contour plots in one figure
#'
#' For multiple datasets, this function plots surface plots relating mean
#' temperatures during specified periods to annually recurring variables (e.g.
#' flowering). It produces one panel per dataset and plots them all in one
#' figure. Plots can be produced in R or directly exported as .png files.
#'
#' This function is only useful, if you want to plot several surface plots in
#' the same figure. These must relate to the same weather dataset. Arguably,
#' this function isn't quite ready to be released, but it performs some useful
#' functions that you may be interested in...
#'
#' @param pheno_list a data.frame with the following columns: varieties
#' (contains a character string), Start_chill (the start of the chill period,
#' in Julian days), End_chill (the end of the chill period, in Julian days),
#' Start_heat (the start of the forcing period, in Julian days), End_heat (the
#' end of the forcing period, in Julian days), Link (the complete path to a csv
#' file that contains all the annual observations for the dataset, with columns
#' Year and pheno)
#' @param fixed_weather daily weather, as produced with the fix_weather
#' function
#' @param split_month the month after which to start a new season. Defaults to
#' 6, meaning the new season will start in July.
#' @param outpath Path to the folder where the images should be saved. Should
#' include a trailing "/". The folder must already exists.
#' @param file_name name of the image file to be produced, if image_type='png'.
#' @param image_type Character string indicating the file format that should be
#' output. Image files are only produced for the moment, if this is "png". All
#' other values, as well as the default NA lead to output as an R plot only.
#' @param fonttype The type of font to be used for the figures. Can be 'serif'
#' (default) for a Times New Roman like font, 'sans' for an Arial type font or
#' 'mono' for a typewriter type font.
#' @param percol number of plots to be placed in a column.
#' @param xlabel label for the x-axis (if unhappy with the default).
#' @param ylabel label for the y-axis (if unhappy with the default).
#' @param height_factor height of the resulting png figure (if this is a png)
#' relative to the width of the plot (e.g. 1 or 0.7, defaults to 0.8).
#' @return Only a side effect is produced: either a .png file or an R graphic
#' showing the multi-panel contour figure.
#' @author Eike Luedeling
#' @keywords utility
#' @examples
#'
#' #this example uses arbitrarily modified versions of the KA_bloom dataset, and the starts
#' #end ends of the periods are also arbitraty. So the outputs may not make a lot of sense...
#'
#' weather<-fix_weather(KA_weather[which(KA_weather$Year>2000),])
#' pheno_list<-data.frame(varieties=c("KA1","KA2","KA3","KA4"), Start_chill=c(270,305,315,320),
#' End_chill=c(15,20,35,40), Start_heat=c(17,25,40,45),End_heat=c(90,100,110,115),
#' Link=c("KA1.csv","KA2.csv","KA3.csv","KA4.csv"))
#'
#' # write.csv(KA_bloom,"KA1.csv",row.names=FALSE)
#' KA_bloom$pheno<-as.numeric(as.character(KA_bloom$pheno))+10
#' # write.csv(KA_bloom,"KA2.csv",row.names=FALSE)
#' KA_bloom$pheno<-KA_bloom$pheno+10
#' # write.csv(KA_bloom,"KA3.csv",row.names=FALSE)
#' KA_bloom$pheno<-KA_bloom$pheno+10
#' # write.csv(KA_bloom,"KA4.csv",row.names=FALSE)
#'
#' # make_multi_pheno_trend_plot(pheno_list,weather, split_month=6,
#' # outpath=NA,file_name=NA,image_type="",fonttype="serif",percol=2)
#'
#'
#' @export make_multi_pheno_trend_plot
make_multi_pheno_trend_plot <-
function(pheno_list,fixed_weather,
split_month=6, #last month in same year
outpath=NA,
file_name=NA,
image_type="png",
fonttype="serif",percol=5,xlabel=NA,ylabel=NA,
height_factor=0.8)
{
read_pheno_flex<-function(phfile)
{csv<-read.csv(phfile)
if(length(colnames(csv))>=2&length(grep("Year",colnames(csv),ignore.case=TRUE)>0)&
length(grep("Pheno",colnames(csv),ignore.case=TRUE))>0)
out<-csv else
{csv<-read.table(phfile,sep=";",header=TRUE)
if(length(colnames(csv))>=2&length(grep("Year",colnames(csv),ignore.case=TRUE)>0)&
length(grep("Pheno",colnames(csv),ignore.case=TRUE))>0)
out<-csv else warning(paste(phfile,"is not a recognized pheno file"))}
colnames(out)<-c("Year","pheno")
out$Year<-as.numeric(out$Year)
out$pheno<-as.numeric(as.character(out$pheno))
return(out)
}
divable<-function(x,div=2)
{x/div==floor(x/div)}
weather<-fixed_weather$weather
#if(!use_Tmean) weather[,"Tmean"]<-(weather$Tmax+weather$Tmin)/2
stack<-stack_hourly_temps(fixed_weather,34.71)$hourtemps
stack[which(stack$Month<=split_month),"Season"]<-stack[which(stack$Month<=split_month),"Year"]
stack[which(stack$Month>split_month),"Season"]<-stack[which(stack$Month>split_month),"Year"]+1
stack[,"Date"]<-stack$Month*100+stack$Day
stack[,"JDay"]<-strptime(paste(stack$Month,"/",stack$Day,"/",stack$Year,sep=""),"%m/%d/%Y")$yday+1
sea<-unique(stack$Season)
res<-data.frame(Season=sea,Chill_Tmean=NA,Heat_Tmean=NA,Year_Tmean=NA)
k<-list()
chillmax<-chillmin<-heatmax<-heatmin<-phenomin<-phenomax<-NA
for (cult in 1:nrow(pheno_list))
{
sea<-unique(stack$Season)
res<-data.frame(Season=sea,Chill_Tmean=NA,Heat_Tmean=NA,Year_Tmean=NA)
if(pheno_list$Start_chill[cult]>pheno_list$End_chill[cult])
chill_days<-c(pheno_list$Start_chill[cult]:366,1:pheno_list$End_chill[cult]) else
chill_days<-pheno_list$Start_chill[cult]:pheno_list$End_chill[cult]
if(pheno_list$Start_heat[cult]>pheno_list$End_heat[cult])
heat_days<-c(pheno_list$Start_heat[cult]:366,1:pheno_list$End_heat[cult]) else
heat_days<-pheno_list$Start_heat[cult]:pheno_list$End_heat[cult]
for (s in sea)
{res[which(res$Season==s),"Chill_Tmean"]<-mean(stack[which(stack$Season==s&stack$JDay %in% chill_days),"Temp"])
res[which(res$Season==s),"Heat_Tmean"]<-mean(stack[which(stack$Season==s&stack$JDay %in% heat_days),"Temp"])
res[which(res$Season==s),"Year_Tmean"]<-mean(stack[which(stack$Season==s),"Temp"])}
pheno<-suppressWarnings(read_pheno_flex(as.character(pheno_list$Link[cult])))
for (i in 1:nrow(pheno))
{if(pheno[i,1] %in% res$Season)
{pheno[i,"Chill_Tmean"]<-res[which(res$Season==pheno[i,1]),"Chill_Tmean"]
pheno[i,"Heat_Tmean"]<-res[which(res$Season==pheno[i,1]),"Heat_Tmean"]
pheno[i,"Year_Tmean"]<-res[which(res$Season==pheno[i,1]),"Year_Tmean"]
}}
pheno[,2]<-as.numeric(as.character(pheno[,2]))
pheno<-pheno[which(!is.na(pheno[,"pheno"])),]
pheno<-pheno[which(!is.na(pheno[,"Chill_Tmean"])),]
pheno<-pheno[which(!is.na(pheno[,"Heat_Tmean"])),]
pheno<-pheno[which(!is.na(pheno[,"Year_Tmean"])),]
k[[cult]]<-Krig(x=as.matrix(pheno[,c("Chill_Tmean","Heat_Tmean")]),Y=pheno$pheno)
chillmin<-min(c(chillmin,pheno$Chill_Tmean),na.rm=TRUE)
chillmax<-max(c(chillmax,pheno$Chill_Tmean),na.rm=TRUE)
heatmin<-min(c(heatmin,pheno$Heat_Tmean),na.rm=TRUE)
heatmax<-max(c(heatmax,pheno$Heat_Tmean),na.rm=TRUE)
phenomax<-max(c(phenomax,pheno$pheno),na.rm=TRUE)
phenomin<-min(c(phenomin,pheno$pheno),na.rm=TRUE)
}
if(is.na(xlabel)) xlabel<-paste("Mean temperature during the chilling phase (deg. C)",sep="")
if(is.na(ylabel)) ylabel<-paste("Mean temperature during the forcing phase (deg. C)",sep="")
if(image_type=="png")
{png(paste(outpath,file_name,".png",sep=""),width=1200,height=1200*height_factor)
lwds<-2;cexs<-3} else {lwds<-1;cexs<-1}
par(family=fonttype)
perrow<-ceiling(length(k)/percol)
par(mfcol=c(percol,perrow))
if(image_type=="png") par(oma=c(8,8,8,22)) else par(oma=c(4,5,4,14))
par(mar=c(0,0,0,0))
for(i in 1:length(k))
{par(mar=c(0,0,0,0))
image(predictSurface(k[[i]]),xlim=c(chillmin,chillmax),ylim=c(heatmin-(heatmax-heatmin)/6,heatmax),zlim=c(phenomin,phenomax),
xaxs="i",axes=FALSE,col = tim.colors(),asp=1)
if(image_type=="png") contour(predictSurface(k[[i]]),add=TRUE,nlevels=5,labcex=1.5) else
contour(predictSurface(k[[i]]),add=TRUE,nlevels=5)
box("plot",lwd=lwds)
if(image_type=="png") mtext(pheno_list$varieties[i],line=-2,side=1,cex=1.6) else
mtext(pheno_list$varieties[i],line=-1.2,side=1,cex=0.8)
if(!divable(length(k),percol)) #if last column isn't complete
{if(divable(ceiling(length(k)/percol))) # if number of columns is even
{ #label x-axis such that the last column has the label on top
if (divable(i-1,percol)&!divable(floor(i/percol))) #
axis(3,tck=-0.03,lwd=lwds,cex.axis=cexs)
if (divable(i,percol)&!divable(floor(i/percol)))
axis(1,tck=-0.03,lwd=lwds,cex.axis=cexs,padj=1)
} else
{ #label x-axis such that the last column has the label on top
if (divable(i,percol)&divable(floor(i/percol))) #
axis(1,tck=-0.03,lwd=lwds,cex.axis=cexs,padj=1)
if (divable(i-1,percol)&divable(floor(i/percol)))
axis(3,tck=-0.03,lwd=lwds,cex.axis=cexs) }
perrow<-ceiling(length(k)/percol)
if(divable(percol*perrow-length(k)))
{if(divable(percol)) if ((i<=percol)&!divable(i)) axis(2,tck=-0.03,lwd=lwds,cex.axis=cexs)
if(!divable(percol)) if ((i<=percol)&divable(i)) axis(2,tck=-0.03,lwd=lwds,cex.axis=cexs)
if(divable(percol)) if ((i>length(k)-percol)&divable(i)) axis(4,tck=-0.03,lwd=lwds,cex.axis=cexs,padj=1)
if(!divable(percol)) if ((i>length(k)-percol)&!divable(i-floor((i-1)/percol)*percol)) axis(4,tck=-0.03,lwd=lwds,cex.axis=cexs,padj=1)
} else
{if(divable(percol)) if ((i<=percol)&divable(i)) axis(2,tck=-0.03,lwd=lwds,cex.axis=cexs)
if(!divable(percol)) if ((i<=percol)&!divable(i)) axis(2,tck=-0.03,lwd=lwds,cex.axis=cexs)
if(divable(percol)) if ((i>length(k)-percol)&!divable(i)) axis(4,tck=-0.03,lwd=lwds,cex.axis=cexs,padj=1)
if(!divable(percol)) if ((i>length(k)-percol)&divable(i-floor((i-1)/percol)*percol)) axis(4,tck=-0.03,lwd=lwds,cex.axis=cexs,padj=1)}
} else
{ #label x-axis such that the first column has the label at the bottom
if (divable(i,percol)&!divable(i/percol)) axis(1,tck=-0.03,lwd=lwds,cex.axis=cexs,padj=1)
if (divable(i-1,percol)&!divable((i-1)/percol)) axis(3,tck=-0.03,lwd=lwds,cex.axis=cexs)
if ((i<=percol)&!divable(i)) axis(2,tck=-0.03,lwd=lwds,cex.axis=cexs)
if ((i>length(k)-percol)&divable(i-floor((i-1)/percol)*percol)) axis(4,tck=-0.03,lwd=lwds,cex.axis=cexs,padj=1)
}
axis(1,labels=FALSE,tck=0.02,lwd=lwds)
axis(2,labels=FALSE,tck=0.03,lwd=lwds)
axis(3,labels=FALSE,tck=0.03,lwd=lwds)
axis(4,labels=FALSE,tck=0.03,lwd=lwds)
}
if(image_type=="png")
{mtext(xlabel,side=1,outer=TRUE,line=6,cex=2.5)
mtext(xlabel,side=3,outer=TRUE,line=4,cex=2.5)
mtext(ylabel,side=2,outer=TRUE,line=5,cex=2.5)
mtext(ylabel,side=4,outer=TRUE,line=6,cex=2.5)} else
{mtext(xlabel,side=1,outer=TRUE,line=2.5,cex=1.3*cexs)
mtext(xlabel,side=3,outer=TRUE,line=2.5,cex=1.3*cexs)
mtext(ylabel,side=2,outer=TRUE,line=3,cex=1.3*cexs)
mtext(ylabel,side=4,outer=TRUE,line=3,cex=1.3*cexs)}
opar <- par(fig=c(0, 1, 0, 1), oma=c(0, 0, 0, 0), mar=c(0, 0, 0, 0), new=TRUE)
plot(0, 0, type='n', bty='n', xaxt='n', yaxt='n')
if(image_type=="png") {par(oma=c(12,0,12,4))
image.plot(predictSurface(k[[i]]),xlim=c(chillmin,chillmax),ylim=c(heatmin,heatmax),zlim=c(phenomin,phenomax),
xaxs="i",axes=FALSE,col = tim.colors(),legend.only=TRUE,axis.args=list(cex.axis=cexs),
legend.lab="Bloom date (day of the year)",legend.cex=2.5,legend.line=6,lwd=lwds)} else
{par(oma=c(4,0,4,0))
image.plot(predictSurface(k[[i]]),xlim=c(chillmin,chillmax),ylim=c(heatmin,heatmax),zlim=c(phenomin,phenomax),
xaxs="i",axes=FALSE,col = tim.colors(),legend.only=TRUE,axis.args=list(cex.axis=cexs),
legend.lab="Bloom date (day of the year)",legend.cex=1.3,legend.line=3,lwd=lwds)}
par(opar)
if(image_type=="png") dev.off()
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.