Hot events in Malaysia

Key points of post

• On Penang Island, the number of particularly hot days at or above 35⁰C has reduced since 1975, while near Kuala Lumpur, in the Klang Valley (Subang Weather Station) the trend has been strongly upwards.
• Long-term changes in average temperatures are not necessarily related to the number of hot events.

Commonly in reports and papers about climate change it is stated that, ‘climate variability and the number extreme events (very hot days, typhoons) has increased in the past and will continue to do so in the future’ but we wonder whether or not it is actually true everywhere.  

Jason recently chanced upon a post published on RPubs by James Elsner that produced analyses to examine this question. We decided to adopt this approach to examine whether or not the frequency of exceptionally hot days on Penang Island has actually gone up over the last few decades.  We know from our previous post that average air temperatures have increased and I (Doug) naively assumed that this must be true also for extremely hot days. Not so, however. In the picture below, using data from Bayan Lepas airport, we plot the number of days when air temperature was at or above 35⁰C at least once (top), the average duration of those hot periods (middle), and the number of such hot events with a line summarizing long-term trend (bottom). Clearly the number of such extreme events has gone down, contrary to popular wisdom about the effects of climate change. Interestingly the middle plot which summarizes successive days which had maximum temperatures of 35⁰C or more show that the average is only 1 day, the exception occurring in 2002 which had an average of 2. Surprisingly the average temperatures as shown in our post from February, were unexceptional.

But have the numbers of hot events decreased everywhere in Malaysia? At the Subang Weather Station, near Kuala Lumpur the answer is an emphatic, “No”. In 1998 which was a hot year (also see average air temperatures) there were 86 days when temperatures (at least once) were at or above 35⁰C when the average contiguous duration of hot days was ~ 2.5.  We don’t know the reason for these differences. The point is to be aware that these sort of differences exist. Subang located in the Klang Valley (refer to location map above) is inland, and is much less influenced by the moderating effects of the sea than Penang Island.  Differential ‘urban heat island’ effects may also be a factor. Or it could even be climate change.

I know where I’d rather live!

The R code used to produce the plots are shown below and the hourly weather data was again sourced from NOAA's National Climatic Data Centre (NCDC).

# Load required packages
 
library(ggplot2)
library(scales)
library(plyr)
library(dplyr)
library(reshape2)
library(lubridate)
library(grid)
library(zoo)
library(gridExtra)
 
# Setting work directory
 
setwd("d:\\ClimData")
 
# Reading and reformatting raw daily data downloaded from NCDC
 
dat<-read.table("CDO2812586929956.txt",header=F,skip=1)
 
colnames(dat)<-c("stn","wban","yearmoda","temp","tempc","dewp","dewpc","slp","slpc","stp","stpc","visib","visibc","wdsp","wdspc","mxspd","gust","maxtemp","mintemp","prcp","sndp","frshtt")
 
dat$yearmoda <- strptime(dat$yearmoda,format="%Y%m%d")
 
# Create date columns
 
dat$dates <- as.Date(dat$yearmoda)
 
dat$year <- as.numeric(format(dat$yearmoda,"%Y"))
dat$month <- as.numeric(format(dat$yearmoda,"%m"))
dat$day <- as.numeric(format(dat$yearmoda,"%d"))
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)
 
# Remove erronous values & convert units
 
dat$maxtemp[dat$maxtemp == 9999.9 ] <- NA
dat$mintemp[dat$mintemp == 9999.90 ] <- NA
 
dat$tempdc <- (dat$temp-32) * (5/9)
 
dat$maxtemp <- gsub("[*]","", dat$maxtemp) 
dat$mintemp <- gsub("[*]","", dat$mintemp) 
 
dat$maxtemp <- as.numeric(dat$maxtemp)
dat$mintemp <- as.numeric(dat$mintemp)
 
dat$maxtempc <- (dat$maxtemp-32) * (5/9)
dat$mintempc <- (dat$mintemp-32) * (5/9)
 
dat$maxtempc[is.na(dat$maxtempc)] <- as.numeric(dat$tempdc[is.na(dat$maxtempc)])
 
dat$mintempc[is.na(dat$mintempc)] <- as.numeric(dat$tempdc[is.na(dat$mintempc)])
 
str(dat)
 
dat<- dat[-c(3)]
 
# Subset data from 1975 onwards
 
dat <- subset(dat, year > 1975)
 
## Create a data frame table
 
dat = tbl_df(dat)
 
