This blog is run by Jason Benedict and Doug Beare who are based in Penang, Malaysia to share insights and developments on open source software that can be used to analyze patterns and trends in weather and climate data; with particular emphasis on developing countries. The work featured here was undertaken as part of the CGIAR Research Program on Climate Change, Agriculture and Food Security (CCAFS) - http://ccafs.cgiar.org

However, the views presented here are of the authors only, and do not reflect the official position of CCAFS, CGIAR, its members, partners or donors.

Friday, 8 August 2014

Blowing in the wind - How many times must we keep experiencing the haze before the Governments do something about it?

Key points of post

  • Visibility in Penang is affected by the seasonal monsoon winds
  • Visibility has got worse since 2001.
  • The governments of south-east Asia must urgently solve this potentially serious health issue.
Every year South-East Asians suffer from the effects of fires burning in the forests of Indonesia. The more intense fires are lit deliberately to clear land for palm oil plantations. It is illegal to use fire to clear land in Indonesia, but the practice is unfortunately widespread, leading to habitat destruction and loss of biodiversity. A satellite image captured by the Moderate Resolution Imaging Spectroradiomater (MODIS) from NASA below for March 2014 illustrates when the fires were particularly bad (fires shown in red). 

Source: NASA Earth Observatory (2014) - http://1.usa.gov/1kqfdYA

Most fires occur during the dry season between April and October. The fires in 2014 started in February, and actually forced Indonesia to declare a state of emergency because of poor air quality. 

But what about the situation for us in Penang ? Our island is subject to two seasonal monsoon winds: the north-east and the south-west. These patterns are clear from the wind-rose diagram below which shows the average monthly wind-speeds and directions for the period 1949-2014. 


Given the location of Indonesia relative to Penang, it follows that visibility will be worst during the south-west monsoon which blows hardest, and most consistently, between May and August each year. Conversely we might expect visibility to be best during the north-east monsoon. This is indeed exactly what we observe. Average visibility readings from Bayan Lepas International Airport by week and year 2001-2014 are plotted below. The white squares represent good visibility and cleaner air; the dark squares poor visibility and high levels of air pollution.


