This blog is run by Jason Benedict and Doug Beare who are based in Penang, Malaysia to share insights and developments on open source software that can be used to analyze patterns and trends in weather and climate data; with particular emphasis on developing countries. The work featured here was undertaken as part of the CGIAR Research Program on Climate Change, Agriculture and Food Security (CCAFS) - http://ccafs.cgiar.org

However, the views presented here are of the authors only, and do not reflect the official position of CCAFS, CGIAR, its members, partners or donors.

Monday, 20 October 2014

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

Tuesday, 23 September 2014

Climate variability and change: hot air from Bangladesh?

Key points of post

  • Here we analyse air temperature from 3 weather stations along a north-south transect in Bangladesh (Rangpur (1982-2013), Dhaka (1998-2013) and Barisal (1998-2013).
  • Air temperatures risen at all 3 stations;  Dhaka had the most substantial increases in temperature rising 0.64°C overall between 1998 and 2013.
  • Seasonality in temperature varies between years but differences overall are slight and there seems to be no systematic changes in shape.

In this blog post, we ask whether air temperatures are rising in Bangladesh using data available to us from the internet. 

Bangladesh is a beautiful, semi-tropical country with India to the north and west and Myanmar to the east. The country is intersected by vast rivers, fed by the mountains of the Himalayas, which meander through the countryside. The normal seasonal cycles of flooding of these rivers regulates food production in many ways, e.g. irrigating rice or triggering fish production. Bangladesh is mostly low-lying and potentially vulnerable to sea-level rise and cyclones. The country is developing rapidly but much of its population remains poor, clinging to precarious livelihoods around rivers, lakes, and the sea.

Location of Rangpur, Dhaka and Barisal Weather Stations - Bangladesh

Bangladesh is a focal country for our CCAFS program which researches the impact of changing climate and weather on food production in developing countries; and particularly on the consequences for the poorest members of society.  Temperature is a ‘master’ variable, ultimately controlling all life in our biosphere including food production, species distributions and pathogen activity.

Here we examine long-term trends in temperature from the 3 following stations in Bangladesh: Rangpur; Dhaka; and Barisal (see map above). Unfortunately, compared to many other countries, data available for Bangladesh are rather sparse. Even for the capital city of Dhaka, we could only find continuous data from 1998, and there are also substantial gaps. The longest series available to us was from Rangpur which extended 32 years in total from 1982 to 2013. Air temperatures for the 3 stations are plotted below both as simple time-series and 3D, ‘heat-map’ plots, the latter being useful for showing how seasonality changes over longer term time.






The red lines in the 3 time-series plots represent linear models fitted to the data (for all 3 stations) as functions of year. The value of the gradient term (amount of rise/amount of tread) is indicated at the top-left. The gradients were positive in all 3 stations, meaning that temperatures indeed rose but the magnitude of increase is extremely slight compared to other locations we have examined, e.g. in Penang, Malaysia air temperature increases have been comparatively substantial. At Rangpur temperatures rose at 0.0117°C per year over 32 years which means an overall increase of 0.37°C overall (0.0117 x 32). Data were available for shorter periods at Dhaka and Barisal which recorded overall rises of 0.33°C and 0.64°C respectively. In Rangpur, extreme high temperatures (ie. > 35°C) occurred more frequently in the past, between 1982 and 1995, than would be the case for the present (see time-series above). Seasonality indeed varies from year to year (heat-map plots) but there is little evidence of the systematic changes over time that we have seen in other locations

The shape of the seasonal cycle seemed to vary only slightly between locations (see below). There is a clear Northern Hemisphere winter, summer cycle with lowest temperatures in December and January. In 2013, Dhaka was on average the warmest of all 3 locations, with peak temperatures occurring in June and July. Barisal was slightly warmer than Rangpur in the early part of the year after when Barisal fell relative to Rangpur. Peak temperatures in all 3 locations occurred (on average between June and July.


The results from these analyses are curious and possibly at odds with the view that the impacts of climate change on Bangladesh are especially worrying.  A rise of 0.64°C between 1998 and 2013 which we see at Dhaka is substantial, whereas further north at Rangpur we only saw a 0.37°C overall increase between 1982 and 2013.  In our view it is unlikely that these rises are even perceptible to the average human, and are also unlikely to have had much impact on food production.

Close examination of actual data can often challenge ‘accepted wisdom’, which is a particular purpose of these articles.  Of course, temperature is just one variable and Bangladesh’s climate future will also depend on the interaction between temperature and other variables such as rainfall and wind speed. Sea-level rise will also be a factor to consider. These will from the subject of future investigations.

The R code used to produce the plots in this blog post is provided below and the datasets used here are from NOAA National Climatic Data Centre's (NCDC) Integrated Surface Database.


# Load required libraries
 
library(ggplot2)
library(scales)
library(grid)
library(plyr)
library(lubridate)
library(zoo)
 
# Set working directory
 
setwd("D:/ClimData/BGD")
 
################################
## Plot Air Temps for Rangpur ##
################################
 
dat <- read.table("826316390028dat.txt",header=TRUE,fill=TRUE,
       na.strings=c("*","**","***","****","*****","******","0.00T*****","0.00T**","0.00T"))
 
colnames(dat)<-tolower(colnames(dat))
 
Sys.setenv(TZ = "UTC")
dat$dates <- as.POSIXct(strptime(dat$yr..modahrmn,format="%Y%m%d%H%M"))  + 6 * 60 * 60
 
dat$tempc <- (dat$temp-32) * (5/9)
 
dat$date <- as.Date(dat$dates,format="%Y%m%d")
 
dat$year <- as.numeric(format(dat$dates,"%Y"))
 
dat$month <- as.numeric(format(dat$dates,"%m"))
 
dat$hour  <- as.numeric(format(dat$dates,"%H"))
 
dat$week  <- as.numeric(format(dat$dates,"%W"))
 
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)
 
# Chuck out data that is < 0 and > 45
 
dat <- dat[dat$tempc > 0 & dat$tempc < 50, ]
 
datb<-subset(dat, year >= 1982 & year <= 2013)
 
plot.title = 'Daily Air Temperatures in Rangpur (1982-2013)'
plot.subtitle = 'Data source: Integrated Surface Database (ISD) - National Climatic Data Centre (NCDC)'
 
b <- ggplot(data=datb, aes(x=dates, y=tempc)) +
     geom_jitter(colour="grey40",size=1.5,alpha=0.4) +
     scale_y_continuous(breaks=c(10,15,20,25,30,35,40)) +
     theme_bw()+
     scale_x_datetime(breaks = "5 year", minor_breaks = "5 year", labels=date_format("%Y")) +
     xlab("") + ylab(as.expression(expression( paste("Temperature (", degree,"C)") ))) +
     ggtitle(bquote(atop(.(plot.title), atop(italic(.(plot.subtitle)), "")))) +
     theme(plot.title = element_text(face = "bold",size = 16,colour="black")) +
     geom_smooth(aes(group = 1), method = "lm",size=0.75,colour="red")
 
b
 
# Get linear equation value and gradient term
 
m <- lm(tempc~year, data=datb )
ms <- summary(m)
 
slope <- coef(m)[2]
 
lg <- list(slope = format(slope, digits=3))
 
eqg <- substitute(italic(Gradient)==slope,lg)
 
eqgstr <-as.character(as.expression(eqg))
 
# Add gradient term value to plot
 
b <- b +annotate(geom="text",as.POSIXct(-Inf, origin = '1970-01-01'), y = Inf, 
                      hjust = -0.1, vjust = 2, label = eqgstr,parse = TRUE,size=3.5)
 
b
 
col<-c("blue","green","yellow","orange","red")
 
wt <- ggplot(data=datb,aes(x=year,y=week,fill=tempc))+ geom_tile(colour=NA,size=0.65)+
      theme_bw()+
      scale_fill_gradientn(colours=col,name=as.expression(expression( paste("Temperature (", degree,"C)"))))+
      coord_equal(ratio=0.2)+
      ylab("WEEK OF YEAR\n")+
      xlab("\nYEAR\n\nSource: Integrated Surface Database (ISD) - National Climatic Data Centre (NCDC)")+
      scale_x_continuous(expand = c(0,0),breaks = seq(1982, 2015, 2)) +
      scale_y_discrete(expand = c(0,0),breaks = seq(0,52,4))+
      ggtitle("Average weekly air temperatures in Rangpur\n")+
      theme(panel.background=element_rect(fill="grey90"),
      panel.border=element_rect(colour="white"),
      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(),
      panel.grid.minor = element_blank(),
      legend.key.width=unit(c(0.1,0.1),"in"))
 
wt
 
# Exporting file to png
 
ggsave(b, file="Rangpur_AT_Plot.png", width=10, height=6,dpi=400,unit="in",type="cairo")
 
ggsave(wt, file="Rangpur_AT_HeatMap_Plot.png", width=15, height=6,dpi=400,unit="in",type="cairo")
 
 
##############################
## Plot Air Temps for Dhaka ##
##############################
 
dat <-read.table("820756390032dat.txt",header=TRUE,fill=TRUE,
      na.strings=c("*","**","***","****","*****","******","0.00T*****","0.00T**","0.00T"))
 
colnames(dat)<-tolower(colnames(dat))
 
Sys.setenv(TZ = "UTC")
dat$dates <- as.POSIXct(strptime(dat$yr..modahrmn,format="%Y%m%d%H%M"))  + 6 * 60 * 60
 
dat$tempc <- (dat$temp-32) * (5/9)
 
dat$year <- as.numeric(format(dat$dates,"%Y"))
 
dat$month <- as.numeric(format(dat$dates,"%m"))
 
dat$hour  <- as.numeric(format(dat$dates,"%H"))
 
dat$week  <- as.numeric(format(dat$dates,"%W"))
 
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)
 
