The commuter raindance?

A recently published scientific paper by Haberlie et al. spotted by Jason (http://www.sciencedaily.com/releases/2015/02/150218123625.htm) showed that urbanization and patterns of commuting had an impact on local prevalence of thunderstorms. Areas that have been heavily built-on had more storms of greater intensity than those left, either under agricultural production, or more naturally vegetated. The researchers also reported a connection between commuting activity and extreme weather with more storms occurring during weekdays. The study mentioned focused on the sub-tropical city of Atlanta, Georgia in the US and we wondered if the same findings would apply to tropical Malaysia. 

We started by looking at the 'event' data that are part of the weather data collected by Weather Underground that can be pulled into R using the 'weatherData' package for Penang and noticed that thunderstorm events have actually decreased since 2001, see below.

Apart from the downward long-term trend between 2001 and 2015, the strong bimodal pattern of seasonality is also clear with peaks in April/May and September/October being typical each year. 

We then looked at the data for the capital city of Malaysia, Kuala Lumpur, using data from the International airport; an area which has undergone tremendous land use change over the last 10-15 years or so. Interestingly the data for Kuala Lumpur showed a significant increase in thunderstorm activity which is opposite the that observed in Penang.

We do not know the reason for the differences between Penang and Kuala Lumpur. Penang has also undergone substantial urbanization and land-use change but it is of course an island and the surrounding ocean presumably moderates its climate to some degree. 

In addition to the long-term changes there are also seasonal and daily cycles in the storm event data. The seasonality at KL is summarised below and the peaks in April and October/November reflect the south-west and north-east monsoons.

The daily cycles we observed in the thunderstorm event data are more interesting and perhaps puzzling. The circle plot below summarises the frequencies of storm events by time of day between 2001 and 2015 at KL International Airport. There are many more storms in the afternoon and early evening between 14:00 and 20:00 with a big peak at 17:00, ie. Rush hour. However if rush hour is the cause of the peak in storm events in the evening then why do you not see a similar pattern in the morning ?

We also investigate whether storms were more prevalent at KL during week days. The storm event data are 'counts' so we used a Generalised Linear Model from the Poisson family as follows. First we calculated a 'trend' vector and appended it to the data so January 2001=1, January 2002 = 13 and so on. We then fitted the following series of 'nested' models to the data:

z0 <- glm(no_event~1,data=storm0,family=quasipoisson) # The NULL model

z1 <- glm(no_event~trend,data=storm0,family=quasipoisson) # Is trend sig?

z2 <- glm(no_event~trend+monthf,data=storm0,family=quasipoisson)

z3 <- glm(no_event~trend+monthf+day,data=storm0,family=quasipoisson)

z4 <- glm(no_event~trend*monthf,data=storm0,family=quasipoisson)

z5 <- glm(no_event~trend*monthf+day,data=storm0,family=quasipoisson)

z0 is the NULL model against which the others are tested. Z1 tests whether trend is signficant and so on. Z3 and z5 test the effect of 'day of week' given that trend and monthf (month as a 'factor') are in the model. Z5 is only different from z3 in the sense that it 'allows' trend and monthf to interact or covary with one another. In English that means that the dependence on month can vary each year.

The results (below) show that there is no statistically significant effect of storm-event, ie. The number of storm events is the same on every day of the week. According to our analysis, therefore, commuting activity does not appear to be affecting thunderstorm initiation in Kuala Lumpur and the data do not support the findings of Haberlie et al.

anova(z0,z1,z2,z3,z4,z5,test='Chi')

Analysis of Deviance Table

Model 1: no_event ~ 1
Model 2: no_event ~ year
Model 3: no_event ~ trend + monthf
Model 4: no_event ~ trend + monthf + day
Model 5: no_event ~ trend * monthf
Model 6: no_event ~ trend * monthf + day
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 5098 6324.4
2 5097 6201.6 1 122.789 < 2.2e-16 ***
3 5086 5967.8 11 233.824 < 2.2e-16 ***
4 5080 5965.0 6 2.767 0.8558 NOT SIG
5 5075 5934.6 5 30.379 2.694e-05 ***
6 5069 5931.8 6 2.781 0.8543 NOT SIG


# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The R code used to produce the plots shown above are also provided below if you wish to replicate this at another location.

# Load required libraries
 
library(weatherData)
library(ggplot2)
library(data.table)
library(lubridate)
library(dplyr)
library(MASS)
library(scales)
library(extrafont)
library(reshape2)
library(plyr)
library(grid)
library(chron)
 
# Load fonts
 
loadfonts()
 
# Get weather data from Weather Underground using weatherData library
 
we <- getWeatherForDate("WMKK", "2001-01-01","2014-12-31", opt_detailed=T, opt_custom_columns=T, custom_columns=c(11))
 
# Convert to character
 
we$Events <- as.character(we$Events) 
 
# Assign values of '1' to Thunderstorm Event and '0' to Non-Thunderstorm Event
 
we$Events[we$Events == "" ] <- NA
we$Events[we$Events == "Rain"] <- 0
we$Events[we$Events == "Rain-Thunderstorm"] <- 1
we$Events[we$Events == "Thunderstorm"] <- 1
we$Events[is.na(we$Events)] <- 0
 
# Convert to numeric
 
we$Events<-as.numeric(we$Events)
 
# Calculate number of events
 
No_Events <- tapply(diff(c(0, we$Events)) == 1, as.Date((we$Time),tz="Asia/Kuala_Lumpur"), sum)
 
we.df <- data.frame(Date = as.Date(names(No_Events)), No_Event = unname(No_Events))
 
# Rename columns
 
colnames(we.df) <- c("date","no_event")
 
# Add date columns
 
we.df$month <- as.numeric(as.POSIXlt(we.df$date)$mon+1)
we.df$year <- strftime(we.df$date, "%Y")
we.df$week <- as.Date(cut(we.df$date,breaks = "week",start.on.monday = TRUE))
we.df$monthf <- factor(we.df$month,levels=as.character(1:12),
                labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),ordered=TRUE)
