Boat Lights from the VIIRS: Shedding some light on night-time fishing activity in the Indo-Malayan waters

Light emitted at night due to human activities can easily be seen from space and detected by satellite. Many compelling pictures of the Earth at night have been taken revealing cities, roads, and conurbations. Light emitted at sea has, however, received less attention. Applications relating to fisheries and their management might prove particularly useful. NOAA/NASA and partners have recently developed a boat detection product derived from satellite imagery at high temporal and spatial resolution and the data can be downloaded for free from here as either kml/csv files and processed in R. 

For many years (since 1972 to be precise), global low light night imaging was collected by the US Air Force's Defence Meteorological Satellite Programmes Operational Linescan System (DMSP-OLS). While night-time lights collected from DMSP-OLS has major shortcomings, there were no reasonable alternative products for mapping night-lights of lit fishing boats. This all changed in 2011 with the launch of NASA/NOAA Suomi National Polar Partnership satellite equipped with the Visible Infrared Radiometer Suite (VIIRS) as the primary imaging sensor. You can find out many of the advantages the new VIIRS sensor has over the DMPS-OLS here. One of the most important improvements is that the 'Day/Night band' on the VIIRS that is used for night-time light detection has a detection footprint 45 times smaller than that on the DMSP, which ultimately means it can detect light sources a lot more smaller than what the DMPS-OLS could! The VIIRS also has a better dynamic range and lower detection limits which means it can identify and discern more sources of light at night with less saturation when it gets too bright with surrounding lights. These two improvements are key in enabling better detection of boat lights at night.

Commercial fishermen the world over use very bright light to attract squid and small pelagics to their nets, for example, and these are visible from space. Here we use the NOAA boat light data to animate nightly boat lights during February 2015 (cumulative for entire night) for a part of south-east Asia including Indonesia, Malaysia and the northern tip of Australia. We think that many of these lights reflect fishing activity, but more work needs to be done to detect and isolate shipping etc. A concerted fishers' knowledge/consultation process would be essential for enabling us to understand these data properly. We think it would it would be worth the effort, though, and could yield fascinating insights into fleet behavior and the species targeted.

As usual, we provide the R code below on how we produced the above animation of the boat lights detection product.

##### Code to read and produce animation of boat night light detection data from NGDC site #####
 
# Load required libraries
 
library(RCurl)
library(rvest)
library(dplyr)
library(lubridate)
library(zoo)
library(ggplot2)
library(rgdal)
library(grid)
library(animation)
library(gridExtra)
library(extrafont)
 
# Set working directory
 
setwd("D:/Projects/BoatLights")
 
# Load fonts
 
loadfonts()
 
# Set date range for csv files to be downloaded from NGDC site
 
date.range <- seq.Date(from=as.Date('2015-02-01'),as.Date('2015-02-28'),by='1 day')
 
# Get dates into relevant format
 
date <- as.character(date.range)
date1 <- gsub("-","",date)
 
 
# Get urls to scrape boat light csv files off  NGDC site
 
urls <- sprintf((paste("http://mapserver.ngdc.noaa.gov/viirs_data/indo_boats///",c(date1),"/VNB_npp_d",c(date1),
        "_ssec_v12.csv",sep="")))
 
# Read boat light data files in csv format and bind them
 
boats <- do.call(rbind,lapply(urls,read.csv))
 
# Extract useful columns
 
boats.df <- boats[-c(18,19)]
 
# Reformat date and add other date and time columns
 
boats.df$Date <- as.POSIXct(boats.df$Date_Mscan, format="%Y-%m-%d %H:%M:%S",tz="GMT")
boats.df$Dates <- as.Date(boats.df$Date)
boats.df$Year <- as.numeric(as.POSIXlt(boats.df$Date)$year+1900)
boats.df$Month <- as.numeric(as.POSIXlt(boats.df$Date)$mon+1)
boats.df$Monthf <- factor(boats.df$Month,levels=as.character(1:12),
                 labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),ordered=TRUE)
