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/

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

No comments:

Post a Comment