we.df$day <- strftime(we.df$date,'%A')
we.df$weekend = chron::is.weekend(we.df$date)
 
# Aggregate no events by month
 
month.agg <- aggregate(no_event ~ month + year, we.df, FUN = sum)
month.agg$date <- as.POSIXct(paste(month.agg$year, month.agg$month, "01", sep = "-"))
month.agg$monthf <- factor(month.agg$month,levels=as.character(01:12),
                    labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),ordered=TRUE)
 
# Plot number of events by month
 
ts <- ggplot(month.agg,aes(x=date,y=no_event,fill=no_event)) +
      geom_bar(stat="identity") +
      stat_smooth(method = "glm.nb",formula = y ~ x, data = month.agg, se = TRUE,colour="blue",size=0.5)+
      scale_fill_continuous(low="lightgoldenrod",high="red")+
      scale_x_datetime(labels = date_format("%Y"),breaks = "1 year")+
      ggtitle("Number of Monthly Thunderstorm Events from 2001 - 2014\n(Observations at KLIA, Selangor)\n")+
      theme_bw()+ ylab("No of events\n") + xlab("")+
      theme(axis.text.x  = element_text(size=11,family="Clear Sans",face="bold"), legend.position="none",
      axis.text.y  = element_text(size=11,family="Clear Sans",face="bold"), legend.position="none",
      panel.background=element_blank(),
      axis.title.y=element_text(size=12,colour="grey20",face="bold",family="Clear Sans"),
      panel.grid.major = element_line(colour = "grey",size=0.25,linetype='longdash'),
      plot.title = element_text(lineheight=1, face="bold",size = 15, colour = "grey20",family="Clear Sans"))
 