boats.df$Weekday = as.numeric(format(as.POSIXlt(boats.df$Date),"%u"))
boats.df$Yearmonth <- as.yearmon(boats.df$Date)
boats.df$Yearmonthf <- factor(boats.df$Yearmonth)
boats.df$Week <- as.numeric(format(as.POSIXlt(boats.df$Date),"%W"))
boats.df$Day <- strftime(boats.df$Date,'%A')
boats.df$Jday <- strptime(boats.df$Date, "%Y-%m-%d")$yday+1
boats.df$Hour <- as.numeric(format(strptime(boats.df$Date, format = "%Y-%m-%d %H:%M"),format = "%H"))
 
# Extract only detections that are classified as 'boat' from quality flag for nightboat detection algorithm
 
boats.df1 <- subset(boats.df,QF_Detect==1)
 
# Import shapefiles
 
world <- readOGR(dsn="D:/GeoSpatial/Data/World_Boundaries", layer="TM_WORLD_BORDERS-0.3") # World boundary shp's download from www.thematicmapping.org
 
eez <- readOGR(dsn="D:/Projects/BoatLights", layer="EEZs") # EEZ boundary shp's downloaded from www.marineregions.org
 
mpa <- readOGR(dsn="D:/Projects/BoatLights", layer="MPAs") # MPA boundary shp's downloaded from CT Atlas website - www.ctlas.reefbase.org
 
# Fortify layers for ggplot mapping
 
world.df <- fortify(world)
 
eez.df <- fortify(eez)
 
mpa.df <- fortify(mpa)
 
# Set theme options for map
 
theme_opts <- list(theme(panel.grid.minor = element_blank(),
                         panel.grid.major = element_blank(),
                         panel.background = element_rect(fill="black"),
                         panel.border = element_rect(colour="black"),
                         legend.text = element_text(size = 8,face="bold",family="Myriad Pro"),
                         legend.position = "right",
                         legend.key.width=unit(0.75,"line"),
                         legend.key.height=unit(0.75,"line"),
                         legend.direction="horizontal",
                         axis.line = element_blank(),
                         axis.text.x = element_blank(),
                         axis.text.y = element_blank(),
                         axis.ticks = element_blank(),
                         axis.title.x = element_blank(),
                         plot.title = element_text(size=20,face="bold",family="Graph Black",hjust=0,vjust=-0.5),
                         axis.title.y = element_blank()))
 
 
# Run loop to produce boat light maps for each night in February 2015
 
for (i in date)
 
{
 
 
b <- ggplot(subset(boats.df1,Dates ==i),aes(Lon_DNB,Lat_DNB))+
     geom_polygon(data=mpa.df, aes(long,lat, group=group,fill="MPA"),colour=NA,alpha=0.5)+
     geom_point(color = "gold", size = 0.5, alpha = 0.5) +
     geom_polygon(data=eez.df, aes(long,lat, group=group,colour="EEZ"),fill=NA,size=0.25,linetype="dashed")+
     geom_polygon(data=world.df, aes(long,lat, group=group),fill="black",colour="darkgray",size=0.35)+
     theme_bw()+ ylab("")+ xlab("")+
     ggtitle((paste("Boat Lights detected from VIIRS on ",i,"\n"))) +
     scale_fill_manual(name="",values=c(MPA="forestgreen"))+
     scale_colour_manual(name="",values=c(EEZ="dimgray"))+
     coord_equal(xlim=c(95,140),ylim=c(-15,9))+
     theme_opts
 
b
 
bg <- arrangeGrob(b, sub = textGrob("Source: VIIRS Boat Detection (VBD) Data (PROTOTYPE) - NGDC (2015)\n", hjust = -0.33,
      vjust=-0.3,gp = gpar(fontface = "bold", fontsize = 8, fontfamily="Tw Cen MT Condensed")))
 
bg
 
ggsave(bg,file=paste("D:/Projects/BoatLights/boatlights_",i,".png",sep=""),width=8,height=5,unit="in",type="cairo-png")
 
}
 
# Set path for Image Magick convert program 
 
path.to.convert <- paste0(shortPathName("C:/Program Files/ImageMagick-6.9.0-Q16)/"),"convert.exe")
ani.options(convert=path.to.convert)
 
 
# Convert png files to GIF animation
 