# Chuck out data that is < 0 and > 45
 
dat <- dat[dat$tempc > 0 & dat$tempc < 50, ]
 
datb <- subset(dat, year >= 1998 & year <= 2013)
 
plot.title = 'Daily Air Temperatures in Dhaka (1998-2013)'
plot.subtitle = 'Data source: Integrated Surface Database (ISD) - National Climatic Data Centre (NCDC)'
 
b <- ggplot(data=datb, aes(x=dates, y=tempc)) +
     geom_jitter(colour="grey40",size=1.5,alpha=0.4) +
     scale_y_continuous(breaks=c(10,15,20,25,30,35,40)) +
     theme_bw()+
     scale_x_datetime(breaks = "1 year", minor_breaks = "1 year", labels=date_format("%Y"))  +
     xlab("") + ylab(as.expression(expression( paste("Temperature (", degree,"C)") ))) +
     ggtitle(bquote(atop(.(plot.title), atop(italic(.(plot.subtitle)), "")))) +
     theme(plot.title = element_text(face = "bold",size = 16,colour="black")) +
     geom_smooth(aes(group = 1), method = "lm",size=0.75,colour="red")
 
b
 
# Get linear equation value and gradient term
 
m <- lm(tempc~year, data=datb )
ms <- summary(m)
 