ts
 
ggsave(ts,file="D:/Thunderstorm_Events_KLIA_2000-2014.png",dpi=300,w=12,h=6,type="cairo-png")
 
month.mean <- ddply(month.agg,.(monthf),summarize,mean_events=mean(no_event,na.rm=T))
 
boxplot.title = 'Number of Thunderstorm Events by Month (2001 - 2014)'
boxplot.subtitle = 'Data source : Weather Underground (2015)'
 
# Box plot number of thunderstorm events by month
 
b <- ggplot(month.agg,aes(monthf,no_event,col=year)) +
     geom_boxplot(outlier.shape = NA,fill="ivory2",range=0.5,col="black")+
     geom_point(size=4,alpha=1/2,width=0.25,height=0.25,shape=18)+
     ggtitle("Number of Thunderstorm Events by Month at KLIA, Selangor (2001 - 2014) \n Source: Weather Underground (2015)\n")+
     ylab("No of events\n")+ xlab("") + stat_boxplot(geom ='errorbar',col="black")+
     scale_y_continuous(breaks=c(0,20,40,60),expand=c(0.1,0))+
     theme_bw()+
     geom_text(aes(label=year),hjust=-0.25, vjust=0,size=3,face="bold",family="Segoe UI Semibold")+
     theme(axis.text.x  = element_text(size=11,family="Clear Sans",face="bold"), legend.position="none",
     axis.text.y  = element_text(size=11,family="Clear Sans",face="bold"), legend.position="none",
     panel.background=element_blank(),
     axis.title.y=element_text(size=12,colour="grey20",face="bold",family="Clear Sans"),
     panel.grid.major = element_blank(),
     panel.grid.minor = element_blank(),
     plot.title = element_text(lineheight=1, face="bold",size = 15, colour = "grey20",family="Clear Sans"))
 
b
 
ggsave(b,file="D:/Thunderstorm_Events_byMonth_KLIA_2000-2014.png",dpi=300,w=12,h=6,type="cairo-png")
 
 
# Reformat weather event data from Weather Underground to get count of hourly thunderstorm events
 
we$hour <- as.numeric(format(strptime(we$Time, format = "%Y-%m-%d %H:%M"),format = "%H"))
we$min <- as.numeric(format(strptime(we$Time, format = "%Y-%m-%d %H:%M"),format = "%M"))
 
# Subset only events on the hour
 
we2<- subset(we, min == 0)
 
we3 <- ddply(we2,.(hour),summarize, event = sum(Events,na.rm=T))
 
# Plot number of thunderstorm events by hour of day
 