files = sprintf('boatlights_%s.png', c(date))
im.convert(files, output = 'BoatLights_Feb2015.gif')

Created by Pretty R at inside-R.org

Mapping Typhoon Haiyan/Yolanda with rNOMADS

The purpose of this article is to:

• Promote the utility of the recently published rNOMADS software package which links directly to the climate model archive maintained by NOAA.
• Demonstrate the use of rNOMADS for plotting global surface temperatures and wind-speeds.
• Show how we used rNOMADS to map the extent and track of Typoon Yolanda/Haiyan.

rNOMADS allows R users to connect directly with the National Oceanic and Atmospheric Administration's (NOAA) Operational Model Archive and Distribution System (NOMADS). It allows rapid and efficient downloading of a vast range of different types of meteorological data from both global and regional weather models. Data for the past are available in addition to data for all sorts of forecasts.

Global surface temperatures from the NOAA archives (details of actual model used) are shown in the map below for 18:00 GMT on 29th October 2014. The band of high sea surface temperatures associated with the tropics is clear. It is also interesting to note the influence of the ‘upwelling zones’ on the western coasts of the southern American and Africa continents which tend to force cold water northwards. On land the relatively cool areas associated with the Andes and the Ethiopian Plateau are also visible.

We follow up global temperatures with the map of global windspeeds at ground level below. The ‘Roaring Forties’, ‘Furious Fifties’ and ‘Screaming Sixties’ are very clearly discernible. So too the ‘Doldrums’, or areas of calm wind which hug the equator . You can also make out Cyclone Nilofar in the north Indian Ocean which was about to strike the west coast of India.

The gif animation below utilizes the same data source to summarize the path of Typhoon Haiyan/Yolanda which hit the Philippines in November 2013 causing great devastation and loss of life.  The animation depicts wind-speed (km/hr) at 3 time points: 15:00 GMT and 21:00 GMT on 7th November 2013 and then early morning (03:00 GMT) on 8th November. Each raster grid is roughly 50 km's by 50 km's (0.5 degrees by 0.5 degrees) in extent, so the whole storm was around 650 km's across.  At its height wind-speeds exceeded 150 km/hr and the ‘eye of the storm’ is also clearly visible as the blue dot in the very centre.  

As usual, the R code is provided below if you are interested in producing similar plots for an area or time of interest. You can get plenty of information and code examples for the rNOMADS package on CRAN (link provided earlier in the post) as well as at the developers' blog posts on the package at the following link - https://bovineaerospace.wordpress.com/tag/rnomads/

# Load libraries
 
library(rNOMADS) 
library(GEOmap)
library(aqfig)
library(rgdal)
 
# Define levels
 
levels <- c("surface", "2 m above ground", "10 m above ground", "300 mb")
 
# Define variables - temperature and wind
 
variables <- c("TMP", "UGRD", "VGRD")
 
# Define additional model inputs
 
abbrev <- "gfsanl"
model.date <- 20141029
model.run <- 06
pred <- 3
 
# Get the data
 
grib.info <- ArchiveGribGrab(abbrev, model.date,model.run, pred, file.type = "grib2")
 
# Read data into R
 
grib.data <- ReadGrib(grib.info$file.name, levels, variables)
 
resolution <- c(0.5, 0.5) #Resolution of the model
 
# Make an array for easier manipulation
 
atmos <- ModelGrid(grib.data, resolution)
 
# Set up display
 
# Plot temperature at 2 m above ground
 
li <- which(atmos$levels == "2 m above ground")
vi <- which(atmos$variables == "TMP")
colormap <- rev(rainbow(500, start = 0 , end = 5/6))
 
# Read world boundaries
 
world <- readOGR(dsn="D:/Data/ne_10m_admin_0_countries", layer="ne_10m_admin_0_countries")
 
png(file = "world_surface_temp.png", width = 9, height = 6, res=400, type="cairo-png",units="in", bg="white",restoreConsole = TRUE)
 
image(atmos$x , sort(atmos$y), atmos$z[li,vi,,], col = colormap,
      xlab = "Longitude", ylab = "Latitude",
      main = paste("World Temperature at Ground Level (deg C):", atmos$fcst.date ,"GMT"))
 
