Will I choke? – when is the best time for outdoor activities on Penang Island

Key points of post

• To assess whether air pollution index (API) levels in Penang during 2014 are unhealthy
• To assess how API varies seasonally and diurnally

Last year our blog post on air quality on Penang Island was particularly popular. In that post we used a comparatively crude ‘visibility index’ collected at Penang International Airport since 2001.  The data suggested that pollution (visibility) had gotten worse since 2001 and that it tended to be worse in June and July.

We discovered recently, however, that the Malaysian Government’s Department of the Environment (DOE) collects much more detailed data on air pollution at 52 stations across the country, see http://bit.ly/1wUMxYG.  Conveniently, one of the stations is located at Universiti Sains Malaysia (USM) and we've managed to acquire this data scraped off the DOE site by Keith Rozario (nice work Keith!) using python and MySQL. He provides a nice summary of the nationwide API data from August 2013 to February 2015 and the link to all the scripts he used (on github) to scrape the data on his site - please see http://bit.ly/1EnwTLs. He was also generous enough to make a zipped copy of the dataset available for the public to use under the Creative Commons Attribution 4.0 license.

Inspired by an analysis (http://bit.ly/1B6o2eM) by ‘Typhoon Tony’ from Taiwan, we sought to investigate the magnitude of air pollution on Penang Island, and how it varies seasonally and diurnally. 

Air pollution is typically quantified using the air pollution index (API) which is calculated from the following five major air pollutants: sulphur dioxide, nitrogen dioxide, carbon monoxide, particulate matter, and ozone, for details see http://bit.ly/1wUMxYG. APIs between 101 and 200 are considered to be ‘unhealthy’ especially for those with existing heart and/or lung problems.

So how does Penang stack up? We’re always hearing how terrible pollution is on the island, but is it actually true ? In the ‘calendar’ plot below, we merged the API data with hourly observations on wind-speed and direction from Weather Underground.  The seasonality in API is clear. We had high API values in February and March 2014, but the highest were in June and July which corresponds to the south-west monsoon transporting smoke from forest fires in Sumatera.  

We summarize diurnal variability in API (and to a certain extent seasonality too) in the ‘rose’ plot below.  It suggests that only rarely can air pollution on the island be considered ‘unhealthy’ since it is typically below 100.  If you are, however, concerned your best bet for strenuous outdoor activities is during the night or in the morning up until 1pm when levels of API start to increase.  The least healthy time of day to be out is rush hour between 4pm and 6pm. You may also want to avoid February, March, June and July!

Do note, however, that the API monitoring station at USM is relatively isolated from the island’s urban pell-mell and may not be a good overall reflection/summary of air pollution on the island. 

The R code used to produce the plots shown above are provided below. We are currently working on trying to scrape the API data off the Malaysian DOE website using R (which would make things a lot more convenient!) and once we've managed to do that, we will update the code accordingly.

# Load required libraries
 
library(ggplot2)
library(reshape)
library(lubridate)
library(zoo)
library(openair)
library(scales)
library(grid)
library(extrafont)
library(weatherData)
 
loadfonts()
 
# Set working directory
 
setwd("D:/Projects/AP/ReadingsByRegion/")
 
# Read data from csv file 
 
df <- read.table('USM.csv',header=F,sep=":")
 
# Rename columns
 
colnames(df) <- c("station_id","station_name","dates","hour","api","unk","unk1")
 
# Remove last 2 columns - not relevant information
 
df <- df[-c(6,7)]
 
# Add minute column
 
df$minute <- as.character("00")
 
df$time <- as.character(paste(df$hour,df$minute,sep = ":" ))
 
df$date <- as.POSIXct(paste(df$dates,df$time, format = "%Y-%m-%d %H:%M" ))
 
# Convert 0 to NA
 
df$api[df$api == 0] <- NA
 
# Date and time columns
 
df$year<-as.numeric(as.POSIXlt(df$date)$year+1900)
df$month<-as.numeric(as.POSIXlt(df$date)$mon+1)
df$monthf<-factor(df$month,levels=as.character(1:12),labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),ordered=TRUE)
df$weekday = as.numeric(format(as.POSIXlt(df$date),"%u"))
df$yearmonth<-as.yearmon(df$dates)
df$yearmonthf<-factor(df$yearmonth)
df$week <- as.numeric(format(as.POSIXlt(df$date),"%W"))
df$day = strftime(df$date,'%A')
df$jday <- strptime(df$date, "%Y-%m-%d")$yday+1
df$hour <- as.numeric(format(strptime(df$date, format = "%Y-%m-%d %H:%M"),format = "%H"))
 
 
##### Creating roseplot of pollutant magnitude over 24 hours #####
 