hd = dat %>% 
  group_by(year) %>%
  summarize(N35 = sum(maxtempc >= 35, na.rm = TRUE),
            N34 = sum(maxtempc >= 35, na.rm = TRUE),
            avgTmaxC = mean(maxtempc, na.rm = TRUE))
 
c1 <- ggplot(hd, aes(x = year, y = N35, fill = N35)) + 
      theme_bw() + 
      geom_bar(stat='identity') + 
      scale_fill_continuous(low='orange', high='red') +
      geom_text(aes(label = N35), vjust = 1.5, size = 3) +
      scale_x_continuous(breaks = seq(1975, 2015, 5)) +
      ylab("Days\n") +
      xlab("") +
      ggtitle("Number of days at or above 35°C\n")+
      theme(axis.text.x  = element_text(size=11), legend.position="none",
      panel.background=element_blank(),
      axis.title.x=element_text(size=7,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 = 12, colour = "grey20"))
 
c1
 
# Determine hot events
 
hotEvents = rle(dat$maxtempc >= 35)
eventLength = hotEvents$lengths[hotEvents$values]
eventNo = rep(1:length(eventLength), eventLength)
Events.df = subset(dat, maxtempc >= 35)
Events.df$eventNo = eventNo
 
# Determine the number of days between successive 35C days. Add this as another column.
 
t1 = Events.df$dates[-length(Events.df$dates)]
t2 = Events.df$dates[-1]
dd = difftime(t2, t1, units = "days")
Events.df$dbe = c(NA, dd)
Events.df = Events.df %>% select(maxtempc, mintempc, dates, year, month, eventNo, dbe)
 
# Length and intensity of events
 
LI.df = Events.df %>%
  group_by(eventNo) %>%
  summarize(eventLength = length(maxtempc),
            avgEventT = mean(maxtempc),
            maxEventT = max(maxtempc),
            whenMaxEvT = which.max(maxtempc),
            Year = year[1])
 
cor(LI.df$eventLength, LI.df$maxEventT)
 
ggplot(LI.df, aes(x = eventLength, y = whenMaxEvT)) +
  geom_point() +
  geom_smooth(method = "lm") +
  xlab("Event Length (days)") +
  ylab("Day of Event When Maximum Occurs") +
  scale_x_continuous(breaks = 1:10) +
  theme_bw()
 
# Average length of hot events
 
LI.df2 = LI.df %>%
  group_by(Year) %>%
  summarize(count = length(Year),
            avgEL = mean(eventLength))
 
AllYears = data.frame(Year = 1975:2014)
LI.df3 = merge(AllYears, LI.df2, by = "Year", all.x = TRUE)
LI.df3$count[is.na(LI.df3$count)] = 0
 
# Number of hot events

suppressMessages(library(MASS))

ne <- ggplot(LI.df3, aes(x = Year, y = count)) +
      geom_bar(stat = "identity") +
      ylab("Count\n") +
      scale_x_continuous(breaks = seq(1975, 2015, 5)) +
      ggtitle("Number of Hot Events in Penang since 1975\n")+
      stat_smooth(method = "glm.nb",formula = y ~ x, data = LI.df3, se = TRUE) + 
      theme_bw()+
      theme(legend.position = "none",
      axis.text.x  = element_text(size=11), legend.position="none",
      panel.background=element_blank(),
      axis.title.x=element_text(size=12,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 = 12, colour = "grey20"))
 
ne
 
var(LI.df3$count)/mean(LI.df3$count)
 
summary(glm.nb(count ~ Year, data = LI.df3))
 
## Export file to png
 
png(file="Hot_Events_Penang_1975-2014.png",width = 12, height = 10, units = "in",
    bg = "white", res = 500, family = "", restoreConsole = TRUE,type = "cairo")
 
