Use Data DF21

library(data.table)
#library(lubridate)
library(tidyverse) 
library(knitr)
library(readr)

Upload data

df21 <- read_csv("D:/FREELANCE-USA/Additional resources for project/dataset/DF21.csv")
## 
## -- Column specification --------------------------------------------------------
## cols(
##   .default = col_double()
## )
## i Use `spec()` for the full column specifications.
dt <- na.omit(df21)

Data preparation and cleaning variable

dt$Cool.Demand.Status <- NULL
dt$SurfaceAzimuthAngle_south <- NULL 
dt$SurfaceAzimuthAngle_north <- NULL 
dt$SurfaceAzimuthAngle_west <- NULL
dt$SurfaceAzimuthAngle_east <- NULL
dim(dt)
## [1] 14125    28
dt$solar_incidence_angle_north<- NULL
dt$solar_incidence_angle_south<- NULL
dt$solar_incidence_angle_east<- NULL

dt$beam_solar_radiation_south<- NULL
dt$beam_solar_radiation_east<- NULL
dt$beam_solar_radiation_north<- NULL

dt$diffuse_solar_radiation_south<- NULL
dt$diffuse_solar_radiation_north<- NULL
dt$diffuse_solar_radiation_east<- NULL

str(dt)
## tibble [14,125 x 19] (S3: tbl_df/tbl/data.frame)
##  $ Year                        : num [1:14125] 2018 2018 2018 2018 2018 ...
##  $ Month                       : num [1:14125] 6 6 6 6 6 6 6 6 6 6 ...
##  $ Day                         : num [1:14125] 5 5 5 5 5 5 5 5 5 5 ...
##  $ Hour                        : num [1:14125] 14 14 14 14 14 14 14 14 14 14 ...
##  $ Minutes                     : num [1:14125] 40 42 44 46 48 50 52 54 56 58 ...
##  $ Temp.Outdoor.Avg            : num [1:14125] 67 67 67 67 67 67 67 67 67 68 ...
##  $ Temp.Dewpoint.Avg           : num [1:14125] 52 52 52 52 52 52 52 52 52 54 ...
##  $ Temp                        : num [1:14125] 72.5 72.5 72.5 72.5 72.5 72.5 72.5 72.5 72.5 72.5 ...
##  $ Humidity                    : num [1:14125] 49 49 49 49 49 49 49 49 49 49 ...
##  $ Cool.Setpoint               : num [1:14125] 73 72 72 72 72 72 72 72 72 72 ...
##  $ cloud.cover                 : num [1:14125] 89 89 89 89 89 89 89 89 89 89 ...
##  $ GHI                         : num [1:14125] 0.496 0.494 0.491 0.489 0.486 0.483 0.481 0.478 0.475 0.472 ...
##  $ SolarTime                   : num [1:14125] 14.1 14.1 14.2 14.2 14.2 ...
##  $ declination_angle           : num [1:14125] 22.5 22.5 22.5 22.5 22.5 ...
##  $ solar_altitude_angle        : num [1:14125] 58.4 58.1 57.7 57.4 57 ...
##  $ solar_azimuth_angle         : num [1:14125] 66.3 66.9 67.5 68 68.5 ...
##  $ solar_incidence_angle_west  : num [1:14125] 61.3 60.9 60.4 60 59.5 ...
##  $ beam_solar_radiation_west   : num [1:14125] 0.21 0.213 0.215 0.218 0.221 0.223 0.226 0.228 0.23 0.233 ...
##  $ diffuse_solar_radiation_west: num [1:14125] 0.123 0.122 0.122 0.122 0.121 0.121 0.12 0.12 0.119 0.118 ...
##  - attr(*, "na.action")= 'omit' Named int [1:6] 4590 6026 6743 8180 11048 13920
##   ..- attr(*, "names")= chr [1:6] "4590" "6026" "6743" "8180" ...

Convert date variable in one colom

library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:data.table':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
date <- with(dt, ymd_h(paste(Year, Month, Day, Hour, sep= ' ')))
date <- data.frame(date)  
dn <- cbind(dt,date)
             
dnn <- dn %>% select("date","Temp","Cool.Setpoint")
LongFormatReadings <- function(building){
 
  colnames(building) <- c("date","temp","relh")  
  return(building)
}
west <- LongFormatReadings(dnn)
head(west,3)
##                  date temp relh
## 1 2018-06-05 14:00:00 72.5   73
## 2 2018-06-05 14:00:00 72.5   72
## 3 2018-06-05 14:00:00 72.5   72

Setting parameter

AddOtherVariableAssumptions <- function(building){
 building$radtemp <- building$temp
  building$vel <- 0.1
  building$met <- 1.1
  building$clo <- 0.5
  building$wme <- 0
  return(building)
}

west <- AddOtherVariableAssumptions(west)
head(west,3)
##                  date temp relh radtemp vel met clo wme
## 1 2018-06-05 14:00:00 72.5   73    72.5 0.1 1.1 0.5   0
## 2 2018-06-05 14:00:00 72.5   72    72.5 0.1 1.1 0.5   0
## 3 2018-06-05 14:00:00 72.5   72    72.5 0.1 1.1 0.5   0

TCA Termal Comfort analysis

#Load Google's Javscript Engine V8 (See https://cran.r-project.org/web/packages/V8/vignettes/v8_intro.html) 
library(V8)
## Using V8 engine 6.2.414.50
#Create a new context
ct <- v8()

#Load Javascript Library for forEach function
ct$source(system.file("js/underscore.js", package="V8"))