# Subset data for 2014
 
df2014 <- subset(df,year==2014)
 
# Set Color Table
 
colortable = c("#99FFFF", "#00FFFF", "#00FF00", "#CCFF33", "#FFFF00",
               "#FFCC00", "#FF6600", "#FF3333", "#FF33CC", "#660033")
 
# Cut data into ten parts
 
STN = "USM"
Time_of_Day = df2014$hour[df2014$station_name==STN]
mag = cut_number(round(df2014$api,100)[df2014$station_name==STN],n = 10)
rosedata = data.frame(dir=Time_of_Day,mag=mag)
 
# Plot rose chart
 
rc <- ggplot(rosedata,aes(x=Time_of_Day,fill=mag))+  geom_bar(binwidth = 1, drop = TRUE) +
      coord_polar() + xlab("") + ylab("") + 
      scale_x_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")) +
      ggtitle("\nDaily API readings in 2014") +  scale_fill_manual(values=colortable,name="API Reading")+
      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"),
      legend.title=element_text(size=13,colour="grey20",face="bold",family="Myriad Pro"),
      legend.text=element_text(size=11,colour="grey20",face="bold",family="Myriad Pro"),
      panel.grid=element_blank(),
      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, face="bold",size = 20, colour = "grey20",family="Graphik-Black"),
      plot.margin = unit(c(-0.5,0.25,-0.25,0.25), "in"),
      legend.key.width=unit(c(0.2,0.2),"in"))
 
rc
 
# Save plot to png
 
ggsave(file="D:/rc.png",dpi=300,w=10,h=8,unit="in",type="cairo-png")
 
 
### Reading weather data from Weather Underground using weatherData library ####
 
# Get weather data for 2014
 
we <- getWeatherForDate("WMKP", "2014-01-01","2015-12-31", opt_detailed=T, opt_custom_columns=T, custom_columns=c(1,6,7,8,9,13,14))
 
# Rename columns
 
colnames(we) <- c("date","date_myt","vis_km","wda","ws","gspd_kmh","wd","date_utc")
 
# Create date and time columns
 
we$date<-as.POSIXct(we$date)
we$hour <- as.numeric(format(strptime(we$date, format = "%Y-%m-%d %H:%M"),format = "%H"))
we$min <- as.numeric(format(strptime(we$date, format = "%Y-%m-%d %H:%M"),format = "%M"))
 
# Only use data on the hour
 
we_df<- subset(we, min == 0)
 
# Remove duplicate data
 
we_df$date[duplicated(we_df$date)]
 
# Merge API and weather data
 
we_api <- merge(we_df, df2014, by=c("date","date")) 
 
# Reformat yearmonth column for use in openair library
 
we_api$yearmonth <- as.factor(we_api$yearmonth)
 
# Reformat wind speed column
 
we_api$ws[we_api$ws == "Calm"] <- 0
we_api$ws <- as.numeric(we_api$ws)
 
# Plot 'calendar plot' and save as png
 
png(filename = "D:/Projects/AP/USM_API_2014_CalendarPlot.png",height=8,width=10,
    bg = "white",units='in', res = 400, family = "",  type = "cairo-png")
 
calendarPlot(we_api, pollutant = "api",year="2014",
             main = "Daily average API readings at USM\n
             with wind-speed scaled wind direction overlay (2014)\n",
             key.header = "API\n(Air Pollutant Index)",
             cols="heat", annotate = "ws")
 
dev.off()

Created by Pretty R at inside-R.org

Why are the British tweeting about climate change at night?

Points of post:

•  To demonstrate how the library twitteR can be combined with mapping capabilities of ggplot

• To reveal the locations of people tweeting about climate change at two points in time.

• To show that climate change tweeters are mostly from the developed world.

• To show that some (sad) tweeters are twittering about climate change in the middle of the night.