# Plot coastlines
 
plot(world, border = "black", add = TRUE, MAPcol = NA)
 
# Plot legend, convert to Celsius
 
vertical.image.legend(col=colormap, zlim = range(atmos$z[li,vi,,] - 273.15))
dev.off()
 
# Plot wind magnitude at 10 m above ground
 
li <- which(atmos$levels == "10 m above ground")
vi.zonal <- which(atmos$variables == "UGRD") # East-West wind
vi.merid <- which(atmos$variables == "VGRD") # North-South wind
 
wind.mag <- sqrt(atmos$z[li,vi.zonal,,]^(2) + atmos$z[li,vi.merid,,]^(2))
colormap <- rev(rainbow(500, start = 0 , end = 5/6))
 
png(file = "world_surface_wind.png", width = 9, height = 6, res=400, type="cairo-png",units="in", bg="white",restoreConsole = TRUE)
 
image(atmos$x, sort(atmos$y), wind.mag, col = colormap,
      xlab = "Longitude", ylab = "Latitude",xlim=c(100,150),ylim=c(-5,20),
      main = paste("World Winds at Ground Level (km/hr):", atmos$fcst.date, "GMT"))
 
# Plot coastlines
 
plot(world, border = "black", add = TRUE, MAPcol = NA)
 
# Plot legend, convert to km/hr
 
vertical.image.legend(col=colormap, zlim = range(wind.mag * 3.6))
 
dev.off()

Created by Pretty R at inside-R.org

Examining seasonal & daily change in weather on Penang Island using circular plots

Key points of post

• Circular plots are useful for summarizing weather data simultaneously exposing seasonal and daily changes;

• The most uncomfortable (heat index) time of year for people living in Penang are the months of April and May although January and February are the hottest in absolute terms.

In this blog we demonstrate the use of circular plots which are handy for summarizing how any one variable is affected simultaneously by two others. The four circular plots below show how changes in average (based on the years 1980-2014) air temperature, wind-speed, relative humidity and heat index typically depend on both month and time of day. The circles are divided into 12 segments for each month, and into 24 concentric circles for each hour of the day, midnight being the center. The magnitudes of each variable are coded by color. The warmest months are February, March, and April, temperatures peaking between around 1pm and 4pm every day. The nights after about 11pm are clearly much cooler. Penang is only around 5 degrees north of the equator but there is nevertheless some variation around day-length and this is reflected in the relatively large light blue band on the outer most rings of the plot in October, November, and December. 

The patterns for wind speed are rather curious, presumably reflecting the Monsoons which typically switch between the North-East (December to February) and the South-West (April to August). Wind-speeds in Penang are apparently highest between January and March with a smaller peak seen in July and August.  They are also strongest between late morning and early afternoon in all seasons/months. Each day winds typically ‘get up’ around 11am blowing relatively strongly until around 5pm coinciding with the hottest part of the day.

Circular plots for selected weather variables showing how average (1980-2014) daily and seasonal patterns tend to covary

Typical seasonal daily changes in relative humidity are described by the third circular plot. Relative humidity is based on the difference between ‘dry bulb’ and ‘wet bulb’ air temperatures.  It describes how easily water can evaporate from any surface.  Since humans cool by sweating (the change in state of water as it evaporates from human skin causes a change in state from liquid to gas which requires a considerable amount of energy/heat, link to article on latent heat of evaporation) relative humidity is an important variable to consider.  High levels of relative humidity indicate that the air is already saturated with water vapour, and absorbing more will be difficult.  Similarly water will evaporate into air with low relative humidity more readily leading to improved function of the human air-conditioning system.  Air movement (or wind) also encourages evaporative loss from human skin improving cooling efficacy as anyone who has switched on a fan will know.