slope <- coef(m)[2]
 
lg <- list(slope = format(slope, digits=3))
 
eqg <- substitute(italic(Gradient)==slope,lg)
 
eqgstr <-as.character(as.expression(eqg))
 
# Add gradient term value to plot
 
b <- b +annotate(geom="text",as.POSIXct(-Inf, origin = '1970-01-01'), y = Inf, 
                 hjust = -0.1, vjust = 2, label = eqgstr,parse = TRUE,size=3.5)
 
b
 
col<-c("blue","green","yellow","orange","red")
 
wt <- ggplot(data=datb,aes(x=year,y=week,fill=tempc))+ geom_tile(colour=NA,size=0.65)+
      theme_bw()+
      scale_fill_gradientn(colours=col,name=as.expression(expression( paste("Temperature (", degree,"C)"))))+
      coord_equal(ratio=0.15)+
      ylab("WEEK OF YEAR\n")+
      xlab("\nYEAR\n\nSource: Integrated Surface Database (ISD) - National Climatic Data Centre (NCDC)")+
      scale_x_continuous(expand = c(0,0),breaks = seq(1998, 2015,1)) +
      scale_y_discrete(expand = c(0,0),breaks = seq(0,52,4))+
      ggtitle("Average weekly air temperatures in Dhaka\n")+
      theme(panel.background=element_rect(fill="grey90"),
      panel.border=element_rect(colour="white"),
      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(),
      panel.grid.minor = element_blank(),
      legend.key.width=unit(c(0.1,0.1),"in"))
 