grid.arrange(c1,el,ne,ncol=1,main=textGrob("\nHot Events in Penang (Bayan Weather Station) from 1975 - 2014",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

A Days’ Weather in Penang: the utility of Weather Underground

Key points of post

• Weather Underground combines a website, database and a global network of airport and Personal Weather Stations maintained largely by amateurs (http://www.wunderground.com/). 

• The network is growing fast, and the data can be accessed freely.

• The utility of Weather Underground is illustrated below using an example from one day’s weather in Penang, Malaysia.

Recently we’ve discovered Weather Underground (http://www.wunderground.com/) where you can access data from a global network of weather stations. Weather Underground is a ‘citizen science’ website specializing in the collection and dissemination of meteorological data.  It promotes a Personal Weather Station network whereby subscribers buy their own weather station, install it in their garden or wherever, and then plug into the global network via the internet. The data they collect are then automatically stored in the ‘cloud’ where they can be accessed, analyzed and visualized by anyone with a computer. Most Personal Weather Stations are in the United States, but other countries are catching up fast. There are three on Penang Island so far; one at Dalat International School, another in Tanjung Tokong, and the International Airport at Bayan Lepas is also connected. 

The data collected via the Personal Weather Stations can be imported automatically into our software of choice, R, (http://cran.r-project.org/) where they can be analyzed using the extensive statistical and data-visualisation tools available. 

As an example of the type of data available on Weather Underground, we downloaded and plotted data for a single day, collected at Bayan Lepas International Airport on 2^nd June 2014.  Data for six important meteorological parameters  were selected and have been plotted below at half-hourly intervals. 

Sometimes it can seem that weather in Penang is rather uniform throughout the day (ie. very hot), but this is actually not the case.  Starting at midnight, Air Temperatures were ~ 27°C until 8am when they rose rapidly to peak between 1.30 and 3.30pm (best time for a siesta).  After that temperature fell slowly, changing little between 6pm and 10pm remaining at 29°C.  In the afternoon of 2^nd June 2014 the wind speed increased from the south-west, after blowing very lightly from the north all night.  The rise in wind-speed during the day was related to a fall in barometric pressure, and this is of course the ‘South West Monsoon’ wind.  By dusk the wind speed fell away again.  ‘Heat-index’ is a measure of the comfort a human is likely to experience, combining measures of temperature and humidity. On 2^nd June this was clearly, directly related to the overall ambient temperature (top graph versus bottom) but this is not always the case.

 

Note: as usual the R-code for accessing Weather Underground and creating the graphs is given below. Readers can use this to plot data for their own location, dependent on the proximity of a Personal Weather Station.

We have used the 'weatherData' package developed by Ram Narasimhan to fetch the data from Weather Underground. The package can be installed directly from the CRAN site of from his site at http://ram-n.github.io/weatherData/

You can also get more information about the package from the github site that goes into further detail with examples and some of the other additional functions that come with it.  

 

## Load required libraries
 
library(weatherData)
library(ggplot2)
library(mgcv)
library(scales)
library(plyr)
library(reshape2)
library(gridExtra)
library(lubridate)
library(weathermetrics)
 
### Find required 4 letter station code using function below
 
# Example getStationCode("Penang") or getStationCode("Honiara")
 
getStationCode("Penang")
 
station.id="WMKP"
date="2014-06-02"
 
s<-getStationCode(station.id)
 
s1<-(strsplit(s,split= " "))[[1]]
 
station.name<-paste(s1[4],s1[5])
 
# Get detailed weather station data
 
wd<-getDetailedWeather(station.id,date, opt_all_columns=T)
 
str(wd)
 
# Rename columns
 
colnames(wd)<-c("time","mytime","tempc","dewpc","hum","slp","vsb","wd","wspd","guspd","prcp","events",
                "conditions","wdd","date_utc")
 
 
## 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=3)+
  theme_bw() +
  ggtitle(paste('Plot of weather variables for',station.name,"(",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_point(col="seagreen",size=3)+
  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,30))+
  stat_smooth(aes(group = 1), col="seagreen4",method = "loess",span=0.3,se=T,size=1.2)+
  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 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.25, "cm")), size = 0.75) + 
  geom_point(data = subset(dw, !is.na(u)), alpha = 0.5,size=3) +
  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=3) + 
  stat_smooth(span = 0.3,method="loess",col="purple",size=1.2)+
  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_step(colour="royalblue",size=0.5) +
  geom_point(colour="royalblue4",size=3) +
  xlab("Time") + ylab("Relative Humidity (%)\n") + 
  ggtitle("Relative humidity\n")+
  geom_smooth(aes(group = 1), col="royalblue",method = "loess",span=0.3,se=T,size=1.2)+
  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
 
## Calculate Heat Index using weathermetrics package
 
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=3)+
  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
 
plot.name1 <- paste(station.id,"_WU_",date,".png",sep="")
 
tiff(plot.name1,width = 13, height = 18, units = "in",
     compression = "lzw",bg = "white", res = 600, family = "", restoreConsole = TRUE,
     type = "cairo")
 
grid.draw(rbind(ggplotGrob(dt_p),ggplotGrob(dh_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