This blog is run by Jason Jon Benedict and Doug Beare to share insights and developments on open source software that can be used to analyze patterns and trends in in all types of data from the natural world. Jason currently works as a geospatial professional based in Malaysia. Doug lives in the United Kingdom and is currently Director of Globefish Consultancy Services which provides scientific advice to organisations that currently include the STECF [Scientific, Technical and Economic Committee for Fisheries, https://stecf.jrc.europe.eu/] and ICCAT, https://www.iccat.int/en/

Showing posts with label Malaysia. Show all posts
Showing posts with label Malaysia. Show all posts

Monday, 13 April 2015

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

Friday, 30 January 2015

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 17th and 24th 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

Thursday, 28 August 2014

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 26th 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