library(ggplot2) library(grid) library(gridExtra) library(comf) #weather="HND_Catacamas.787140_SWERA.epw" #weather="CHN_Hong.Kong.SAR.450070_CityUHK.epw" #weather='RUS_Moscow.276120_IWEC.epw' weather = "~/weather-datasets/GBR_London-Heathrow.AP.1989_2020High90pct_CIBSETM49.epw" epw = read.csv(weather, skip=9) colnames(epw) = cbind('year','month','day','hour','minute','dataSource','dryBulb','dewPoint','relHum','atmosPressure','exthorrad','extdirrad','horirsky','glohorrad','dirnorrad','difhorrad','glohorillum','dirnorillum','difhorillum','zenlum','winddir','windspd','totskycvr','opaqskycvr','visibility','ceiling_hgt','presweathobs','presweathcodes','precip_wtr','aerosol_opt_depth','snowdepth','days_last_snow','col1','col2','col3') #epw$dateTime <- paste(epw$month,epw$day,epw$hour) #epw$dataTime = as.POSIXlt(epw$dateTime,format="%m %d %H") for (i in c(90,0,-90)){ temps = data.frame(month =aggregate(dryBulb ~ hour + month,epw,mean)$month, hour =aggregate(dryBulb ~ hour + month,epw,mean)$hour, dryBulb_mean=aggregate(dryBulb ~ hour + month,epw,mean)$dryBulb, difhorrad =aggregate(difhorrad ~ hour + month,epw,mean)$difhorrad, dirnorrad =aggregate(dirnorrad ~ hour + month,epw,mean)$dirnorrad) # doy =as.numeric(format(ISOdate(2000, temps$month, 15), "%j")), # dayAng =NISTdegTOradian((360)*((temps$doy-1)/365)), # EOT =2.2918*(0.00075+0.1868*cos(temps$dayAng)-3.2077*sin(temps$dayAng)-1.4615*cos(2*temps$dayAng)-4.089*sin(2*temps$dayAng)), # solDec =NISTdegTOradian(23.45*sin((NISTdegTOradian(360*((temps$doy+284)/365))))), # AST =temps$hour+(temps$EOT/60)+((0-(15*0))/15), # hrAng =NISTdegTOradian(15*(temps$AST-12)), # solAlt =asin((sin(Lat))*sin(temps$solDec)+cos(Lat)*cos(temps$solDec)*cos(temps$hrAng)), # solAz =asin(sin(temps$hrAng)*(cos(temps$solDec)/cos(temps$solAlt))), # ssAZ =NISTdegTOradian(0)-temps$solAz, # Etb =temps$dirnorrad*(cos(temps$solAlt)*cos(temps$ssAZ))) } temps$number = as.numeric(format(ISOdate(2000, temps$month, 15), "%j")) months = c(1:12) E <- read.csv("edges.csv", header=TRUE) baseArea=abs(polyarea(E$V1x,E$V1y)) HEIGHT=2.8 WIN_U=1 WALL_U=1.0 WIN_G=0.5 diff_G=0.5 E$LENGTH=(((E$V2x-E$V1x)^2)+((E$V2y-E$V1y)^2))^(1/2) E$AREA=(E$LENGTH)*HEIGHT E$WIN_AREA=(E$AREA)*(E$WIN_frac) E$WALL_AREA=(E$AREA)*(1-E$WIN_fra) WIN_UA=(sum(E$WIN_AREA))*WIN_U WALL_UA=(sum(E$WALL_AREA))*WALL_U BASE_UA=(baseArea*2)*WALL_U Htr=(WIN_UA)+(WALL_UA)+(BASE_UA) E$Azimuth=ifelse( E$Or==180,180, ifelse (E$Or==0,0, ifelse(E$Or==-90,-(atan(((-(E$V2x-E$V1x))-((E$V2x-E$V1x)))/(((E$V2y-E$V1y))-(-(E$V2y-E$V1y))))*(180/pi)+90), ifelse (E$Or==90,(180-(atan(((-(E$V2x-E$V1x))-((E$V2x-E$V1x)))/(((E$V2y-E$V1y))-(-(E$V2y-E$V1y))))*(180/pi)+90)), 0)))) L=NISTdegTOradian(DATA$V7) hours=c(0:23) for (j in months){ for (i in hours) { n=as.numeric(j)*29.6 Rdc=NISTdegTOradian((360)*((n-1)/365)) ET=2.2918*(0.00075+0.1868*cos(Rdc)-3.2077*sin(Rdc)-1.4615*cos(2*Rdc)-4.089*sin(2*Rdc)) DEC=NISTdegTOradian(23.45*sin((NISTdegTOradian(360*((n+284)/365))))) # solar declination AST=i+(ET/60)+((DATA$V8-(15*DATA$V9))/15) # apparent solar time H=NISTdegTOradian(15*(AST-12)) # hour angle B=asin((sin(L))*sin(DEC)+cos(L)*cos(DEC)*cos(H)) # solar altitude angle AZ=asin(sin(H)*(cos(DEC)/cos(B))) ssAZ=NISTdegTOradian(E$Azimuth)-AZ solar=mean(subset(subset(EPW, HR==i), format(DT,'%m')==j)$V15) diffuse=E$WIN_AREA*(mean(subset(subset(EPW, HR==i), format(DT,'%m')==j)$V16)) Etb=solar*(cos(B)*cos(ssAZ)) Etb[Etb<0]=0 SOL=Etb*E$WIN_AREA*WIN_G assign(paste0("S",j,"",i), (SOL+(diff_G*diffuse))/Htr) } } plots <- list() for (i in months) { tsix=temps[temps$month == i, ]$dryBulb_mean[6] if (tsix < -5){ clo = 1.0 } else { if(tsix>-55 && tsix<5){ clo = 0.818-(0.0364*tsix) } else { if(tsix>5 && tsix<26){ clo = 10^(-0.1635-(0.0066*tsix)) } else { clo =0.46 } } } results=data.frame(temps=runif(100, 15, 35)) results$pmv=sapply(results$temps, function(x) abs(calcPMV(x,x,0.1,50,clo,1.5,0,58.15)) ) comftemp=results$temps[which.min(results$pmv)] plots[[i]] <- ggplot() + geom_line(data=temps[temps$month == i, ],aes(x=hour,y = dryBulb_mean)) + # geom_line(data=temps[temps$month == i, ],aes(x=hour,y = dryBulb_max)) + # geom_line(data=temps[temps$month == i, ],aes(x=hour,y = dryBulb_min)) + scale_y_continuous(limits = c(-10, 40),breaks=seq(from = -10, to = 40, by = 5),expand = c(0, 0)) + scale_x_continuous(breaks=c(0,6,12,18,24),expand = c(0, 0))+ geom_hline(yintercept=comftemp)+ geom_ribbon(data=temps,aes(ymin=18,ymax=22,x=hour,alpha=0.3,fill='blue'))+ geom_ribbon(data=temps,aes(ymin=22,ymax=24,x=hour,alpha=0.3,fill='red'))+ theme( panel.background = element_rect(fill = 'white'), panel.border = element_rect(colour = "black", fill=NA, size=1), panel.grid.major = element_line(colour = "#B0B0B0"), panel.grid.minor.x = element_blank(), panel.grid.minor.y = element_blank())+ ggtitle(format(ISOdate(2000, i, 1), "%B")) + labs(x='Date', y= expression("Temperature (" *degree * "C)")) # axis.title.x = element_text(vjust=-3,size=10), # axis.title.y = element_text(vjust=+3,size=10), # panel.border = element_rect(fill=NA,colour = "#AAAAAA",size=2), # plot.margin=unit(c(1,1,1.5,1.2),"cm") } for (i in months) { pdf(paste("/home/martin/www/month_", i, ".pdf", sep=""),width = 3, height = 5) print(plots[[i]]) dev.off() }