In a previous blog post, we used twitteR and wordcloud to summarize which other words were occurring in tweets (combined from around the world) together with the keywords, ‘climate-change’ and ‘Bangladesh’ at two arbitrarily selected time-points. 

The geo-locations of tweets are, however, also often available and are potentially very interesting and revealing. 

Fifteen-hundred (1500) is the maximum number of tweets that can be captured using the twitteR library with one call.  There probably are ways to get more data but we guess you probably have then to spend money. 

Here we used #climatechange because it coincided with the last days of COP20 in Lima, Peru which ran from 1^st to 12^th December 2014.  It is an extremely important global forum at which nations can meet and discuss their options for reducing carbon emissions, see http://unfccc.int/meetings/lima_dec_2014/meeting/8141.php)

The first ‘tweet map’ (below) we produced is based on approximately 1500 geo-located tweets that contained the hash-tag, #climatechange, and which were ‘tweeted’ at about 10am GMT on the 11^th December 2014. It shows that #climatechange tweets were coming from 4 main areas: North America, Europe, India and Australia. There didn’t appear either to be too many tweets coming out of Lima which surprised us. Maybe the delegates were too busy enjoying the South-American hospitality, and catching up with old mates to take much interest in Climate Change!

Geo-located tweets with #climatechange tweeted at around 10am GMT 11th December 2014

The second ‘tweet-map’ (below), also based on approximately 1500 geo-located #climatechange tweets, is for a snapshot that took place at 5 hours later at around 3am GMT on the final day of the conference (12^th December 2014).  The overall pattern between the maps remains the same but the relative frequency of #climatechange tweeters from Europe, as compared to North America, has increased. People in the United Kingdom were particularly keen, twittering like mad about climate change at 3am. Why? We don’t know.

Geo-located tweets with #climatechange tweeted at around 3am GMT 12th December 2014

Note that tweets are geo-located, either by exploiting the users’ location as defined in their profile, or by ascertaining the exact location directly  if allowed by the user. This can be effected, either from GPS-enabled software which many people have installed on their smart-phones, or by using an IP-address.  This means that not all tweets can be geo-located with any great precision. Some are only geo-located at the National and/or regional levels, as evident from the large circle in the middle of Australia. That’s to say these cautious tweeters only gave ‘Australia’ as their location.

As we have explained in an earlier blog post on worldclouds and twitteR, to pull the data from Twitter using its API, you will need to have a Twitter account and carry out a 'twitter authentication'.  The R code to perform a search on twitter for the selected 'term(s)' and mapping them out is detailed below.

# Load required libraries
 
library(RCurl)
library(maps)
library(stringr)
library(tm)
library(twitteR)
library(streamR)
library(grid)
library(ggplot2)
library(rgdal)
library(ggmap)
 
# Set working directory
 
setwd("D:/ClimData/")
 
#### Fonts on Windows ####
windowsFonts(ClearSans="TT Clear Sans")
 
# Load Credentials
 