The combination of (dry bulb) air temperatures with data for relative humidity, and wind-speed facilitates the calculation of a ‘heat index’ which is more directly related to ‘human comfort’ than any one of its constituent variables. When the heat index is high the weather will be uncomfortable and, when low, more bearable.  In a previous blog post, we examined average air temperatures and rainfall in Penang suggesting that January and February were the worst times of year to visit Penang since they are usually the hottest and October the best since it is usually coolest. The heat index plotted in the circular plot suggests a slightly different story, however.  April and May have the highest heat indices and are probably potentially the most uncomfortable times of year to visit Penang while the entire period between July and December have similar heat indices.

To produce the plot above, we have again sourced the data from NOAA's National Climatic Data Centre with the hourly dataset coming from the Integrated Surface Database.

We also provide the R code as usual (below) if you would like to produce similar plots for your area of interest.

# Load package libraries
 
library(ggplot2)
library(scales)
library(plyr)
library(dplyr)
library(reshape2)
library(circular)
library(lubridate)
library(grid)
library(zoo)
library(gridExtra)
library(weathermetrics)
library(RColorBrewer)
 
# Setting work directory
 
setwd("d:\\ClimData")
 
# Reading and reformatting raw hourly data downloaded from NCDC
 
dat<-read.table("831677118578dat.txt",header=TRUE,fill=TRUE,na.strings=c("*","**","***","****","*****","******","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"))  + 8 * 60 * 60
 
dat$year <- as.numeric(format(dat$dates,"%Y"))
 
dat$month <- as.numeric(format(dat$dates,"%m"))
 
dat$hour<-substring(as.character(dat$dates),12,13)
 
dat$min<-substr(dat$dates,15,16)
 
dat$time<-paste(dat$hour,dat$min,sep=":")
 
dat$tempc <- (dat$temp-32) * (5/9)
 
dat$dewpc <- (dat$dewp-32) * (5/9)
 
dat$tempc[dat$tempc<=10] <- NA
 
dat$tempc[dat$tempc>=40] <- NA
 
dat$dir[dat$dir == 990.0] <- NA 
 
dat$wspd <- (dat$spd)*0.44704
 
# Convert precipitation from inches to mms
 
dat$rain  <- dat$pcp24*25.4
 
# Calculate relative humidity & heat index using weathermetrics package
 
dat$rh <- dewpoint.to.humidity(t = dat$tempc, dp = dat$dewpc, temperature.metric = "celsius")
 
dat$hi <- heat.index(t = dat$tempc,rh = dat$rh,temperature.metric = "celsius",output.metric = "celsius",round = 2)
 
# Commands to reformat dates
 
dat$year <- as.numeric(as.POSIXlt(dat$dates)$year+1900)
dat$month <- as.numeric(as.POSIXlt(dat$dates)$mon+1)
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)
dat$weekday <- as.POSIXlt(dat$dates)$wday
dat$weekdayf <- factor(dat$weekday,levels=rev(0:6),labels=rev(c("Mon","Tue","Wed","Thu","Fri","Sat","Sun")),ordered=TRUE)
dat$yearmonth <- as.yearmon(dat$dates)
dat$yearmonthf <- factor(dat$yearmonth)
dat$week <- as.numeric(format(as.Date(dat$dates),"%W"))
dat$hour <- as.numeric(format(strptime(dat$dates, format = "%Y-%m-%d %H:%M"),format = "%H"))
 
dat <- ddply(dat,.(yearmonthf),transform,monthweek=1+week-min(week))
 
# Extract data from 1980 onwards
 
dat1 <- subset(dat, year >= 1980 )
 
# Summarize data for weather variables by hour and month
 
dat2 <- ddply(dat1,.(monthf,hour),summarize, wspd = mean(wspd,na.rm=T))
 
dat3 <- ddply(dat1,.(monthf,hour),summarize, temp = mean(tempc,na.rm=T))
 
dat4 <- ddply(dat1,.(monthf,hour),summarize, hi = mean(hi,na.rm=T))
 
dat6 <- ddply(dat1,.(monthf,hour),summarize, rh = mean(rh,na.rm=T))
 
## Plot Temperature circular Plot
 