ct$source("D:/PMV/thermal_comfort-master/pmv/comfortmodels.js") # change with your directory
## [1] "function (ta, tr, runningMean, vel) {\n    var to = (ta + tr) / 2;\n    var coolingEffect = 0;\n    if (vel >= 0.2 && to > 25) {\n        // calculate cooling effect of elevated air speed\n        // when top > 25 degC.\n        var coolingEffect = 1.7856 * Math.log(vel) + 2.9835;\n    }\n    var tComf = 0.33 * runningMean + 18.8;\n    if(runningMean > 15){\n        var tComfILower = tComf - 2;\n        var tComfIUpper = tComf + 2 + coolingEffect;\n        var tComfIILower = tComf - 3;\n        var tComfIIUpper = tComf + 3 + coolingEffect;\n        var tComfIIILower = tComf - 4;\n        var tComfIIIUpper = tComf + 4 + coolingEffect;\n    } else if (12.73 < runningMean && runningMean < 15){\n        var tComfLow = 0.33 * 15 + 18.8;\n        var tComfILower = tComfLow - 2;\n        var tComfIUpper = tComf + 2 + coolingEffect;\n        var tComfIILower = tComfLow - 3;\n        var tComfIIUpper = tComf + 3 + coolingEffect;\n        var tComfIIILower = tComfLow - 4;\n        var tComfIIIUpper = tComf + 4 + coolingEffect;\n    } else {\n        var tComfLow = 0.33 * 15 + 18.8;\n        var tComfILower = tComfLow - 2;\n        var tComfIUpper = tComf + 2;\n        var tComfIILower = tComfLow - 3;\n        var tComfIIUpper = tComf + 3 + coolingEffect;\n        var tComfIIILower = tComfLow - 4;\n        var tComfIIIUpper = tComf + 4 + coolingEffect;\n    }\n    var acceptabilityI, acceptabilityII, acceptabilityIII;\n\n    if (comf.between(to, tComfILower, tComfIUpper)) {\n        // compliance at all levels\n        acceptabilityI = acceptabilityII = acceptabilityIII = true;\n    } else if (comf.between(to, tComfIILower, tComfIIUpper)) {\n        // compliance at II and III only\n        acceptabilityII = acceptabilityIII = true;\n        acceptabilityI = false;\n    } else if (comf.between(to, tComfIIILower, tComfIIIUpper)) {\n        // compliance at III only\n        acceptabilityIII = true;\n        acceptabilityI = acceptabilityII = false;\n    } else {\n        // neither\n        acceptabilityI = acceptabilityII = acceptabilityIII = false;\n    }\n    r = {}\n    r.acceptabilityI = acceptabilityI\n    r.acceptabilityII = acceptabilityII\n    r.acceptabilityIII = acceptabilityIII\n    r.tComfILower = tComfILower\n    r.tComfIILower = tComfIILower\n    r.tComfIIILower = tComfIIILower\n    r.tComfIUpper = tComfIUpper\n    r.tComfIIUpper = tComfIIUpper\n    r.tComfIIIUpper = tComfIIIUpper\n    return r\n    return [[acceptabilityIII, tComfIIILower, tComfIIIUpper],\n            [acceptabilityII, tComfIILower, tComfIIUpper],\n            [acceptabilityI, tComfILower, tComfIUpper]];\n}"
ct$source("D:/PMV/thermal_comfort-master/pmv/util.js") # change with your directory
## [1] "function (x) {\n    return (x - 32) * 5 / 9;\n}"
ct$source("D:/PMV/thermal_comfort-master/pmv/psychrometrics.js") # change with your directory
## [1] "function (ta, vel, tglobe, diameter, emissivity) {\n    pow = Math.pow;\n    return pow(pow(tglobe + 273, 4) + (1.1 * pow(10, 8) * pow(vel, 0.6)) / (emissivity * pow(diameter, 0.4)) * (tglobe - ta), 0.25) - 273;\n}"

Apply the function over all the table will produce PMV table

west.elevated.air.pmv <- data.table(ct$call("_.map",
west, JS("function(x){return(comf.pmvElevatedAirspeed(util.FtoC(x.temp),
         util.FtoC(x.radtemp),x.vel/196.85,x.relh,x.clo,x.met,x.wme))}")))

west.pmv<- data.table(ct$call("_.map", west, JS("function(x)
  {return(comf.pmv(util.FtoC(x.temp),util.FtoC(x.radtemp),x.vel/196.85,x.relh,x.clo,x.met,x.wme))}")))

west.pmv
##              pmv      ppd
##     1: -3.764946 99.99487
##     2: -3.777924 99.99542
##     3: -3.777924 99.99542
##     4: -3.777924 99.99542
##     5: -3.777924 99.99542
##    ---                   
## 14121: -3.472400 99.94757
## 14122: -3.472400 99.94757
## 14123: -3.472400 99.94757
## 14124: -3.625261 99.98345
## 14125: -3.625261 99.98345

Comfort threshold

west.elevated.air.pmv$iscomfortable <- ifelse(west.elevated.air.pmv$pmv > -3.764946 & 
               west.elevated.air.pmv$pmv < 3.764946,"comfortable","uncomfortable")
kable(west.elevated.air.pmv[,.(frequency=length(pmv)),by="iscomfortable"])
iscomfortable frequency
uncomfortable 5944
comfortable 8181

Plot confort zone

library(dygraphs)
west.elevated.air.pmv$date <- west$date
dygraph(west.elevated.air.pmv[,.(date,pmv)], main = "Comfort Level") %>% dyRangeSelector()