load("D:/ClimData/Twitter/twitter authentification.Rdata")
registerTwitterOAuth(twitCred)
options(RCurlOptions = list(cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl")))
 
# Search term on twitter
 
searchTerm <- "#climatechange"
searchResults <- searchTwitter(searchTerm,n=1500,since='2014-12-11', until='2014-12-12')  
tweetFrame <- twListToDF(searchResults) 
 
userInfo <- lookupUsers(tweetFrame$screenName)  
userFrame <- twListToDF(userInfo)
 
locatedUsers <- !is.na(userFrame$location)
 
# Geocode locations using 'ggpmap' library
 
locations <- geocode(userFrame$location[locatedUsers])
 
locations_robin <- project(as.matrix(locations), "+proj=robin")
 
locations_robin_df <- as.data.frame(locations_robin)
 
# Import world boundaries
 
world <- readOGR(dsn="D:/Data/ne_10m_admin_0_countries", layer="ne_10m_admin_0_countries")
 
world_robin <- spTransform(world, CRS("+proj=robin"))
 
world_robin_df <- fortify(world_robin)
 
counts <- aggregate(locations_robin_df$V1,by=list(x=locations_robin_df$V1,y=locations_robin_df$V2),length)
names(counts)[3] <- "count"
 
# Theme options for Map
 
theme_opts <- list(theme(panel.grid.minor = element_blank(),
                         panel.grid.major = element_blank(),
                         panel.background = element_blank(),
                         panel.border = element_blank(),
                         plot.background = element_blank(),
                         axis.line = element_blank(),
                         axis.text.x = element_blank(),
                         axis.text.y = element_blank(),
                         axis.ticks = element_blank(),
                         axis.title.x = element_blank(),
                         axis.title.y = element_blank(),
                         legend.position = "bottom",
                         legend.key = element_blank(),
                         legend.title = element_text(colour="black", size=12, face="bold",family="Clear Sans"),
                         legend.text = element_text(colour="black", size=10, face="bold",family="Clear Sans"),
                         plot.title = element_text(size=15,face="bold",lineheight=0.5,family="Clear Sans")))
 
# Plot map and tweet counts 
 
tp <- ggplot(world_robin_df)+
      geom_polygon(aes(x = long, y = lat, group = group), fill = "grey20")+
      geom_path(aes(x = long, y = lat, group = group),colour = "grey40", lwd = 0.2)+
      geom_point(data= counts, aes(x=x,y=y,size=count),color="#32caf6", alpha=I(8/10))+
      scale_size_continuous(name="Number of tweets")+
      ggtitle("Twitter Map of #climatechange\n")+
      xlab("")+ ylab("")+
      coord_equal()+
      theme_bw() + 
      guides(size = guide_legend(title.position = "top",title.hjust =0.5))+
      theme_opts
 
tp
 
# Save to png
 
ggsave(tp,file="D:/Twitter_ClimateChange_Map.png",dpi=500,w=10,h=6,unit="in",type="cairo-png")

Created by Pretty R at inside-R.org

The recent deluge in Malaysia - using raincpc to map extreme rainfall events

Key point of post

• To describe an application of raincpc to map the rainfall that led to the recent (December 2014) floods in Malaysia

Raincpc (http://cran.r-project.org/web/packages/raincpc/vignettes/raincpc_demo.html) is a new library for R that exploits the Climate Prediction Center’s (CPC, www.cpc.ncep.noaa.gov) global (1979 to present, 50km resolution) datasets for precipitation. It renders CPC’s rainfall data more readily available for plotting and visualization, allowing any user to conveniently side-step problems relating to changing data formats, file-naming conventions etc. And all this free of charge!

We thought it would be fun to demonstrate the use of raincpc, focusing on Malaysia which experienced devastating floods over the Christmas and New Year Period, leading to the evacuation of 1000s of people. The damage has been estimated to have cost ~ 2 billion RM. Please see the following links - http://en.wikipedia.org/wiki/2014%E2%80%9315_Malaysia_floods and http://reliefweb.int/report/malaysia/asean-flash-update-northeast-monsoon-flood-24-december-2014.

In the plot below we used raincpc to show the amount of rain that fell over south-east Asia between 17^th and 24^th December 2014. It confirms that rainfall was indeed particularly heavy along the east coast of peninsular Malaysia; but also over northern Sumatera. Penang was certainly wet during December but the island had nothing like the amount of rainfall endured by communities on Malaysia’s east coast.

We do not know what caused the extreme rainfall that led to the flooding. Meteorologists think that it is related to the 'Madden-Julian Oscillation' (http://www.themalaysianinsider.com/malaysia/article/malaysia-could-see-more-severe-floods-like-in-kelantan-say-experts) and it's interaction with the north-east Monsoon. Very heavy rain is of course common in the tropics, but it doesn't neccesarily lead to flooding if drainage is adequate. Some experts think that rampant deforestation in Malaysia has led to more siltation of rivers, in effect blocking Malaysia's drains, and this exacerbates the impact of rainfall events (https://www.youtube.com/watch?v=r_eZUxgoxCw)

As usual the R-code for producing the map is outlined below.

## Load package libraries
 
library(raincpc)
library(SDMTools)
library(raster)
library(ggplot2)
library(rgdal)
library(grid)
library(maptools)
 
# Set working directory
 
setwd("D:/ClimData/")
 
## Get raw CPC rain data - data has a 2 day lag
 
cpc_get_rawdata(2014,12,17,2014,12,24) 
 
## Read raw CPC rain data into raster grids
 
rain1 <- cpc_read_rawdata(2014, 12, 17)
rain2 <- cpc_read_rawdata(2014, 12, 18)
rain3 <- cpc_read_rawdata(2014, 12, 19)
rain4 <- cpc_read_rawdata(2014, 12, 20)
rain5 <- cpc_read_rawdata(2014, 12, 21)
rain6 <- cpc_read_rawdata(2014, 12, 22)
rain7 <- cpc_read_rawdata(2014, 12, 23)
rain8 <- cpc_read_rawdata(2014, 12, 24)
 
# Combine multiple day rasters
 
rain_tot <- rain1 + rain2 + rain4 + rain5 + rain6 + rain7 + rain8
 
# Get summary of raster grid
 
print(rain_tot)
 
raster_ggplot <- function(rastx) {
 
require(SDMTools)
 
  stopifnot(class(rastx) == "RasterLayer")
 
  gfx_data <- getXYcoords(rastx)
  # lats need to be flipped
  gfx_data <- expand.grid(lons = gfx_data$x, lats = rev(gfx_data$y), 
                          stringsAsFactors = FALSE, KEEP.OUT.ATTRS = FALSE)
  gfx_data$rain <- rastx@data@values
 
  return (gfx_data)
}
 
rain_gg <- raster_ggplot(rain_tot)
 
# Read shapefile of country boundaries (shapefiles can be downloaded from http://thematicmapping.org/downloads/world_borders.php)
 
bounds <- readOGR(dsn="D:/Data/World_Borders", layer="TM_WORLD_BORDERS-0.3")
 
## Extents of ggplot map
 
xmin<-95
xmax<-120
ymin<--10
ymax<-15
 
interval <-(xmax-xmin)/5
 
lon_vals <- seq(xmin, xmax, 0.5)
lat_vals <- seq(ymin, ymax, 0.5)
 
 
# Set theme options
 
theme_opts <- list(theme(panel.grid.minor = element_blank(),
                         panel.grid.major = element_blank(),
                         panel.background = element_rect(fill="grey95"),
                         panel.border = element_rect(colour="black"),
                         axis.line = element_blank(),
                         axis.text.x = element_blank(),
                         axis.text.y = element_blank(),
                         axis.ticks = element_blank(),
                         axis.title.x = element_blank(),
                         legend.key.size=unit(0.35,"in"),
                         legend.key.width=unit(0.15,"in"),
                         legend.text=element_text(size=14,family="Myriad Pro Cond"),
                         legend.title=element_text(size=16,family="Myriad Pro Cond"),
                         plot.title = element_text(size=23,face="bold",vjust = -10,hjust=0.96,family="Graph Black"),
                         legend.position = c(0.17, 0), 
                         legend.justification = c(1, 0), 
                         legend.background = element_blank(),
                         axis.title.y = element_blank()))
 
# Plot rainfall map
 
rf <-  ggplot()+
       geom_raster(data=rain_gg,aes(x=lons,y=lats,fill=rain),alpha=0.8) +
       scale_fill_gradientn(colours=c("#FFFFFF","#ADD8E6","#95BBE9","#7E9EEC","#6781F0","#5064F3","#3948F6","#222BFA","#0B0EFD","#0A02FE","#1F06FC","#350AFA","#4A0EF8","#6013F6","#7517F3"),limits=c(0,1200),na.value=NA, name="Rainfall (mm)\n")+
       geom_polygon(data=bounds, aes(long,lat, group=group),colour="grey30",fill="transparent",size=0.35)+
       coord_equal(xlim=c(xmin,xmax),ylim=c(ymin,ymax)) + theme_bw()+
       ggtitle("Extreme rainfall event over Malaysia\n(17th to 24th of December 2014)\n")+
       xlab("") + ylab("")+ theme_opts +
       annotate("text", x = 115.00, y = -9.5, label = "Data source: Climate Prediction Center - NOAA (2014)",size=5,family="Myriad Pro Cond") 
 
plot(rf)
 
# Save rainfall map to png file
 
ggsave(rf,file="D:/ClimData/CPC_Extreme_Rainfall_Event_MYS_Dec2014.png",dpi=500,w=10,h=10,unit="in",type="cairo-png")

Created by Pretty R at inside-R.org