hr <- ggplot(we3,aes(x=hour,y=event))+
      geom_bar(breaks=seq(0,24),width = 1,colour="black",stat="identity",fill="blue",drop=TRUE)+ 
      coord_polar(start=0)+
      theme_minimal()+
      ylab("Count\n")+
      xlab("Source: Weather Underground (2015)")+
      scale_x_continuous(limits=c(0,24),breaks=seq(0,24),labels=seq(0,24))+
      ggtitle("Thunderstorm Event Frequency by Hour\nobserved at KLIA, Selangor\n")+
      guides(colour = guide_legend(show = FALSE)) +
      theme(panel.background=element_blank(),
      axis.title.y=element_text(size=12,hjust=0.65,colour="grey20",family="Clear Sans"),
      axis.title.x=element_text(size=10,colour="grey20",family="Clear Sans"),
      axis.ticks=element_blank(),
      axis.text.y=element_text(size=10,colour="grey20",family="Clear Sans"),
      axis.text.x=element_text(size=10,colour="grey20",face="bold",family="Clear Sans"),
      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"))
 
hr
 
ggsave(hr,file="D:/Thunderstorm_Events_byHourofDay_KLIA_2000-2014.png",dpi=300,w=8,h=8,type="cairo-png")

Created by Pretty R at inside-R.org

Rainfall event frequencies recorded by the Worldfish PWS in September 2014

Key points of post

·         Wednesday and Saturdays were far wetter.
·         Far more rain fell between 6 and 6.59 am and between 11 and 11.59pm

It has been just over a month and a half since we got the Weather Station here at WorldFish up and running and we thought it would fun to have a quick look at our first complete month of data, ie. data collected at 5 minute intervals for all of September (2014) typically one of the wettest months in Penang. Remember that a couple of months ago we looked at the frequency distribution of rainfall events by hour of day based on data from Bayan Lepas International airport.

After aggregating of the rainfall data, we looked at the frequency distribution of rainfall events by hour of day, and also by the day of week and we noticed some features we thought were worth sharing. One might be forgiven for assuming that there would not be too much variation in the amount of rain that falls either during a specific hour of the day, or a specific day of the week. However, the rainfall data captured here at the WorldFish office at Batu Maung show that certain times of  day and days of the week saw significantly more rainfall. In September 2014, we saw much more rain between 6 and 6.59am (65.9 mm) and between 11 and 11.59pm (48.5 mm). The rain falling at these times contributed over  30% of the rainfall for that entire month (354.5 mm) with 19% of the total monthly rainfall falling between 6 and 6.59am. Bizarrely it was also much more likely to rain on Wednesdays (102.6mm) and Saturdays (106.2mm's). To put this into perspective over 60% of September 2014’s  rain fell on either Saturday or Wednesday!

The plots also show that there was very little rainfall on Fridays during September 2014 (5.6 mm). Interestingly north-east Penang Island experienced a once-in-40-year rainfall event recently (137 mm's in a few hours) on FRIDAY night (3 October 2014). At the same time, however, very little rain was recorded at our PWS at Worldfish. Rainfall events can thus be extremely localized indicating the importance of establishing well-distributed and densely placed weather stations for studying weather and climate.

Granted, the rainfall data we looked at only covers a short period of time and the large variations could be pure coincidence. Conversely, these could also indicate consistent patterns or trends with scientific explanations possible (lunar periodicity, seasonal effects?). We intend to keep an eye on these features of the data and it will be interesting to see what emerges when time periods are examined.

As usual, the R code used the produce the plots in this post is provided below. The data is pulled from Weather Underground using the 'weatherData' package. 

## Load required libraries
 
library(weatherData)
library(ggplot2)
library(scales)
library(plyr)
library(reshape2)
library(lubridate)
library(zoo)
 
# Get data for PWS using weatherData package
 
pws <- getWeatherForDate("IPENANGB2", "2014-09-01","2014-09-30", station_type = "id",opt_detailed=T, opt_custom_columns=T, custom_columns=c(1,2,6,7,10,13))
 
# Rename columns
 
colnames(pws)<-c("time","time1","tempc","wdd","wspd","prcp","rain")
 
## Adding date columns
 
pws$time<-as.POSIXct(pws$time1,format="%Y-%m-%d %H:%M:%S",tz="Australia/Perth")
pws$year <- as.numeric(format(pws$time,"%Y"))
pws$date <-as.Date(pws$time,format="%Y-%m-%d",tz="Australia/Perth")
pws$year <- as.numeric(as.POSIXlt(pws$date)$year+1900)
pws$month <- as.numeric(as.POSIXlt(pws$date)$mon+1)
pws$monthf <- factor(pws$month,levels=as.character(1:12),labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),ordered=TRUE)
pws$weekday <- as.POSIXlt(pws$date)$wday
pws$weekdayf <- factor(pws$weekday,levels=rev(0:6),labels=rev(c("Mon","Tue","Wed","Thu","Fri","Sat","Sun")))
pws$yearmonth <- as.yearmon(pws$date)
pws$yearmonthf <- factor(pws$yearmonth)
pws$week <- as.numeric(format(as.Date(pws$date),"%W"))
pws$weekf<- factor(pws$week)
pws$jday<-yday(pws$date)
pws$hour <- as.numeric(format(strptime(pws$time, format = "%Y-%m-%d %H:%M"),format = "%H"))
pws$min <- as.numeric(format(strptime(pws$time, format = "%Y-%m-%d %H:%M"),format = "%M"))
 
# Remove duplicate values
 
# Function to compute consecutive differences 
dif <- function(x) c(diff(x), NA) 
 
# Apply the function to each individual 
 
pws<-ddply(pws, 'jday', transform, actrain = dif(rain)) 
 
pws$actrain[pws$actrain<0]<- 0
 
pws$actrain[is.na(pws$actrain)] <- 0
 
# Summarize rainfall data by hour and by day
 
rd <- ddply(pws,.(weekday,weekdayf),summarize, raintot = sum(actrain,na.rm=T))
 
rh <- ddply(pws,.(hour),summarize, raintot = sum(actrain,na.rm=T))
 
# Plot rainfall by hour
 
h <-  ggplot(rh,aes(hour,raintot,fill=raintot))+geom_bar(stat="identity")+
      geom_text(aes(label=round(raintot,3)), vjust = -0.4, size = 2.5)+
      ggtitle("Total rainfall by hour for September 2014 - WorldFish Weather Station\n")+
      scale_fill_continuous(low='grey90', high='steelblue') + theme_bw()+
      xlab("\nHour of Day")+
      ylab("Total Precipitation (mm)\n")+
      theme(axis.text.x  = element_text(size=10), legend.position="none",
      panel.background=element_blank(),
      axis.title.x=element_text(size=10,colour="grey20"),
      axis.title.y=element_text(size=10,colour="grey20"),
      panel.grid.major = element_line(colour = "grey",size=0.25,linetype='longdash'),
      axis.text.x=element_text(size=10,colour="grey20",face="bold"),
      plot.title = element_text(lineheight=1, face="bold",size = 13, colour = "grey20"))
 
h
 
# Plot rainfall by Week Day
 
d <- ggplot(rd,aes(weekday,raintot,fill=raintot))+geom_bar(stat="identity")+
     geom_text(aes(label=round(raintot,3)), vjust = -0.4, size = 2.5)+
     ggtitle("Total rainfall by day for September 2014 - WorldFish Weather Station\n")+
     scale_fill_continuous(low='grey90', high='steelblue') + theme_bw()+
     xlab("\nDay of Week")+
     ylab("Total Precipitation (mm)\n")+
     scale_x_continuous(breaks=c(0,1,2,3,4,5,6),labels=c("Mon","Tue","Wed","Thu","Fri","Sat","Sun"))+
     theme(axis.text.x  = element_text(size=10), legend.position="none",
     panel.background=element_blank(),
     axis.title.x=element_text(size=10,colour="grey20"),
     axis.title.y=element_text(size=10,colour="grey20"),
     panel.grid.major = element_line(colour = "grey",size=0.25,linetype='longdash'),
     axis.text.x=element_text(size=10,colour="grey20",face="bold"),
     plot.title = element_text(lineheight=1, face="bold",size = 13, colour = "grey20"))
 
d
 
# Save plots to png
 
ggsave(h,file="WF-PWS-RFbyHr.png",width=7,height=5,dpi=400,type="cairo")
ggsave(d,file="WF-PWS-RFbyDay.png",width=7,height=5,dpi=400,type="cairo")

Created by Pretty R at inside-R.org

Weather station at WorldFish HQ goes live

Worldfish HQ in Penang now has a weather monitoring station which is transmitting data live every five minutes to the internet where it is stored and served by Weather Underground (http://www.wunderground.com/personal-weather-station/dashboard?ID=IPENANGB2#history). 

Data and graphical output can be downloaded and examined by anyone in the world with an internet connection for free.   This is a really brilliant Citizen Science project, useful for both professionals and amateur enthusiasts, and we have been surprised at how simple it is to set up.  So far the majority of the weather stations on the Weather Underground network are in developed countries, but with the help of the CGIAR (Consultative Group on International Agricultural Research) of which we are a part, we think the network can expand increasingly into developing countries.  The next step for us is to establish weather stations at WorldFish offices worldwide (Zambia, Australia, Malawi, Bangladesh, Solomon Islands, Cambodia, Philippines, Myanmar and Timor-Leste).  From our headquarters in these countries we hope to expand the network out into the farming and fishing communities that we serve.

One substantial barrier for understanding the impact of climate variability and change on food production in developing countries is the lack of fine-scale spatial data, especially out in the sticks (boondocks). We think that Citizen Science projects like this can really help food producers, from farmers to fishers, smooth out risk from weather variability (perhaps using finance products like Index-Based Insurance).

If you are interested in setting up a weather station and linking it to Weather Underground please don’t hesitate to get in touch with us if you need any advice or further information.

The weather station we are using was bought in Penang (http://www.lacrossetechnology.com/2810/index.php) and installed very quickly by local contractors.  It consists of a temperature and humidity sensor, a self-emptying rain gauge and a solar powered anemometer (measures wind speed & direction).  All are linked wirelessly to a central control panel installed conveniently in our office.  Registration with Weather Underground is quick and straightforward allowing users to upload their weather data to the world-wide-web for free.

If you want to examine the data we are collecting at Worldfish HQ in more detail (see figure which shows weather variability at our office on August 26^th 2014)  you can use R (http://cran.r-project.org/) as outlined below.

Anemometer, wind vane and the tipping bucket self-emptying rain gauge

Temperature and Humidity Sensor (within radiation shield) outside the WorldFish NRM office

Control panel inside the WorldFish NRM office that collects
data from the sensors and transmits to a CPU via a USB dongle

# Code to plot weather variables from WorldFish Penang PWS linked via Weather Underground (Station: IPENANGB2)
 
# Load required libraries
 
library(weatherData)
library(ggplot2)
library(mgcv)
library(scales)
library(plyr)
library(reshape2)
library(gridExtra)
library(lubridate)
library(weathermetrics)
 
station.id="IPENANGB2"
date="2014-08-26"
 
# Get detailed weather station data
 
wd<-getDetailedWeather(station.id,date,station_type = "id", opt_all_columns=T)
 
str(wd)
 
# Rename columns
 
colnames(wd)<-c("time","time1","tempc","dewpc","slp","wd","wdd","wspd","guspd","hum","prcp","conditions",
                 "cc","rain","software","date_utc","unk")
 
wd$time<-as.POSIXct(wd$time1,format="%Y-%m-%d %H:%M:%S")
 
wd$rain[wd$rain < 0 ] <- NA
 
# Plot temperatures
 
dt_p <- ggplot(wd, aes(time, tempc)) + 
        xlab("Time") + ylab(as.expression(expression( paste("Temperature (", degree,"C)","\n")))) + 
        geom_line(colour="red",size=0.5) +
        geom_point(colour="red",size=1)+
        theme_bw() +
        ggtitle(paste('Plot of weather variables for',station.id, 'on',strftime(date,"%d-%b-%Y"),'\n\n',"Temperature\n")) +
        scale_x_datetime(labels = date_format("%I:%M:%p"),breaks = date_breaks("4 hours"))+
        theme(plot.title = element_text(lineheight=1.2, face="bold",size = 16, colour = "grey20"),
        panel.border = element_rect(colour = "black",fill=F,size=1),
        panel.grid.major = element_line(colour = "grey",size=0.25,linetype='longdash'),
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "ivory",colour = "black"))
 
dt_p
 
# Plot windspeed
 
wd$wspd[wd$wspd=="Calm"] <- 0
 
wd$wspd<-as.numeric(wd$wspd)
 
winds <- subset(melt(wd[,c("time","wspd")], id = "time"), value > 0)
 
dws_p <- ggplot(winds, aes(time, value))+
         geom_bar(stat="identity",col="limegreen",width=.5)+
         geom_point(data=wd,aes(time,guspd),col="seagreen4",size=1.5,shape=1)+        
         scale_x_datetime(labels = date_format("%I:%M:%p"),breaks = date_breaks("4 hours"))+
         ylab("km/hour\n")+
         xlab("Time")+
         scale_y_continuous(limits=c(0,20))+
         theme_bw()+
         theme(plot.title = element_text(lineheight=1.2, face="bold",size = 16, colour = "grey20"),
         panel.border = element_rect(colour = "black",fill=F,size=1),
         panel.grid.major = element_line(colour = "grey",size=0.25,linetype='longdash'),
         panel.grid.minor = element_blank(),
         panel.background = element_rect(fill = "ivory",colour = "black"))+        
         ggtitle("Wind and gust speed\n")
 
dws_p
 
# Plot wind vectors
 
wd$u <- (-1 * wd$wspd) * sin((wd$wdd) * pi / 180.0)
wd$v <- (-1 * wd$wspd) * cos((wd$wdd) * pi / 180.0)
dw = subset(wd, u != 0 & v != 0)
 
v_breaks = pretty_breaks(n = 5)(min(dw$v):max(dw$v))
v_labels = abs(v_breaks)
 
dwd_p <-  ggplot(data = dw, aes(x = time, y = 0)) +
          theme_bw() +
          theme(plot.margin = unit(c(0, 1, 0.5, 0.5), 'lines')) +
          geom_segment(aes(xend = time + u*360, yend = v), arrow = arrow(length = unit(0.15, "cm")), size = 0.5) + 
          geom_point(data = subset(dw, !is.na(u)), alpha = 0.5,size=1) +
          scale_x_datetime(name="Time",labels = date_format("%I:%M:%p"),breaks = date_breaks("4 hours")) +
          ggtitle("Wind vectors\n")+
          scale_y_continuous(name = "Wind vectors (km/h)\n", labels = v_labels, breaks = v_breaks)+
          theme(plot.title = element_text(lineheight=1.2, face="bold",size = 16, colour = "grey20"),
          panel.border = element_rect(colour = "black",fill=F,size=1),
          panel.grid.major = element_line(colour = "grey",size=0.25,linetype='longdash'),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "ivory",colour = "black"))
 
dwd_p
 
# Plot Sea Level Pressure
 
wd$slp <- as.numeric(as.character(wd$slp))
 
dslp_p <- ggplot(wd, aes(time, slp)) + 
          geom_point(col="purple4",size=1) + 
          geom_step(col="purple",size=0.5)+
          scale_x_datetime(labels = date_format("%I:%M:%p"),breaks = date_breaks("4 hours"))+
          ggtitle("Sea Level Pressure\n")+
          ylab("hPa\n\n")+
          xlab("Time")+
          theme_bw()+ 
          theme(plot.title = element_text(lineheight=1.2, face="bold",size = 16, colour = "grey20"),
          panel.border = element_rect(colour = "black",fill=F,size=1),
          panel.grid.major = element_line(colour = "grey",size=0.25,linetype='longdash'),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "ivory",colour = "black"))
 
dslp_p
 
 
# Plot Relative Humidity
 
dh_p <-  ggplot(wd, aes(time, hum)) + 
         geom_point(colour="royalblue4",size=1) +
         geom_step(col="royalblue",size=0.5)+
         xlab("Time") + ylab("Relative Humidity (%)\n") + 
         ggtitle("Relative humidity\n")+
         theme_bw() +
         scale_x_datetime(labels = date_format("%I:%M:%p"),breaks = date_breaks("4 hours"))+
         theme(plot.title = element_text(lineheight=1.2, face="bold",size = 16, colour = "grey20"),
         panel.border = element_rect(colour = "black",fill=F,size=1),
         panel.grid.major = element_line(colour = "grey",size=0.25,linetype='longdash'),
         panel.grid.minor = element_blank(),
         panel.background = element_rect(fill = "ivory",colour = "black"))
 
dh_p
 
# Plot Precipitation
 
dr_p <- ggplot(wd, aes(time, prcp)) + 
        xlab("Time") + ylab("Rainfall (mm)") + 
        geom_bar(stat="identity",colour="darkcyan",size=0.5,fill="darkcyan") +
        #geom_point(colour="darkcyan",size=2)+
        theme_bw() +
        ggtitle("Precipitation\n")+
        scale_x_datetime(labels = date_format("%I:%M:%p"),breaks = date_breaks("4 hours"))+
        theme(plot.title = element_text(lineheight=1.2, face="bold",size = 16, colour = "grey20"),
        panel.border = element_rect(colour = "black",fill=F,size=1),
        panel.grid.major = element_line(colour = "grey",size=0.25,linetype='longdash'),
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "ivory",colour = "black"))
 
dr_p
 
# Plot Cumulative Precipitation
 
dcr_p <- ggplot(wd, aes(time, rain)) + 
         xlab("Time") + ylab("Rainfall (mm)") + 
         geom_line(colour="darkcyan",size=0.5) +
         geom_point(colour="darkcyan",size=1)+
         theme_bw() +
         ggtitle("Cumulative Precipitation\n")+
         scale_x_datetime(labels = date_format("%I:%M:%p"),breaks = date_breaks("4 hours"))+
         theme(plot.title = element_text(lineheight=1.2, face="bold",size = 16, colour = "grey20"),
         panel.border = element_rect(colour = "black",fill=F,size=1),
         panel.grid.major = element_line(colour = "grey",size=0.25,linetype='longdash'),
         panel.grid.minor = element_blank(),
         panel.background = element_rect(fill = "ivory",colour = "black"))
 
dcr_p
 
 
# Calculate Heat Index using weathermetrics
 
wd$hi <- heat.index(t = wd$tempc,rh = wd$hum,temperature.metric = "celsius",output.metric = "celsius",round = 2)
 
# Plot Heat Index
 
dhi_p <- ggplot(wd, aes(time, hi)) + 
         xlab("Time") + ylab(as.expression(expression(paste("Temperature (", degree,"C)",)))) + 
         geom_line(colour="red",size=0.5) +
         geom_point(colour="red",size=1)+
         theme_bw() +
         ggtitle("Heat Index\n") +
         scale_x_datetime(labels = date_format("%I:%M:%p"),breaks = date_breaks("4 hours"))+
         scale_y_continuous(breaks=c(25,30,35,40,45))+
         theme(plot.title = element_text(lineheight=1.2, face="bold",size = 16, colour = "grey20"),
         panel.border = element_rect(colour = "black",fill=F,size=1),
         panel.grid.major = element_line(colour = "grey",size=0.25,linetype='longdash'),
         panel.grid.minor = element_blank(),
         panel.background = element_rect(fill = "ivory",colour = "black"))
 
dhi_p
 
# Export weather variables plot to png file
 
plot.name1 <- paste(station.id,"_WU_",date,".png",sep="")
 
png(plot.name1,width = 12, height = 25, units = "in",
     bg = "white", res = 500, family = "", restoreConsole = TRUE,
     type = "cairo")
 
grid.draw(rbind(ggplotGrob(dt_p),ggplotGrob(dh_p),ggplotGrob(dr_p),ggplotGrob(dcr_p),ggplotGrob(dslp_p),ggplotGrob(dws_p),ggplotGrob(dwd_p),ggplotGrob(dhi_p),size="first"))
 
dev.off()

Created by Pretty R at inside-R.org