path<-"/clima2/users/denis/monitmensal/chuva_ref/" pathfig<-"/clima2/users/denis/monitmensal/figuras/figurasnovas/" for (i in 1:124){ infile<-paste(path,"chuva_ref_",as.character(i),".txt",sep="") chuva<-read.table(infile,header=T) if(chuva <0){ chuva<-is.na(chuva)} serie<-as.matrix(chuva[1:30,2:13],na.rm=T) meses<-as.matrix(c("janeiro","fevereiro","marco","abril","maio","junho","julho","agosto","setembro","outubro","novembro","dezembro")) for(k in 1:12){ ##Densidade de probabilidade cdf<-ecdf(serie[,k]) r<-range(serie[,k]) den<-density(serie[,k]) ##Parametros da distribuicao gama alpha<-mean(serie[,k])^2/var(serie[,k]) beta<-var(serie[,k])/mean(serie[,k]) rate<-1/beta x<-sort(serie[,k]) yy<-sort(rank(serie[,k])/length(serie[,k])) #a<-rgamma(100,shape=alpha,rate=1/beta) #x<-seq(0,1000,0.1) desvio<-sd(serie[,k]) media<-mean(serie[,k]) ################## #-----Curvas-----# ################## pgamma<-curve(pgamma(x,shape=alpha,scale=beta),from=0,to=r[2]+250,n=3001) pgamma2<-curve(pgamma(x,shape=alpha,scale=beta),from=0,to=r[2]+150,n=3001) pgamma3<-curve(pgamma(x,shape=alpha,scale=beta),from=0,to=r[2]+50,n=3001) gammaexc<-curve(1-pgamma(x,shape=alpha,scale=beta),from=0,to=r[2]+150) gammaexc2<-curve(1-pgamma(x,shape=alpha,scale=beta),from=0,to=r[2]+150) gammaexc3<-curve(1-pgamma(x,shape=alpha,scale=beta),from=0,to=r[2]+100) dgamma<-curve(dgamma(x,shape=alpha,scale=beta),from=0,to=r[2]+150,n=6001) exdf<-curve(1-cdf(x),from=0,to=r[2]+200) dnormal<-curve(dnorm(x,sd=desvio,mean=media),from=r[1]-150,to=r[2]+150,n=3001) minden<-min(dgamma$y) maxden<-range(dgamma$y) maxden2<-maxden[2]/1.2 if(maxden[1] <0){ maxden[1]<-0.2 } ########### ##Quantis## ########### quantil5<-quantile(serie[,k],probs=.05) quantil15<-quantile(serie[,k],probs=.15) quantil33<-quantile(serie[,k],probs=.33) quantil66<-quantile(serie[,k],probs=.66) quantil85<-quantile(serie[,k],probs=.85) quantil95<-quantile(serie[,k],probs=.95) quantil87<-quantile(serie[,k],probs=.875) quantil62<-quantile(serie[,k],probs=.625) quantil37<-quantile(serie[,k],probs=.375) quantil12<-quantile(serie[,k],probs=.125) quantil75<-quantile(serie[,k],probs=.75) quantil25<-quantile(serie[,k],probs=.25) tercilinferior<-0.3333333333333333 tercilsuperior<-0.6666666666666666 #assimetria<-momentskew(serie[,k]) curtose<-((quantil87-quantil62)+(quantil37-quantil12))/(quantil75-quantil25) #Parametros da serie de dados parametros<-list(N=length(serie[,k]),mediana=median(serie[,k]),media=mean(serie[,k]),dp=sd(serie[,k]),min=min(serie[,k]),max=max(serie[,k])) if(parametros$min <0){ parametros$min<-12.0} #------------------------------------------ # DISTRIBUICAO ACUMULADA DE PROBABILIDADE # #------------------------------------------ # png(paste(pathfig,"distrib_acumulada_regiao_",as.character(i),"_",as.character(meses[k,]),".png",sep=""),width = 640, height = 660) # par(las=1,mar=c(4,4,4,7.5)) # plot(pgamma,type="l",from=0,to=r2+200,col="black",cex.main=1.5,font.lab=2,cex.lab=1.4,cex.axis=1.4,lwd=2.2,ylim=c(0.0,1.0),main=paste("Distribuição acumulada \nno mês de março"," (1981-2010): Região ",as.character(i),sep=""),xlab="Precipitação acumulada (mm) ",ylab="Probabilidade acumulada",xaxs="i") # lines(cdf,type="l",col="red",xlim=r,ylim=c(0.0,1.0),lwd=2,do.p=F,verticals=T) # grid(col=80,lty=3) # } # mtext(side=1,cex=1,adj=1,padj=4,text=" (Produto:CPTEC/INPE)") # legend("topleft",cex=1.1,xpd=TRUE,legend=c("Curva ajustada (função Gama)",paste("Curva empírica",sep="")),lty=1,col=c(1,2),lwd=2,bty="o",ncol=1) # mtext(at=0.95,side=4,col="dimgray",font=2,text=" Parâmetros",cex=1) # mtext(at=0.90,side=4,col="dimgray",font=2,text=" da série (mm)",cex=1) # mtext(at=0.85,side=4,col="dimgray",text=c(paste(" Média: ",sprintf("%.1f",parametros$media),sep="")),cex=1) # mtext(at=0.80,side=4,col="dimgray",text=c(paste(" Mediana: ",sprintf("%.1f",parametros$mediana),sep="")),cex=1) # mtext(at=0.75,side=4,col="dimgray",text=c(paste(" Desvio \n padrão: ",sprintf("%.1f",parametros$dp),sep="")),cex=1) # mtext(at=0.70,side=4,col="dimgray",text=c(paste(" Mín: ",sprintf("%.1f",parametros$min),sep="")),cex=1) # mtext(at=0.65,side=4,col="dimgray",text=c(paste(" Máx: ",sprintf("%.1f",parametros$max),sep="")),cex=1) # mtext(at=0.50,side=4,font=2,text=" Percentis(mm)",col="dimgray",cex=1) # mtext(at=0.45,side=4,col="dimgray",text=c(paste(" 05% - ",sprintf("%.1f",quantil5),sep="")),cex=1) # mtext(at=0.40,side=4,col="dimgray",text=c(paste(" 15% - ",sprintf("%.1f",quantil15),sep="")),cex=1) # mtext(at=0.35,side=4,col="dimgray",text=c(paste(" 33% - ",sprintf("%.1f",quantil33),sep="")),cex=1) # mtext(at=0.30,side=4,col="dimgray",text=c(paste(" 66% - ",sprintf("%.1f",quantil66),sep="")),cex=1) # mtext(at=0.25,side=4,col="dimgray",text=c(paste(" 85% - ",sprintf("%.1f",quantil85),sep="")),cex=1) # mtext(at=0.20,side=4,col="dimgray",text=c(paste(" 95% - ",sprintf("%.1f",quantil95),sep="")),cex=1) # mtext(at=0.10,side=4,col="dimgray",text=c(paste(" Curtose: ",sprintf("%.1f",curtose),sep="")),cex=1) # mtext(at=0.05,side=4,col="dimgray",text=c(paste(" Assimetria: ",sprintf("%.1f",assimetria),sep="")),cex=1) #Texto Eixo Y # mtext(c(0,0.15),cex=1.2,font=2,side=2,line=0.5,at=0.05,text="05%",col="dimgray") # mtext(c(0,0.15),cex=1.2,font=2,side=2,line=0.5,at=0.15,text="15%",col="dimgray") # mtext(c(0,0.3333333333333333),cex=1.2,font=2,side=2,line=0.5,at=tercilinferior,text="33%",col="dimgray") # mtext(c(0,0.6666666666666666),cex=1.2,font=2,side=2,line=0.5,at=tercilsuperior,text="66%",col="dimgray") # mtext(c(0,0.85),cex=1.2,font=2,side=2,line=0.5,at=0.85,text="85%",col="dimgray") # mtext(c(0,0.95),cex=1.2,font=2,side=2,line=0.5,at=0.95,text="95%",col="dimgray") #Linhas verticais # vector<-pgamma(x,shape=alpha,rate=1/beta) # index <- (1:length(vector))[vector > 0.65 & vector < 0.66] # if(length(index)>1){segments(x0=mean(x[index]),y0=0,x1=mean(x[index]),y1=0.66,lty=2,lwd=1.5)} # if(length(index)==1){segments(x0=x[index[1]],y0=0,x1=x[index[1]],y1=0.66,lty=2,lwd=1.5)} # text(x=mean(x[index]),cex=1.2,y=-0.01,col="dimgray",font=2,label="66%") # segments(x0=0,y0=2/3,x1=mean(x[index]),y1=2/3,lty=2,lwd=1.5) # vector<-pgamma(x,shape=alpha,rate=1/beta) # index <- (1:length(vector))[vector > 0.31 & vector < 0.35] # if(length(index)>1){segments(x0=mean(x[index]),y0=0,x1=mean(x[index]),y1=1/3,lty=2,lwd=1.5)} # if(length(index)==1){segments(x0=x[index[1]],y0=0,x1=x[index[1]],y1=1/3,lty=2,lwd=1.5)} # text(x=mean(x[index]),cex=1.2,y=-0.01,col="dimgray",font=2,label="33%") # segments(x0=0,y0=1/3,x1=mean(x[index]),y1=1/3,lty=2,lwd=1.5) # vector<-pgamma(x,shape=alpha,rate=1/beta) # index <- (1:length(vector))[vector > 0.03 & vector < 0.07] # if(length(index)>1){segments(x0=mean(x[index]),y0=0,x1=mean(x[index]),y1=0.05,lty=2,lwd=1.5)} # if(length(index)==1){segments(x0=x[index[1]],y0=0,x1=x[index[1]],y1=0.05,lty=2,lwd=1.5)} # text(x=mean(x[index]),cex=1.2,y=-0.01,col="dimgray",font=2,label="05%") # segments(x0=0,y0=0.05,x1=mean(x[index]),y1=0.05,lty=2,lwd=1.5) # vector<-pgamma(x,shape=alpha,rate=1/beta) # index <- (1:length(vector))[vector > 0.14 & vector < 0.16] # if(length(index)>1){segments(x0=mean(x[index]),y0=0,x1=mean(x[index]),y1=0.15,lty=2,lwd=1.5)} # if(length(index)==1){segments(x0=x[index[1]],y0=0,x1=x[index[1]],y1=0.15,lty=2,lwd=1.5)} # text(x=mean(x[index]),cex=1.2,y=-0.01,col="dimgray",font=2,label="15%") # segments(x0=0,y0=0.15,x1=mean(x[index]),y1=0.15,lty=2,lwd=1.5) # vector<-pgamma(x,shape=alpha,rate=1/beta) # index <- (1:length(vector))[vector > 0.84 & vector < 0.86] # if(length(index)>1){segments(x0=mean(x[index]),y0=0,x1=mean(x[index]),y1=0.85,lty=2,lwd=1.5)} # if(length(index)==1){segments(x0=x[index[1]],y0=0,x1=x[index[1]],y1=0.85,lty=2,lwd=1.5)} # text(x=mean(x[index]),cex=1.2,y=-0.01,col="dimgray",font=2,label="85%") # segments(x0=0,y0=0.85,x1=mean(x[index]),y1=0.85,lty=2,lwd=1.5) # vector<-pgamma(x,shape=alpha,rate=1/beta) # index <- (1:length(vector))[vector > 0.94 & vector < 0.96] # if(length(index)>1){segments(x0=mean(x[index]),y0=0,x1=mean(x[index]),y1=0.95,lty=2,lwd=1.5)} # if(length(index)==1){segments(x0=x[index[1]],y0=0,x1=x[index[1]],y1=0.95,lty=2,lwd=1.5)} # text(x=mean(x[index]),cex=1.2,y=-0.01,col="dimgray",font=2,label="95%") # segments(x0=0,y0=0.95,x1=mean(x[index]),y1=0.95,lty=2,lwd=1.5) # dev.off() #----------------------------- # DISTRIBUICAO DE EXCEDENCIA # #----------------------------- # png(paste(pathfig,"distrib_excedencia_regiao_",as.character(i),"_",as.character(meses[k,]),".png",sep=""),width = 640, height=660) # par(las=1,mar=c(4,4,4,7.5)) # plot(gammaexc,cex.main=1.5,font.lab=2,cex.lab=1.4,cex.axis=1.4,col="black",lwd=2.2,type="l",ylim=c(0.0,1.0),main=paste("Distribuição de probabilidade de excedência \nno mês de março"," (1981-2010): Região ",as.character(i),sep=""),xaxs="i",xlab="Precipitação acumulada (mm) ",ylab="Probabilidade de excedência") # lines(exdf,lwd=2,type="l",col="red",xlim=r,ylim=c(0.0,1.0),do.p=F,verticals=T) # grid(col=80,lty=3) # mtext("Probabilidade de excedência",font=2,cex=1.4,side=2,las=3,line=4) # if(r[2]<=300){ # plot(gammaexc3,cex.main=1.5,font.lab=2,cex.lab=1.4,cex.axis=1.4,col="black",lwd=2.2,type="l",ylim=c(0.0,1.0),main=paste("Distribuição de probabilidade de excedência \nno mês de março"," (1981-2010): Região ",as.character(i),sep=""),xaxs="i",xlab="Precipitação acumulada (mm) ",ylab="Probabilidade de excedência") # lines(exdf,lwd=2,type="l",col="red",xlim=r,ylim=c(0.0,1.0),do.p=F,verticals=T) # grid(col=80,lty=3) # mtext("Probabilidade de excedência",font=2,cex=1.4,side=2,las=3,line=4) # } # legend("topright",cex=1.1,xpd=TRUE,legend=c("Curva ajustada (função Gama)",paste("Curva empírica",sep="")),lty=1,col=c(1,2),lwd=2,bty="o",ncol=1) # mtext(side=1,cex=1,adj=1,padj=4,text=" (Produto:CPTEC/INPE)") # mtext(at=0.95,side=4,col="dimgray",font=2,text=" Parâmetros",cex=1) # mtext(at=0.90,side=4,col="dimgray",font=2,text=" da série (mm)",cex=1) # mtext(at=0.85,side=4,col="dimgray",text=c(paste(" Média: ",sprintf("%.1f",parametros$media),sep="")),cex=1) # mtext(at=0.80,side=4,col="dimgray",text=c(paste(" Mediana: ",sprintf("%.1f",parametros$mediana),sep="")),cex=1) # mtext(at=0.75,side=4,col="dimgray",text=c(paste(" Desvio \n padrão: ",sprintf("%.1f",parametros$dp),sep="")),cex=1) # mtext(at=0.70,side=4,col="dimgray",text=c(paste(" Mín: ",sprintf("%.1f",parametros$min),sep="")),cex=1) # mtext(at=0.65,side=4,col="dimgray",text=c(paste(" Máx: ",sprintf("%.1f",parametros$max),sep="")),cex=1) # mtext(at=0.50,side=4,font=2,text=" Percentis(mm)",col="dimgray",cex=1) # mtext(at=0.45,side=4,col="dimgray",text=c(paste(" 05% - ",sprintf("%.1f",quantil5),sep="")),cex=1) # mtext(at=0.40,side=4,col="dimgray",text=c(paste(" 15% - ",sprintf("%.1f",quantil15),sep="")),cex=1) # mtext(at=0.35,side=4,col="dimgray",text=c(paste(" 33% - ",sprintf("%.1f",quantil33),sep="")),cex=1) # mtext(at=0.30,side=4,col="dimgray",text=c(paste(" 66% - ",sprintf("%.1f",quantil66),sep="")),cex=1) # mtext(at=0.25,side=4,col="dimgray",text=c(paste(" 85% - ",sprintf("%.1f",quantil85),sep="")),cex=1) # mtext(at=0.20,side=4,col="dimgray",text=c(paste(" 95% - ",sprintf("%.1f",quantil95),sep="")),cex=1) #Texto Eixo Y # mtext(c(0,0.15),cex=1.2,font=2,side=2,line=0.5,at=0.05,text="05%",col="dimgray") # mtext(c(0,0.15),cex=1.2,font=2,side=2,line=0.5,at=0.15,text="15%",col="dimgray") # mtext(c(0,0.3333333333333333),cex=1.2,font=2,side=2,line=0.5,at=tercilinferior,text="33%",col="dimgray") # mtext(c(0,0.6666666666666666),cex=1.2,font=2,side=2,line=0.5,at=tercilsuperior,text="66%",col="dimgray") # mtext(c(0,0.85),cex=1.2,font=2,side=2,line=0.5,at=0.85,text="85%",col="dimgray") # mtext(c(0,0.95),cex=1.2,font=2,side=2,line=0.5,at=0.95,text="95%",col="dimgray") #Linhas verticais # vector<-1-pgamma(x,shape=alpha,scale=beta) # index <- (1:length(vector))[vector > 0.64 & vector < 0.68] # if(length(index)>1){segments(x0=mean(x[index]),y0=0,x1=mean(x[index]),y1=0.66,lty=2,lwd=1)} # if(length(index)==1){segments(x0=x[index[1]],y0=0,x1=x[index[1]],y1=0.66,lty=2,lwd=1)} # text(x=mean(x[index]),cex=1.2,y=-0.01,col="dimgray",font=2,label="33%") # segments(x0=0,y0=2/3,x1=mean(x[index]),y1=2/3,lty=2,lwd=1) # vector<-1-pgamma(x,shape=alpha,scale=beta) # index <- (1:length(vector))[vector > 0.31 & vector < 0.35] # if(length(index)>1){segments(x0=mean(x[index]),y0=0,x1=mean(x[index]),y1=1/3,lty=2,lwd=1)} # if(length(index)==1){segments(x0=x[index[1]],y0=0,x1=x[index[1]],y1=1/3,lty=2,lwd=1)} # text(x=mean(x[index]),cex=1.2,y=-0.01,col="dimgray",font=2,label="66%") # segments(x0=0,y0=1/3,x1=mean(x[index]),y1=1/3,lty=2,lwd=1) # vector<-1-pgamma(x,shape=alpha,scale=beta) # index <- (1:length(vector))[vector > 0.03 & vector < 0.07] # if(length(index)>1){segments(x0=mean(x[index]),y0=0,x1=mean(x[index]),y1=0.05,lty=2,lwd=1)} # if(length(index)==1){segments(x0=x[index[1]],y0=0,x1=x[index[1]],y1=0.05,lty=2,lwd=1)} # text(x=mean(x[index]),cex=1.2,y=-0.01,col="dimgray",font=2,label="95%") # segments(x0=0,y0=0.05,x1=mean(x[index]),y1=0.05,lty=2,lwd=1) # vector<-1-pgamma(x,shape=alpha,scale=beta) # index <- (1:length(vector))[vector > 0.13 & vector < 0.17] # if(length(index)>1){segments(x0=mean(x[index]),y0=0,x1=mean(x[index]),y1=0.15,lty=2,lwd=1)} # if(length(index)==1){segments(x0=x[index[1]],y0=0,x1=x[index[1]],y1=0.15,lty=2,lwd=1)} # text(x=mean(x[index]),cex=1.2,y=-0.01,col="dimgray",font=2,label="85%") # segments(x0=0,y0=0.15,x1=mean(x[index]),y1=0.15,lty=2,lwd=1) # vector<-1-pgamma(x,shape=alpha,scale=beta) # index <- (1:length(vector))[vector > 0.83 & vector < 0.87] # if(length(index)>1){segments(x0=mean(x[index]),y0=0,x1=mean(x[index]),y1=0.85,lty=2,lwd=1)} # if(length(index)==1){segments(x0=x[index[1]],y0=0,x1=x[index[1]],y1=0.85,lty=2,lwd=1)} # text(x=mean(x[index]),cex=1.2,y=-0.01,col="dimgray",font=2,label="15%") # segments(x0=0,y0=0.85,x1=mean(x[index]),y1=0.85,lty=2,lwd=1) # vector<-1-pgamma(x,shape=alpha,scale=beta) # index <- (1:length(vector))[vector > 0.93 & vector < 0.97] # if(length(index)>1){segments(x0=mean(x[index]),y0=0,x1=mean(x[index]),y1=0.95,lty=2,lwd=1)} # if(length(index)==1){segments(x0=x[index[1]],y0=0,x1=x[index[1]],y1=0.95,lty=2,lwd=1)} # text(x=mean(x[index]),cex=1.2,y=-0.01,col="dimgray",font=2,label="05%") # segments(x0=0,y0=0.95,x1=mean(x[index]),y1=0.95,lty=2,lwd=1) # dev.off() #------------------------------------------------ # HISTOGRAMA E CURVA DENSIDADE DE PROBABILIDADE # #------------------------------------------------ png(paste(pathfig,"histograma_regiao_",as.character(i),"_",as.character(meses[k,]),".png",sep=""),width = 640, height = 660) par(las=1,mar=c(4,6,4,7.5)) hist(serie[,k],cex.lab=1.3,cex.main=1.5,font.lab=2,cex.axis=1.1,xlim=c(0,r[2]+50),ylim=c(0,0.005),col="lightcyan3",from=0,to=1000,freq=FALSE,breaks=5,main=paste("Histograma e curva densidade de probabilidade \n no mês de ",meses[k]," (1981-2010): Região ",as.character(i),sep=""),ylab="",xlab="Precipitação acumulada (mm) ") box() lines(dnormal,lwd=2,col="red",xlim=r,xlab="Precipitação acumulada (mm)",ylab="Densidade") grid(col=80,lty=3) axis(side=2,tck=-.015,labels=NA) mtext("Densidade de probabilidade",font=2,cex=1.4,side=2,las=3,line=4) legend("topright",cex=1.1,xpd=TRUE,legend=c("Curva ajustada (função Gama)",paste("Empírica",sep="")),lty=1,col=c(2,1),lwd=2,bty="o",ncol=1) mtext(side=1,cex=1,adj=1,padj=4,text=" (Produto:CPTEC/INPE)") mtext(at=(0.004),side=4,col="dimgray",text=c(paste(" PARÂMETROS \n DA SÉRIE (mm)","\n\n Média: ",sprintf("%.1f",parametros$media),"\n Mediana: ",sprintf("%.1f",parametros$mediana), "\n Desvio \n padrão: ",sprintf("%.1f",parametros$dp)," \n Mín: ",sprintf("%.1f",parametros$min)," \n Máx: ",sprintf("%.1f",parametros$max))),cex=1) mtext(at=(0.004/1.8),side=4,col="dimgray",text=c(paste(" PERCENTIS(mm)","\n\n 05% - ",sprintf("%.1f",quantil5),"\n 15% - ",sprintf("%.1f",quantil15), "\n 33% - ",sprintf("%.1f",quantil33)," \n 66% - ",sprintf("%.1f",quantil66),"\n 85% - ",sprintf("%.1f",quantil85),"\n 95% - ",sprintf("%.1f",quantil95), "\n\n Curtose: ",sprintf("%.1f",curtose),"\n Assimetria: ",sprintf("%.1f",assimetria))),cex=1) vetor<-(dnormal$x) index <- (dnormal$y)[vetor > quantil5-0.8 & vetor < quantil5+0.8] if(index>0){segments(x0=quantil5,y0=0,x1=quantil5,y1=index[1],lty=3,lwd=2,col="dimgray") text(x=quantil5,y=-0.0001,font=2,cex=1,col="dimgray",label="05%")} index <- (dnormal$y)[vetor > quantil15-0.8 & vetor < quantil15+0.8] if(index>0){segments(x0=quantil15,y0=0,x1=quantil15,y1=index[1],lty=3,lwd=2,col="dimgray") text(x=quantil15+5,y=-0.0001,cex=1,col="dimgray",font=2,label="15%")} index <- (dnormal$y)[vetor > quantil33-0.8 & vetor < quantil33+0.8] if(index>0){segments(x0=quantil33,y0=0,x1=quantil33,y1=index[1],lty=3,lwd=2,col="dimgray") text(x=quantil33+5,y=-0.0001,cex=1,col="dimgray",font=2,label="33%")} index <- (dnormal$y)[vetor > quantil66-0.8 & vetor < quantil66+0.8] if(index>0){segments(x0=quantil66,y0=0,x1=quantil66,y1=index[1],lty=3,lwd=2,col="dimgray") text(x=quantil66+5,y=-0.0001,cex=1,col="dimgray",font=2,label="66%")} index <- (dnormal$y)[vetor > quantil85-0.8 & vetor < quantil85+0.8] if(index>0){segments(x0=quantil85,y0=0,x1=quantil85,y1=index[1],lty=3,lwd=2,col="dimgray") text(x=quantil85+5,y=-0.0001,cex=1,col="dimgray",font=2,label="85%")} index <- (dnormal$y)[vetor > quantil95-0.8 & vetor < quantil95+0.8] if(index>0){segments(x0=quantil95,y0=0,x1=quantil95,y1=index[1],lty=3,lwd=2,col="dimgray") text(x=quantil95+7,y=-0.00008,cex=1,col="dimgray",font=2,label="95%")} print(index) if(index == 0){ index<-0.008678319 } dev.off() }}