Temperature_Plot <- function (infile, title) { source("Utilities.R") data = read.table(infile,header=T) attach(data) minall = min(Tmin, reclo_T, avglo_T, minhi_T)-10 maxall = max(Tmax, rechi_T, avghi_T, maxlo_T) RDate = strptime(as.character(Date),format="%Y/%m/%d") plot(RDate, Tmax, ylim=c(minall, maxall), type='n', xlab = "Date", ylab="Temperature in degrees F", main=paste("Daily Temperature Extrema for",title), xaxt='n', tck=1) # polygons first # Create new avghi and avglo which include intersections with # the data values. relDate = RDate - min(RDate) max_hi = intersect.lines(Tmax, avghi_T, relDate) max_lo = intersect.lines(Tmax, avglo_T, relDate) min_hi = intersect.lines(Tmin, avghi_T, relDate) min_lo = intersect.lines(Tmin, avglo_T, relDate) rec_hi = intersect.lines(Tmax, rechi_T, relDate) rec_lo = intersect.lines(Tmin, reclo_T, relDate) # low < loavg and low < freezing o1 = order(c(relDate,min_lo[1,])); ndate = as.POSIXct(min_lo[1,] + min(RDate)) n1RDate = c(RDate, ndate)[o1] nTmin = c(Tmin, min_lo[2,])[o1] n1avglo_T = c(avglo_T, min_lo[2,])[o1] ndate = as.POSIXct(max_lo[1,] + min(RDate)) o2 = rev(order(c(relDate,max_lo[1,]))); n2RDate = c(RDate, ndate)[o2] nTmax = c(Tmax, max_lo[2,])[o2] n2avglo_T = c(avglo_T, max_lo[2,])[o2] doubleDate = c(n1RDate,n2RDate) polygon(doubleDate, c(pmin(nTmin,n1avglo_T),pmin(nTmax,n2avglo_T)), col="blue", border=NA) # below freezing n0 = intersect.lines(Tmin, rep(32,length(Tmin)), relDate) ndate = as.POSIXct(n0[1,] + min(RDate)) o1 = order(c(relDate,n0[1,])); n0RDate = c(RDate, ndate)[o1] nTmin = c(Tmin, n0[2,])[o1] revDate = strptime(rev(strftime(n0RDate)),"%Y-%m-%d") doubleDate = c(n0RDate,revDate) polygon(doubleDate, c(pmin(nTmin,32),rep(32,length(nTmin))), col="blue", border=NA ) # loavg < low < min(hi, hiavg) newpts = c(min_lo[2,], max_lo[2,], min_hi[2,]) o1 = order(c(relDate, min_lo[1,], max_lo[1,], min_hi[1,])); ndate = c(as.POSIXct(min_lo[1,] + min(RDate)), as.POSIXct(max_lo[1,] + min(RDate)) , as.POSIXct(min_hi[1,] + min(RDate))) n1RDate = c(RDate, ndate)[o1] n1Tmin = c(Tmin, newpts)[o1] n1Tmax = c(Tmax, newpts)[o1] n1avglo_T = c(avglo_T, newpts)[o1] n1avghi_T = c(avghi_T, newpts)[o1] newpts = c(max_hi[2,], max_lo[2,], min_hi[2,]) o2 = rev(order(c(relDate,max_hi[1,], max_lo[1,], min_hi[1,]))); ndate = c(as.POSIXct(max_hi[1,] + min(RDate)), as.POSIXct(max_lo[1,] + min(RDate)) , as.POSIXct(min_hi[1,] + min(RDate))) n2RDate = c(RDate, ndate)[o2] n2Tmin = c(Tmin, newpts)[o1] n2Tmax = c(Tmax, newpts)[o2] n2avghi_T = c(avghi_T, newpts)[o2] doubleDate = c(n1RDate,n2RDate) polygon(doubleDate, c(pmin(pmax(n1Tmin,n1avglo_T),n1Tmax,n1avghi_T),pmin(n2Tmax,n2avghi_T)), col="springgreen", border=NA) # hi > hiavg newpts = min_hi[2,] o1 = order(c(relDate, min_hi[1,])); ndate = as.POSIXct(min_hi[1,] + min(RDate)) n1RDate = c(RDate, ndate)[o1] n1Tmin = c(Tmin, newpts)[o1] n1avghi_T = c(avghi_T, newpts)[o1] newpts = max_hi[2,] o2 = rev(order(c(relDate,max_hi[1,]))); ndate = c(as.POSIXct(max_hi[1,] + min(RDate))) n2RDate = c(RDate, ndate)[o2] n2Tmax = c(Tmax, newpts)[o2] n2avghi_T = c(avghi_T, newpts)[o2] doubleDate = c(n1RDate,n2RDate) polygon(doubleDate, c(pmax(n1Tmin,n1avghi_T),pmax(n2Tmax,n2avghi_T)), col="red", border=NA) # hi > rechi n1RDate = RDate n1Tmin = rechi_T newpts = rec_hi[2,] o2 = rev(order(c(relDate,rec_hi[1,]))); ndate = c(as.POSIXct(rec_hi[1,] + min(RDate))) n2RDate = c(RDate, ndate)[o2] n2Tmax = c(Tmax, newpts)[o2] n2rechi_T = c(rechi_T, newpts)[o2] doubleDate = c(n1RDate,n2RDate) polygon(doubleDate, c(n1Tmin,pmax(n2Tmax,n2rechi_T)), col="orange", border=NA) # lo < reclo n1RDate = RDate n1Tmin = reclo_T newpts = rec_lo[2,] o2 = rev(order(c(relDate,rec_lo[1,]))); ndate = c(as.POSIXct(rec_lo[1,] + min(RDate))) n2RDate = c(RDate, ndate)[o2] n2Tmin = c(Tmin, newpts)[o2] n2reclo_T = c(reclo_T, newpts)[o2] doubleDate = c(n1RDate,n2RDate) polygon(doubleDate, c(n1Tmin,pmin(n2Tmin,n2reclo_T)), col="DeepSkyBlue", border=NA) # Lines # hot ones points(RDate, Tmax, col="black",type='l', pch=20) points(RDate, rechi_T, col="red", type='l', lty='solid', lwd=1) points(RDate, avghi_T, col="red", type='l', lty='solid', lwd=2) #points(RDate, minhi_T, col="red", type='l', lty='solid', lwd=1) #points(RDate, rechi_T, col="red", type='l', lty='dashed', lwd=2) #points(RDate, avghi_T, col="red", type='l', lty='solid', lwd=2) #points(RDate, minhi_T, col="red", type='l', lty='dotted', lwd=2) # cold ones points(RDate, Tmin, col="black", type='l', pch=20) #points(RDate, reclo_T, col="blue", type='l', lty='dashed', lwd=2) #points(RDate, avglo_T, col="blue", type='l', lty='solid', lwd=2) #points(RDate, maxlo_T, col="blue", type='l', lty='dotted', lwd=2) points(RDate, reclo_T, col="blue", type='l', lty='solid', lwd=1) points(RDate, avglo_T, col="blue", type='l', lty='solid', lwd=2) #points(RDate, maxlo_T, col="blue", type='l', lty='solid', lwd=1) if (minall < 32) { lines(c(RDate[1],RDate[length(relDate)]), c(32,32), col="blue", type='l', lty='solid') } # add text for major extrema and broken records text(RDate[which(Tmax==max(Tmax))[1]], max(Tmax)+2, labels=paste("Max =",max(Tmax)),xpd=T,font=2) text(RDate[which(Tmin==min(Tmin))[1]], min(Tmin)-2, labels=paste("Min =",min(Tmin)),xpd=T,font=2) # legend legend(RDate[1], minall+3, horiz=T, legend=c("Rec Lo", "Avg Lo", "Avg Hi", "Rec Hi"), col=c("blue","blue","red","red"), lty=c("solid","solid","solid","solid"), lwd=c(1,2,2,1)) #legend(RDate[1], minall+5, horiz=T, legend=c("Rec Lo", "Min Hi", "Avg Lo", "Avg Hi", "Max Lo", "Rec Hi"), col=c("blue","blue","blue","red","red","red"), lty=c("dashed","dotted","solid","solid","dotted","dashed"), lwd=c(2,2,2,2,2,2)) # build horizontal axis mon = c('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec') if (max(RDate)-min(RDate) > 120) { # build a monthly axis - put month name on the 15th of each month xlocs = RDate[RDate$mday==15] xlabs = paste(mon[xlocs$mon+1], ifelse(xlocs$mon==0, xlocs$year+1900,"") ) xlocs = as.numeric(as.POSIXct(xlocs)) axis(1, at=xlocs, labels=xlabs) } }