wt
 
# Exporting file to png
 
ggsave(b, file="Dhaka_AT_Plot.png", width=10, height=6,dpi=400,unit="in",type="cairo")
 
ggsave(wt, file="Dhaka_AT_HeatMap_Plot.png", width=10, height=5,dpi=400,unit="in",type="cairo")
 
 
################################
## Plot Air Temps for Barisal ##
################################
 
dat <- read.table("8998396389910dat.txt",header=TRUE,fill=TRUE,
       na.strings=c("*","**","***","****","*****","******","0.00T*****","0.00T**","0.00T"))
 
colnames(dat)<-tolower(colnames(dat))
 
Sys.setenv(TZ = "UTC")
dat$dates <- as.POSIXct(strptime(dat$yr..modahrmn,format="%Y%m%d%H%M"))  + 6 * 60 * 60
 
dat$tempc <- (dat$temp-32) * (5/9)
 
dat$year <- as.numeric(format(dat$dates,"%Y"))
 
dat$month <- as.numeric(format(dat$dates,"%m"))
 
dat$hour  <- as.numeric(format(dat$dates,"%H"))
 
dat$week  <- as.numeric(format(dat$dates,"%W"))
 
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)
 
# Chuck out data that is < 0 and > 45
 
dat <- dat[dat$tempc > 0 & dat$tempc < 50, ]
 
datb<-subset(dat, year >= 1998 & year <= 2013)
 
plot.title = 'Daily Air Temperatures in Barisal (1998-2013)'
plot.subtitle = 'Data source: Integrated Surface Database (ISD) - National Climatic Data Centre (NCDC)'
 
b <- ggplot(data=datb, aes(x=dates, y=tempc)) +
     geom_jitter(colour="grey40",size=1.5,alpha=0.4) +
     scale_y_continuous(breaks=c(10,15,20,25,30,35,40)) +
     theme_bw()+
     scale_x_datetime(breaks = "1 year", minor_breaks = "1 year", labels=date_format("%Y"))  +
     xlab("") + ylab(as.expression(expression( paste("Temperature (", degree,"C)") ))) +
     ggtitle(bquote(atop(.(plot.title), atop(italic(.(plot.subtitle)), "")))) +
     theme(plot.title = element_text(face = "bold",size = 16,colour="black"))+
     geom_smooth(aes(group = 1), method = "lm",size=0.75,colour="red")
 
b
 
# Get linear equation value and gradient term
 
m <- lm(tempc~year, data=datb )
ms <- summary(m)
 
slope <- coef(m)[2]
 
lg <- list(slope = format(slope, digits=3))
 
eqg <- substitute(italic(Gradient)==slope,lg)
 
eqgstr <-as.character(paste(as.expression(eqg),as.expression(paste("/year"))))
 
# Add gradient term value to plot
 
b <- b +annotate(geom="text",as.POSIXct(-Inf, origin = '1970-01-01'), y = Inf, 
                 hjust = -0.1, vjust = 2, label = eqgstr,parse = TRUE,size=3.5)
 
b
 
col<-c("blue","green","yellow","orange","red")
 
