Mar 28, 2015

Remote sensing rainfall product, part I: CMORPH

Remote sensing in hydrology: 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
plot of chunk unnamed-chunk-5
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))
plot of chunk unnamed-chunk-6
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)))
plot of chunk unnamed-chunk-8