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")
No comments:
Post a Comment