p1 = ggplot(dat3, aes(x=monthf, y=hour, fill=temp)) +
  geom_tile(colour="grey70") +
  scale_fill_gradientn(colours = c("#99CCFF","#81BEF7","#FFFFBD","#FFAE63","#FF6600","#DF0101"),name="Temperature\n(Degree C)\n")+
  scale_y_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")) +
  coord_polar(theta="x") +
  ylab("HOUR OF DAY")+
  xlab("")+
  ggtitle("Temperature")+
  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"),
  panel.grid=element_blank(),
  axis.ticks=element_blank(),
  axis.text.y=element_text(size=5,colour="grey20"),
  axis.text.x=element_text(size=10,colour="grey20",face="bold"),
  plot.title = element_text(lineheight=1.2, face="bold",size = 14, colour = "grey20"),
  plot.margin = unit(c(-0.25,0.1,-1,0.25), "in"),
  legend.key.width=unit(c(0.2,0.2),"in"))
 
p1
 
 
## Plot Wind Speed circular plot
 
p2 = ggplot(dat2, aes(x=monthf, y=hour, fill=wspd)) +
  geom_tile(colour="grey70") +
  scale_fill_gradientn(colours = rev(topo.colors(7)),name="Wind Speed\n(meter/sec)\n")+
  scale_y_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")) +
  coord_polar(theta="x") +
  ylab("HOUR OF DAY")+
  xlab("")+
  ggtitle("Wind Speed")+
  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"),
  panel.grid=element_blank(),
  axis.ticks=element_blank(),
  axis.text.y=element_text(size=5,colour="grey20"),
  axis.text.x=element_text(size=10,colour="grey20",face="bold"),
  plot.title = element_text(lineheight=1.2, face="bold",size = 14, colour = "grey20"),
  plot.margin = unit(c(-0.25,0.25,-1,0.25), "in"),
  legend.key.width=unit(c(0.2,0.2),"in"))
 
p2
 
 
# Plot Heat Index circular plor
 
p3 = ggplot(dat4, aes(x=monthf, y=hour, fill=hi)) +
  geom_tile(colour="grey70") +
  scale_fill_gradientn(colours = c("#99CCFF","#81BEF7","#FFFFBD","#FFAE63","#FF6600","#DF0101"),name="Temperature\n(Degree C)\n")+
  scale_y_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")) +
  coord_polar(theta="x") +
  ylab("HOUR OF DAY")+
  xlab("")+
  ggtitle("Heat Index")+
  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"),
  panel.grid=element_blank(),
  axis.ticks=element_blank(),
  axis.text.y=element_text(size=5,colour="grey20"),
  axis.text.x=element_text(size=10,colour="grey20",face="bold"),
  plot.title = element_text(lineheight=1.2, face="bold",size = 14, colour = "grey20"),
  plot.margin = unit(c(-0.5,0.25,-0.25,0.25), "in"),
  legend.key.width=unit(c(0.2,0.2),"in"))
 
p3
 
 
# Plot Relative Humidity circular plor
 
col<-brewer.pal(11,"Spectral")
 
p4 = ggplot(dat6, aes(x=monthf, y=hour, fill=rh)) +
  geom_tile(colour="grey70") +
  scale_fill_gradientn(colours=col,name="Relative Humidity\n(%)\n")+
  scale_y_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")) +
  coord_polar(theta="x") +
  ylab("HOUR OF DAY")+
  xlab("")+
  ggtitle("Relative Humidity")+
  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"),
  panel.grid=element_blank(),
  axis.ticks=element_blank(),
  axis.text.y=element_text(size=5,colour="grey20"),
  axis.text.x=element_text(size=10,colour="grey20",face="bold"),
  plot.title = element_text(lineheight=1.2, face="bold",size = 14, colour = "grey20"),
  plot.margin = unit(c(-0.5,0.1,-0.25,0.25), "in"),
  legend.key.width=unit(c(0.2,0.2),"in"))
 
p4
 
## Plot and export to png file
 
png(file="Weather_Variables_Circular_Plot.png",width = 12, height = 10, units = "in",
    bg = "white", res = 400, family = "", restoreConsole = TRUE,
    type = "cairo")
 
grid.arrange(p1,p2,p4,p3,nrow=2,main=textGrob("\nWeather Variables by Month and Hour\n(Bayan Lepas Weather Station)",
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