Remote sensing rainfall product, part I: CMORPH

As hydrologist, there is no way to escape from remote sensing. Especially if we are working in developing countries where there is no (or few) ground measurements. Since I started working on the Upper Blue Nile basin (UBN) last month, I slowly learning how to harvest the major remote sensing products. In the area of remote sensing research, there is huge resources that can be used in hydrology. Now, let me start from the rainfall products. I will post how I will analyse other products, and for now, this post is how I tried to use one rainfall product in R: CMORPH products. CMORPH is actually high temporal and spatial resolution product. First load the necessay library

 library(sp)
library(raster)
library(ncdf4)
library(rgdal)
library(rasterVis)
library(latticeExtra)
library(maps)
library(mapdata)
library(lattice)

Gather data (or download it on your Pc) and look at the netCDF file structure

There should be direct way to downloading methods inside R, however, I download CMORPH from NOAA webpage manualy. So I need to read the files inside my folder

setwd("~/UBN_rainfall/CMORPH/Project/") 
fl<-list.files("~/UBN_rainfall/CMORPH/Project/", pattern='*.nc') 

Create a raster brick

We can create raster brick (layer of raster, each corresponds to one time series data). It is very huge dataset, and it took a while to brick all the 4000 file (4000 days of images) with 8 layer(which if for every three hours) together. I chnaged the no value into -9999

CMORPH<-brick(stack(fl)) # I changed  raster stack into brick so that it will be faster for analysis
NAvalue(CMORPH) <- -9999
#dim(CMORPH)

The area downloaded is a little larger than I need so here I have to extract to a particular basin extent,zoom in to the extent of the study area, in my case UBN basin.

extent_UBN <- extent(34.5,39.75,7.75,12.75)
CMORPH_UBN <- crop(CMORPH, extent_UBN)

Visualization

we can do many things on visualizing this brick raster, like animation, but lets do one thing now. select the most rainy time steps and plot it

RF<-cellStats(CMORPH_UBN,stat=mean) # rainfall mean over the basin
cmorph_rainyday<-which(RF==max(RF[1:32000],na.rm=T)) # select max rainfall time steps
Rainyday <-raster(CMORPH_UBN,layer=cmorph_rainyday) # use to single out the rainy time step
levelplot(Rainyday)# plot for the rainy time-step

Lets improve the visualization using new color palette, and lets use lattice (trellis) plot methods. I also want to overlay the shapefile of the basin

UBNshape<- readOGR(dsn="~/UBN_20136/Abbay/cell", layer="mask")  
proj4string(UBNshape) # describes data’s current coordinate reference system # to change to correct projection: 
UBNshape <- spTransform(UBNshape, CRS("+proj=longlat +datum=WGS84"))   #### change to WGS84 

palette_new <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"),  space = "rgb")

spplot(Rainyday, col.regions=palette_new,
       colorkey=list(height=0.3),
       sp.layout=list("sp.polygons", UBNshape, lwd=1))

Now let’s create sptial point object (lets say meteorological station) and create a time series object.

p <- cbind(35.5,11.5) 
sp <- SpatialPoints(p) 
CMORPH_P1 <- extract(CMORPH_UBN, sp)

#create the dataframe with time and CMORPH rainfall column 
CMORPH_p<-data.frame(time=seq(ISOdate(2002,12,7,0), ISOdate(2013,11,18,21), by = "3 hours"),CMORPH=as.vector(CMORPH_P1))

Usually I like to use high level graphics such as ggplot and lattice. lets try the lattice graphics here

xyplot(CMORPH ~ time| equal.count(as.numeric(time),
                                    7, overlap = 0.1),
   data=CMORPH_p, strip = FALSE, type=c("g", "s"),
   scales = list(x = list(relation = "sliced", axs = "i"),
        y = list(relation = "free", alternating = FALSE)))