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 structureThere 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 brickWe 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
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.
CMORPH<-brick(stack(fl)) # I changed raster stack into brick so that it will be faster for analysis NAvalue(CMORPH) <- -9999 #dim(CMORPH)
extent_UBN <- extent(34.5,39.75,7.75,12.75) CMORPH_UBN <- crop(CMORPH, extent_UBN)
Visualizationwe 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.
Usually I like to use high level graphics such as ggplot and lattice. lets try the lattice graphics here
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))
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)))