It is clear that visibility has gotten worse since 2001, ie. the numbers of white squares have decreased. The pattern is also strongly seasonal which reflects: (a) the monsoons; and (b) the fact that fires are used more in Indonesia’s dry season. The relatively long periods of good visibility (>10 km's) that Penangites must have enjoyed especially during January, November and December seem now to be a thing of the past. I’m (Doug) a keen hiker regularly gazing down on the Island from the heights of Tea Station 5, 89, and Penang Hill. I’ve only been here since January 2012 and, judging from these data, I guess I've probably never seen a really good view!

Let’s hope the Government of Malaysia / Penang will continue to apply pressure and work with the Indonesians and other countries in the region to solve this terrible problem.

If you're keen on producing similar plots like the ones we've included in this blog post, you can use the code provided below. The wind-rose plot was produced using the 'openair' package for R which is an open-source air pollution analysis tool. There are a number of useful plots that can be produced using this package to show magnitude and direction of wind / air pollutants which you can further explore on your own. We also provide you with the code to produce the 'heatmap' plot of the visibility data which was pulled into R from Weather Underground (using the 'weatherData' package).  

# Load libraries
 
library(ggplot2)
library(RColorBrewer)
library(openair)
library(weatherData)
library(mgcv)
library(scales)
library(plyr)
library(reshape2)
library(circular)
library(gridExtra)
library(lubridate)
library(weathermetrics)
library(zoo)
 
# Setting work directory
 
setwd("d:\\ClimData")
 
# Reading and reformatting raw daily data downloaded from NCDC
 
dat<-read.table("2422706962434dat.txt",header=TRUE,fill=TRUE,na.strings=c("*","**","***","****","*****","******","0.00T*****"))
 
colnames(dat)<-tolower(colnames(dat))
 
Sys.setenv(TZ = "UTC")
dat$date <- as.POSIXct(strptime(dat$yr..modahrmn,format="%Y%m%d%H%M"))  + 8 * 60 * 60
 
dat$dir[dat$dir == 990.0] <- NA 
 
# Convert windspeed in miles/hour to metres/second
 
dat$wspd <- (dat$spd)*0.44704
 
# Create columns for year, month and dat
 
dat$year <- as.numeric(format(dat$date,"%Y"))
 
dat$month <- as.numeric(format(dat$date,"%m"))
 
dat$day <- as.numeric(format(dat$date,"%d"))
 
# Subset year of interest
 
datsub <- subset(dat, year >= 1949 )
 
# Remove NA's
 
datsub$dir[is.na(datsub$dir)] <- 0
 
# Plot and produce png file of Wind Rose
 
png(filename = "Monthly_WindRose_Penang_1949-2014.png",height=10,width=10,
    bg = "white",units='in', res = 600, family = "", restoreConsole = TRUE,
    type = "cairo-png")
 
windRose(datsub,ws="wspd",wd="dir",bias.corr = TRUE,border=TRUE,type = "month",width = 0.5,grid.line = 10,
         statistic = "prop.mean",offset = 10, paddle =FALSE,cols="hue",annotate=FALSE,auto.text=FALSE,
         ws.int=2,breaks=c(0,2,4,6,8,10,12),key = TRUE, key.footer = "(meter/second)", key.position = "bottom",na.action=NULL,
         key.header = "Wind Speed",main="       Average Monthly Wind Rose plots for Penang (Bayan Station)\n from years 1949 - 2014 \n",font.main=1,
         sub="\nData source: Integrated Surface Database (ISD) - National Climatic Data Centre (NCDC)",font.sub=2)
 
dev.off()
 
###########################################
#### Produce visibility 'heatmap' plot ####
###########################################
 
### Find required 4 letter station code using function below
 
# Example getStationCode("Penang") or getStationCode("Honiara")
 
# Information can also be acquired from the following link - http://weather.rap.ucar.edu/surface/stations.txt
 
getStationCode("Butterworth")
 
# Use station code below to get required plot data and parameters
 
station.id="WMKP"
 
s<-getStationCode(station.id)
 
s1<-(strsplit(s,split= " "))[[1]]
 
station.name<-paste(s1[4],s1[5])
 
### Getting summarized weather data for WU
 
ws2001<-getSummarizedWeather(station.id, "2001-01-01", "2001-12-31", opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
ws2002<-getSummarizedWeather(station.id, "2002-01-01", "2002-12-31", opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
ws2003<-getSummarizedWeather(station.id, "2003-01-01", "2003-12-31", opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
ws2004<-getSummarizedWeather(station.id, "2004-01-01", "2004-12-31", opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
ws2005<-getSummarizedWeather(station.id, "2005-01-01", "2005-12-31", opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
ws2006<-getSummarizedWeather(station.id, "2006-01-01", "2006-12-31", opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
ws2007<-getSummarizedWeather(station.id, "2007-01-01", "2007-12-31", opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
ws2008<-getSummarizedWeather(station.id, "2008-01-01", "2008-12-31", opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
ws2009<-getSummarizedWeather(station.id, "2009-01-01", "2009-12-31", opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
ws2010<-getSummarizedWeather(station.id, "2010-01-01", "2010-12-31", opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
ws2011<-getSummarizedWeather(station.id, "2011-01-01", "2011-12-31", opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
ws2012<-getSummarizedWeather(station.id, "2012-01-01", "2012-12-31", opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
ws2013<-getSummarizedWeather(station.id, "2013-01-01", "2013-12-31", opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
ws2014<-getSummarizedWeather(station.id, "2014-01-01", Sys.Date(), opt_custom_columns=T,
custom_columns=c(1:23), opt_verbose=T)
 
ws<-rbind(ws2001,ws2002,ws2003,ws2004,ws2005,ws2006,ws2007,ws2008,ws2009,ws2010,ws2011,ws2012,ws2013,ws2014) 
 
summary(ws)
 
str(ws)
 
### List of variables from Wunderground ###
 
# [1] "MYT"                        "Max_TemperatureC"           "Mean_TemperatureC"          "Min_TemperatureC"          
# [5] "Dew_PointC"                 "MeanDew_PointC"             "Min_DewpointC"              "Max_Humidity"              
# [9] "Mean_Humidity"              "Min_Humidity"               "Max_Sea_Level_PressurehPa"  "Mean_Sea_Level_PressurehPa"
# [13] "Min_Sea_Level_PressurehPa"  "Max_VisibilityKm"           "Mean_VisibilityKm"          "Min_VisibilitykM"          
# [17] "Max_Wind_SpeedKm_h"         "Mean_Wind_SpeedKm_h"        "Max_Gust_SpeedKm_h"         "Precipitationmm"           
# [21] "CloudCover"                 "Events"                     "WindDirDegrees"            
 
colnames(ws)<-c("date","date1","maxtemp","meantemp","mintemp","dewp","meandewp","maxdewp","maxhum","meanhum","minhum","maxslp","meanslp",
                "minslp","maxvsb","meanvsb","minvsb","maxwspd","meanwspd","maxgust","prcp","cc","events","wd")
 
## Adding date columns

ws$dates <- as.Date(ws$date)
ws$year <- as.numeric(as.POSIXlt(ws$dates)$year+1900) ws$month <- as.numeric(as.POSIXlt(ws$dates)$mon+1) ws$monthf <- factor(ws$month,levels=as.character(1:12),labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),ordered=TRUE) ws$weekday <- as.POSIXlt(ws$dates)$wday ws$weekdayf <- factor(ws$weekday,levels=rev(0:6),labels=rev(c("Mon","Tue","Wed","Thu","Fri","Sat","Sun")),ordered=TRUE) ws$yearmonth <- as.yearmon(ws$date) ws$yearmonthf <- factor(ws$yearmonth) ws$week <- as.numeric(format(as.Date(ws$dates),"%W")) ws$weekf<- factor(ws$week) ws$jday<-yday(ws$dates)   # Define colour palette   col<-c("black","grey10","grey20","grey50","white")   # Plot 'heatmap'   v <- ggplot(data=ws,aes(x=week,y=year,fill=meanvsb))+ geom_tile(colour="black",size=0.65)+ theme_bw()+ scale_fill_gradientn(colours=col,name="Visibility\n(km)\n")+coord_equal(ratio=1)+ ylab("YEAR\n")+ xlab("\nWEEK OF YEAR\n\nSource: Weather Underground (2014)")+ scale_y_continuous(expand = c(0,0),breaks = seq(2000, 2015, 1)) + scale_x_discrete(expand = c(0,0),breaks = seq(0,52,2))+ ggtitle("Average weekly visibility in Penang\n")+ theme(panel.background=element_rect(fill="transparent"), panel.border=element_blank(), axis.title.y=element_text(size=10,colour="grey20"), axis.title.x=element_text(size=10,colour="grey20"), axis.text.y=element_text(size=10,colour="grey20",face="bold"), axis.text.x=element_text(size=10,colour="grey20",face="bold"), plot.title = element_text(lineheight=1.2, face="bold",size = 14, colour = "grey20"), panel.grid.major = element_blank(), legend.key.width=unit(c(0.2,0.2),"in"))   v   # Save plot to png   ggsave(v, file="Penang_Weekly_Average_Visibility.png", width=15, height=5,dpi=400,type = "cairo-png")
Created by Pretty R at inside-R.org

Monday, 4 August 2014

Will I get soaked? - Best time of day to commute on foot or by bike on Penang Island

Key points of post

  • Chances of getting wet are always highest during March, April, October and November;
  • On average 8am-9am is a particularly wet time of day on Penang Island.
  • The wettest hour is 8pm-9pm; especially in October and November.

Unfortunately not nearly enough of us cycle to work in Penang. If you do, however, what is the best time of day to go to avoid getting drenched, and how does this vary over the year ? 

To answer this question we needed to find rainfall data at sufficient resolution (hourly). We searched through a lot of online databases but eventually it was the citizen scientists who came up with the goods via Weather Underground.

From Weather Underground, Jason managed to download hourly rainfall ‘event’ data for Penang for 2002-2013. These data describe simply whether it was raining or not each hour of the day and these numbers can be summed into frequencies.

Unfortunately we do not how much rain fell per hour and it is difficult to acquire such data since rain gauges are usually set up to record daily quantities of rain. Nevertheless we think that the frequency data recording such rainfall ‘events’ are still instructive about whether or not you will get a soaking. 

The ‘event’ data (2002-2013) are plotted in the circular plot below. These types of plot have been discussed in previous blog post. The plot shows that the chance of rain depends on both the month (see previous blogs) and the time of day. As we’ve seen before, January, February, May, June and July are the driest while October and November are the wettest in Penang. In all months there seems to be slightly more chance of getting wet at 8am so it would be best to come into work either earlier or later. 2am and 2pm also seem to be particularly wet times of day. The worst time of day to commute by bike or on foot, however is 8pm and this is particularly true for the months of March, April, October, and November. Luckily most of us are at home by then!


As yet we have no idea what causes these differences but the data certainly suggest that they exist.

The code used to prepare this is very similar to the ones we produced in our previous post here, although you will need the 'weatherData' package to pull the hourly rainfall 'event' data from Weather Underground into R. The hourly data is particularly huge in this case and I would suggest running it over a stable internet connection overnight or getting the datasets year by year and then 'rbind-ing' the data together.  The R code used to produce the above plot is as follows. 

# Load required libraries
 
library(weatherData)
library(zoo)
library(lubridate)
library(plyr)
library(ggplot2)
library(circular)
library(grid)
library(RColorBrewer)
 
# Download hourly weather data from Weather Underground 
 
we <- getWeatherForDate("WMKP", "2002-01-01","2013-12-31", opt_detailed=T, opt_custom_columns=T, custom_columns=c(11))
 
# Convert to characters
 
we$Events<-as.character(we$Events)
 
# Assign values of '1' to Rain Event and 'NA' to Non-Rain Event
 
we$Events[we$Events == "" ] <- NA
 
we$Events[we$Events == "Rain"] <- 1
 
we$Events[we$Events == "Rain-Thunderstorm"] <- 1
 
we$Events[we$Events == "Thunderstorm"] <- NA
 
# Convert to numeric
 
we$Events<-as.numeric(we$Events)
 
# Create date and time columns
 
we$Dates<-as.POSIXct(we$Time)
 
we$year <- as.numeric(as.POSIXlt(we$Dates)$year+1900)
we$month <- as.numeric(as.POSIXlt(we$Dates)$mon+1)
we$monthf <- factor(we$month,levels=as.character(1:12),labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),ordered=TRUE)
we$weekday <- as.POSIXlt(we$Dates)$wday
we$weekdayf <- factor(we$weekday,levels=rev(0:6),labels=rev(c("Mon","Tue","Wed","Thu","Fri","Sat","Sun")),ordered=TRUE)
we$week <- as.numeric(format(as.Date(we$Dates),"%W"))
we1$hour <- as.numeric(format(strptime(we1$Dates, format = "%Y-%m-%d %H:%M"),format = "%H"))
we1$min <- as.numeric(format(strptime(we1$Dates, format = "%Y-%m-%d %H:%M"),format = "%M"))
 
## Use only data on the hour
 
we1<- subset(we, min == 0)
 
we2 <- ddply(we1,.(monthf,hour),summarize, event = sum(Events,na.rm=T))
 
# Define colour palette
 
col<-brewer.pal(9,"Blues")
 
# Plot circular chart of rain event frequency
 
r1 <-  ggplot(we2, aes(x=monthf, y=hour, fill=event)) +
       geom_tile(colour="grey70") +
       scale_fill_gradientn(colours=col,guide=FALSE)+
       scale_y_continuous(breaks = seq(0,23),
       labels=c("12.00am","1:00am","2:00am","3:00am","4:00am","5:00am","6:00am",
       "7:00am","8:00am","9:00am","10:00am","11:00am","12:00pm",
       "1:00pm","2:00pm","3:00pm","4:00pm","5:00pm","6:00pm","7:00pm",
       "8:00pm","9:00pm","10:00pm","11:00pm")) +
       coord_polar(theta="x") +
       ylab("HOUR OF DAY")+
       xlab("Source: Weather Underground (2014)")+
       ggtitle("Rain Event Frequency by Month and Hour\n(Bayan Weather Station)\n")+
       guides(colour = guide_legend(show = FALSE)) +
       theme(panel.background=element_blank(),
       axis.title.y=element_text(size=10,hjust=0.75,colour="grey20"),
       axis.title.x=element_text(size=7,colour="grey20"),
       panel.grid=element_blank(),
       axis.ticks=element_blank(),
       axis.text.y=element_text(size=5,colour="grey20"),
       axis.text.x=element_text(size=10,colour="grey20",face="bold"),
       plot.title = element_text(lineheight=1.2, face="bold",size = 14, colour = "grey20"),
       plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm"))
 
r1
 
# Save to png file
 
ggsave(r1,file="Rain_Event_Frequency_Plot.png", width=6, height=6, dpi=400,type="cairo-png")
Created by Pretty R at inside-R.org

Saturday, 19 July 2014

Examining seasonal & daily change in weather on Penang Island using circular plots


Key points of post

  • Circular plots are useful for summarizing weather data simultaneously exposing seasonal and daily changes;
  • The most uncomfortable (heat index) time of year for people living in Penang are the months of April and May although January and February are the hottest in absolute terms.

In this blog we demonstrate the use of circular plots which are handy for summarizing how any one variable is affected simultaneously by two others. The four circular plots below show how changes in average (based on the years 1980-2014) air temperature, wind-speed, relative humidity and heat index typically depend on both month and time of day. The circles are divided into 12 segments for each month, and into 24 concentric circles for each hour of the day, midnight being the center. The magnitudes of each variable are coded by color. The warmest months are February, March, and April, temperatures peaking between around 1pm and 4pm every day. The nights after about 11pm are clearly much cooler. Penang is only around 5 degrees north of the equator but there is nevertheless some variation around day-length and this is reflected in the relatively large light blue band on the outer most rings of the plot in October, November, and December. 

The patterns for wind speed are rather curious, presumably reflecting the Monsoons which typically switch between the North-East (December to February) and the South-West (April to August). Wind-speeds in Penang are apparently highest between January and March with a smaller peak seen in July and August.  They are also strongest between late morning and early afternoon in all seasons/months. Each day winds typically ‘get up’ around 11am blowing relatively strongly until around 5pm coinciding with the hottest part of the day.

Circular plots for selected weather variables showing how average (1980-2014) daily and seasonal patterns tend to covary


Typical seasonal daily changes in relative humidity are described by the third circular plot. Relative humidity is based on the difference between ‘dry bulb’ and ‘wet bulb’ air temperatures.  It describes how easily water can evaporate from any surface.  Since humans cool by sweating (the change in state of water as it evaporates from human skin causes a change in state from liquid to gas which requires a considerable amount of energy/heat, link to article on latent heat of evaporation) relative humidity is an important variable to consider.  High levels of relative humidity indicate that the air is already saturated with water vapour, and absorbing more will be difficult.  Similarly water will evaporate into air with low relative humidity more readily leading to improved function of the human air-conditioning system.  Air movement (or wind) also encourages evaporative loss from human skin improving cooling efficacy as anyone who has switched on a fan will know.

The combination of (dry bulb) air temperatures with data for relative humidity, and wind-speed facilitates the calculation of a ‘heat index’ which is more directly related to ‘human comfort’ than any one of its constituent variables. When the heat index is high the weather will be uncomfortable and, when low, more bearable.  In a previous blog post, we examined average air temperatures and rainfall in Penang suggesting that January and February were the worst times of year to visit Penang since they are usually the hottest and October the best since it is usually coolest. The heat index plotted in the circular plot suggests a slightly different story, however.  April and May have the highest heat indices and are probably potentially the most uncomfortable times of year to visit Penang while the entire period between July and December have similar heat indices.

To produce the plot above, we have again sourced the data from NOAA's National Climatic Data Centre with the hourly dataset coming from the Integrated Surface Database.

We also provide the R code as usual (below) if you would like to produce similar plots for your area of interest.

# Load package libraries
 
library(ggplot2)
library(scales)
library(plyr)
library(dplyr)
library(reshape2)
library(circular)
library(lubridate)
library(grid)
library(zoo)
library(gridExtra)
library(weathermetrics)
library(RColorBrewer)
 
# Setting work directory
 
setwd("d:\\ClimData")
 
# Reading and reformatting raw hourly data downloaded from NCDC
 
dat<-read.table("831677118578dat.txt",header=TRUE,fill=TRUE,na.strings=c("*","**","***","****","*****","******","0.00T*****"))
 
colnames(dat)<-tolower(colnames(dat))
 
Sys.setenv(TZ = "UTC")
dat$dates <- as.POSIXct(strptime(dat$yr..modahrmn,format="%Y%m%d%H%M"))  + 8 * 60 * 60
 
dat$year <- as.numeric(format(dat$dates,"%Y"))
 
dat$month <- as.numeric(format(dat$dates,"%m"))
 
dat$hour<-substring(as.character(dat$dates),12,13)
 
dat$min<-substr(dat$dates,15,16)
 
dat$time<-paste(dat$hour,dat$min,sep=":")
 
dat$tempc <- (dat$temp-32) * (5/9)
 
dat$dewpc <- (dat$dewp-32) * (5/9)
 
dat$tempc[dat$tempc<=10] <- NA
 
dat$tempc[dat$tempc>=40] <- NA
 
dat$dir[dat$dir == 990.0] <- NA 
 
dat$wspd <- (dat$spd)*0.44704
 
# Convert precipitation from inches to mms
 
dat$rain  <- dat$pcp24*25.4
 
# Calculate relative humidity & heat index using weathermetrics package
 
dat$rh <- dewpoint.to.humidity(t = dat$tempc, dp = dat$dewpc, temperature.metric = "celsius")
 
dat$hi <- heat.index(t = dat$tempc,rh = dat$rh,temperature.metric = "celsius",output.metric = "celsius",round = 2)
 
# Commands to reformat dates
 
dat$year <- as.numeric(as.POSIXlt(dat$dates)$year+1900)
dat$month <- as.numeric(as.POSIXlt(dat$dates)$mon+1)
dat$monthf <- factor(dat$month,levels=as.character(1:12),labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),ordered=TRUE)
dat$weekday <- as.POSIXlt(dat$dates)$wday
dat$weekdayf <- factor(dat$weekday,levels=rev(0:6),labels=rev(c("Mon","Tue","Wed","Thu","Fri","Sat","Sun")),ordered=TRUE)
dat$yearmonth <- as.yearmon(dat$dates)
dat$yearmonthf <- factor(dat$yearmonth)
dat$week <- as.numeric(format(as.Date(dat$dates),"%W"))
dat$hour <- as.numeric(format(strptime(dat$dates, format = "%Y-%m-%d %H:%M"),format = "%H"))
 
dat <- ddply(dat,.(yearmonthf),transform,monthweek=1+week-min(week))
 
# Extract data from 1980 onwards
 
dat1 <- subset(dat, year >= 1980 )
 
# Summarize data for weather variables by hour and month
 
dat2 <- ddply(dat1,.(monthf,hour),summarize, wspd = mean(wspd,na.rm=T))
 
dat3 <- ddply(dat1,.(monthf,hour),summarize, temp = mean(tempc,na.rm=T))
 
dat4 <- ddply(dat1,.(monthf,hour),summarize, hi = mean(hi,na.rm=T))
 
dat6 <- ddply(dat1,.(monthf,hour),summarize, rh = mean(rh,na.rm=T))
 
## Plot Temperature circular Plot
 
p1 = ggplot(dat3, aes(x=monthf, y=hour, fill=temp)) +
  geom_tile(colour="grey70") +
  scale_fill_gradientn(colours = c("#99CCFF","#81BEF7","#FFFFBD","#FFAE63","#FF6600","#DF0101"),name="Temperature\n(Degree C)\n")+
  scale_y_continuous(breaks = seq(0,23),
                     labels=c("12.00am","1:00am","2:00am","3:00am","4:00am","5:00am","6:00am","7:00am","8:00am","9:00am","10:00am","11:00am","12:00pm",
                              "1:00pm","2:00pm","3:00pm","4:00pm","5:00pm","6:00pm","7:00pm","8:00pm","9:00pm","10:00pm","11:00pm")) +
  coord_polar(theta="x") +
  ylab("HOUR OF DAY")+
  xlab("")+
  ggtitle("Temperature")+
  theme(panel.background=element_blank(),
  axis.title.y=element_text(size=10,hjust=0.75,colour="grey20"),
  axis.title.x=element_text(size=7,colour="grey20"),
  panel.grid=element_blank(),
  axis.ticks=element_blank(),
  axis.text.y=element_text(size=5,colour="grey20"),
  axis.text.x=element_text(size=10,colour="grey20",face="bold"),
  plot.title = element_text(lineheight=1.2, face="bold",size = 14, colour = "grey20"),
  plot.margin = unit(c(-0.25,0.1,-1,0.25), "in"),
  legend.key.width=unit(c(0.2,0.2),"in"))
 
p1
 
 
## Plot Wind Speed circular plot
 
p2 = ggplot(dat2, aes(x=monthf, y=hour, fill=wspd)) +
  geom_tile(colour="grey70") +
  scale_fill_gradientn(colours = rev(topo.colors(7)),name="Wind Speed\n(meter/sec)\n")+
  scale_y_continuous(breaks = seq(0,23),
                     labels=c("12.00am","1:00am","2:00am","3:00am","4:00am","5:00am","6:00am","7:00am","8:00am","9:00am","10:00am","11:00am","12:00pm",
                              "1:00pm","2:00pm","3:00pm","4:00pm","5:00pm","6:00pm","7:00pm","8:00pm","9:00pm","10:00pm","11:00pm")) +
  coord_polar(theta="x") +
  ylab("HOUR OF DAY")+
  xlab("")+
  ggtitle("Wind Speed")+
  theme(panel.background=element_blank(),
  axis.title.y=element_text(size=10,hjust=0.75,colour="grey20"),
  axis.title.x=element_text(size=7,colour="grey20"),
  panel.grid=element_blank(),
  axis.ticks=element_blank(),
  axis.text.y=element_text(size=5,colour="grey20"),
  axis.text.x=element_text(size=10,colour="grey20",face="bold"),
  plot.title = element_text(lineheight=1.2, face="bold",size = 14, colour = "grey20"),
  plot.margin = unit(c(-0.25,0.25,-1,0.25), "in"),
  legend.key.width=unit(c(0.2,0.2),"in"))
 
p2
 
 
# Plot Heat Index circular plor
 
p3 = ggplot(dat4, aes(x=monthf, y=hour, fill=hi)) +
  geom_tile(colour="grey70") +
  scale_fill_gradientn(colours = c("#99CCFF","#81BEF7","#FFFFBD","#FFAE63","#FF6600","#DF0101"),name="Temperature\n(Degree C)\n")+
  scale_y_continuous(breaks = seq(0,23),
                     labels=c("12.00am","1:00am","2:00am","3:00am","4:00am","5:00am","6:00am","7:00am","8:00am","9:00am","10:00am","11:00am","12:00pm",
                              "1:00pm","2:00pm","3:00pm","4:00pm","5:00pm","6:00pm","7:00pm","8:00pm","9:00pm","10:00pm","11:00pm")) +
  coord_polar(theta="x") +
  ylab("HOUR OF DAY")+
  xlab("")+
  ggtitle("Heat Index")+
  theme(panel.background=element_blank(),
  axis.title.y=element_text(size=10,hjust=0.75,colour="grey20"),
  axis.title.x=element_text(size=7,colour="grey20"),
  panel.grid=element_blank(),
  axis.ticks=element_blank(),
  axis.text.y=element_text(size=5,colour="grey20"),
  axis.text.x=element_text(size=10,colour="grey20",face="bold"),
  plot.title = element_text(lineheight=1.2, face="bold",size = 14, colour = "grey20"),
  plot.margin = unit(c(-0.5,0.25,-0.25,0.25), "in"),
  legend.key.width=unit(c(0.2,0.2),"in"))
 
p3
 
 
# Plot Relative Humidity circular plor
 
col<-brewer.pal(11,"Spectral")
 
p4 = ggplot(dat6, aes(x=monthf, y=hour, fill=rh)) +
  geom_tile(colour="grey70") +
  scale_fill_gradientn(colours=col,name="Relative Humidity\n(%)\n")+
  scale_y_continuous(breaks = seq(0,23),
                     labels=c("12.00am","1:00am","2:00am","3:00am","4:00am","5:00am","6:00am","7:00am","8:00am","9:00am","10:00am","11:00am","12:00pm",
                              "1:00pm","2:00pm","3:00pm","4:00pm","5:00pm","6:00pm","7:00pm","8:00pm","9:00pm","10:00pm","11:00pm")) +
  coord_polar(theta="x") +
  ylab("HOUR OF DAY")+
  xlab("")+
  ggtitle("Relative Humidity")+
  theme(panel.background=element_blank(),
  axis.title.y=element_text(size=10,hjust=0.75,colour="grey20"),
  axis.title.x=element_text(size=7,colour="grey20"),
  panel.grid=element_blank(),
  axis.ticks=element_blank(),
  axis.text.y=element_text(size=5,colour="grey20"),
  axis.text.x=element_text(size=10,colour="grey20",face="bold"),
  plot.title = element_text(lineheight=1.2, face="bold",size = 14, colour = "grey20"),
  plot.margin = unit(c(-0.5,0.1,-0.25,0.25), "in"),
  legend.key.width=unit(c(0.2,0.2),"in"))
 
p4
 
## Plot and export to png file
 
png(file="Weather_Variables_Circular_Plot.png",width = 12, height = 10, units = "in",
    bg = "white", res = 400, family = "", restoreConsole = TRUE,
    type = "cairo")
 
grid.arrange(p1,p2,p4,p3,nrow=2,main=textGrob("\nWeather Variables by Month and Hour\n(Bayan Lepas Weather Station)",
gp=gpar(fontsize=18,col="grey20",fontface="bold")),sub=textGrob("Source: NOAA National Climatic Data Centre (NCDC)\n",
                                                                                                              gp=gpar(fontsize=9,col="grey20",fontface="italic")))
 
dev.off()
Created by Pretty R at inside-R.org