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

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