wt <- ggplot(data=datb,aes(x=year,y=week,fill=tempc))+ geom_tile(colour=NA,size=0.65)+
      theme_bw()+
      scale_fill_gradientn(colours=col,name=as.expression(expression( paste("Temperature (", degree,"C)"))))+
      coord_equal(ratio=0.15)+
      ylab("WEEK OF YEAR\n")+
      xlab("\nYEAR\n\nSource: Integrated Surface Database (ISD) - National Climatic Data Centre (NCDC)")+
      scale_x_continuous(expand = c(0,0),breaks = seq(1998, 2015,1)) +
      scale_y_discrete(expand = c(0,0),breaks = seq(0,52,4))+
      ggtitle("Average weekly air temperatures in Barisal\n")+
      theme(panel.background=element_rect(fill="grey90"),
      panel.border=element_rect(colour="white"),
      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(),
      panel.grid.minor = element_blank(),
      legend.key.width=unit(c(0.1,0.1),"in"))
 
wt
 
# Exporting file to png
 
ggsave(b, file="Barisal_AT_Plot.png", width=10, height=6,dpi=400,unit="in",type="cairo")
 
ggsave(wt, file="Barisal_AT_HeatMap_Plot.png", width=10, height=6,dpi=400,unit="in",type="cairo")
 
 
############################################
## Plot Air Temps for All Stations (2013) ##
############################################
 
# Read raw data from files
 
datr <- read.table("826316390028dat.txt",header=TRUE,fill=TRUE,
        na.strings=c("*","**","***","****","*****","******","0.00T*****","0.00T**","0.00T"))
 
datd <- read.table("820756390032dat.txt",header=TRUE,fill=TRUE,
        na.strings=c("*","**","***","****","*****","******","0.00T*****","0.00T**","0.00T"))
 
datb <- read.table("8998396389910dat.txt",header=TRUE,fill=TRUE,
        na.strings=c("*","**","***","****","*****","******","0.00T*****","0.00T**","0.00T"))
 
# Bind data from 3 stations
 
bgt<-rbind(datr,datb,datd)
 
colnames(bgt)<-tolower(colnames(bgt))
 
bgt$usaf <- gsub("418590", "Rangpur", bgt$usaf)
bgt$usaf <- gsub("419220", "Dhaka", bgt$usaf)
bgt$usaf <- gsub("419500", "Barisal", bgt$usaf)
 
Sys.setenv(TZ = "UTC")
bgt$dates <- as.POSIXct(strptime(bgt$yr..modahrmn,format="%Y%m%d%H%M"))  + 6 * 60 * 60
 
bgt$tempc <- (bgt$temp-32) * (5/9)
 
bgt$year <- as.numeric(format(bgt$dates,"%Y"))
 
bgt$month <- as.numeric(format(bgt$dates,"%m"))
 
bgt$hour  <- as.numeric(format(bgt$dates,"%H"))
 
bgt$week  <- as.numeric(format(bgt$dates,"%W"))
 
bgt$monthf <- factor(bgt$month,levels=as.character(1:12),
              labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),ordered=TRUE)
 
bgt1 <-subset(bgt,year == 2013 )
 
plot.title = 'Mean Air Temperatures for 2013 (by Month)'
plot.subtitle = 'Source: Integrated Surface Database (ISD) - National Climatic Data Centre (NCDC)'
 
bt <-  ggplot(bgt1, aes(x = monthf, y = tempc,colour=factor(usaf),group=usaf))+
       geom_smooth(method="loess",se=F,size=1.2,span=1.5)+theme_bw()+
       xlab("\nMonth") + ylab(as.expression(expression( paste("Temperature (", degree,"C)"))))+
       ggtitle(bquote(atop(.(plot.title), atop(italic(.(plot.subtitle)), "")))) +
       theme(plot.title = element_text(face = "bold",size = 16,colour="black")) +
       guides(color = guide_legend(title = "Locality", title.position = "top")) +
       theme(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"),
       legend.position="right")
bt
 
# Exporting file to png
 
ggsave(bt, file="MonthlyAirTemps_BGD.png", width=10, height=6,dpi=400,unit="in",type="cairo")
Created by Pretty R at inside-R.org