Plotting Rainfall Map with ggplot2

Tags

In previous post , we post how to retrieve and read the TRMM rainfall format. This post will show you how mapping the rainfall grid to the region boundary.

trmm_map
library(dplyr)
library(ggplot2)
library(reshape2)
library(gridExtra)
library(maps)

inputfile ="3B42_daily.2014.01.02.7.csv"
infile = file(inputfile,"r")
rain = tbl_df(read.csv(infile,stringsAsFactors=FALSE))

#==== CROP BOUNDARY AREA and Data ============
max.lat = 20.75
min.lat = 5.5
max.lon = 106
min.lon = 97.00

rain_thai = rain %>%
 filter(LAT=min.lat,LAT=max.lat,LON=min.lon,LON=max.lon)

#========= Plot Grids Data ============
p = ggplot(data = rain_thai,
    aes(x=LON, y=LAT,group=RAIN,fill=RAIN))
p = p+geom_tile(aes(fill=RAIN)) 
+scale_fill_gradient(name="Rain(mm/day)",low="white",high="blue")
p = p+ggtitle("Daily Rainfall 2014-01-02")

# ====== THAILAND Boundary ====================

world = map_data(world2)
world = subset(world,region=="Thailand")
map = geom_path(size=1,aes(x=long,y=lat,group = group),data=world)
p1 = ggplot(rain_thai,aes(x=LON,y=LAT))
p1 = p1+geom_tile(aes(fill = RAIN))
+scale_fill_gradient(name="Rain(mm/day)",low="white",high="blue")
p1 = p1+ ggtitle("Daily Rainfall 2014-01-02")
p1 = p1+theme_bw()
p1 = p1+map

grid.arrange(p, p1+map, ncol=2, main="examples")
Advertisements

Reading TRMM (Tropical Rainfall Measuring Mission) Rainfall Binary Data

Tags

Tropical Rainfall Measuring Mission (TRMM) is providing the global rainfall and climate. It’s spatial resolution 0.25-degree by 0.25-degree and spatial coverage extend from 50 degrees south to 50 degrees north latitude. You can download the daily TRMM rainfall from this link TRMM Daily 3B42 Product. The TRMM manual give us the fortran and matlab code for reading the trmm 3B42 daily product but we want to use R. So the below code is modified from matlab to R.

inputfile ="3B42_daily.2014.01.02.7.bin"
binfile = file(inputfile,"rb")
rawdata = readBin(binfile, numeric(),size=4,n=576000,endian = "big") 
close(binfile)

dailyRain.trmm =data.frame(Lat=numeric(576000),Lon=numeric(576000),Rain=numeric(576000))
i=1
#========= The below code is modified from       ===========
#========= the reading trmm 3B42 daily manual    ===========
#========= for lat-lon correction in matlab code ===========
 
for( i_lat in 1:400)
{
 for(j_lon in 1:1440)
 {
    lat = 49.875 - 0.25*(i_lat - 1)
    if (j_lon <= 720)
    {
       lon = 0.125 + 0.25*(j_lon - 1) 
    }else{
       lon = 0.125 + 0.25*(j_lon - 1) - 360.0 
    }
    dailyRain.trmm$Lat[i] = lat
    dailyRain.trmm$Lon[i] = lon
    dailyRain.trmm$Rain[i] = rawdata[i] 
    i=i+1
 }
}

The running time of binary reading script (line 1-4) is vary fast but the lat-lon correction is very long. Any comments will be thankful.

Role of R in Regional Water-Economic Model

Tags


Since I have seen the water demand forecasting model based on economic activities/parameters by Input-Output Table (I-O Table), I look forward the R package to solve the I-O Table. Then I see the R|E|A|L (Regional Economic Analysis Laboratory) developed by THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN.

This software is freeware based on python as interface building and R as main matrix calculation. When you downloaded and installed in your PC, there will providing you the source code to matrix calculation in economic analysis such as R.A.S method. So I think it is useful for anyone who interesting in R and water-economic analysis.

Create XY coordinate at the specified distance in the river stream / polyline

Tags

, ,


In using the hydrodynamic software such as InfoWorks RS for flood simulation in the river basin, we hit the problems in creating the river sections and stream link between adjacent sections in the model especially for the large networks model. The general problems are that you have to copy the river section data from excel file, paste data into the model at each section by section, drawing the link between upstream and downstream section and sometime finding the survey cross section point by measuring the distance along the river stream too. However, some software can import the river network from text file as specified format data. So we wrote the some scripts to create the text file input format for InfoWorks RS as following steps:
1. Reading coordinate (x,y) of point in polyline from shape file of river       stream line
2. Generate the new coordinate point at specified distance as survey       point
3. Reading the cross section data file
4. Generate coordinate of each cross section data point where the           new coordinate points were set to be the coordinate of river bed for     each section
5. Import coordinate of river section data into the model
 Because of the input format for coordinate of river section is different in each program. We will only present script for use in step 1-2.

# THIS EXAMPLE CODE IS USED FOR CREATING THE INTERPOLATE POINT/ INTERESTING POINT AT SPECIFIED DISTANCE
# FROM STARTING POINT IN RIVER STREAM
# WRITTEN BY WONGWATANA SOMBUNYING
library("sp")
library("spatial")
library("shapefiles")
#setwd("D:5_R\\49_IW\\")
xy = read.shapefile(".\\inp_shp\\line")$shp$shp[[1]]$points
x  = xy$X
y  = xy$Y

# Create dataframe and segment length
shp.dataframe = as.data.frame(cbind(x,y,d=0,sumd=0))
sumd=0
for (i in 1:(length(x)-1))
{
      dx = x[i]-x[i+1]
      dy = y[i]-y[i+1] 
	  dd = dx^2+dy^2
	  d  = dd^0.5
	  sumd=sumd+d	
	  shp.dataframe$d[i+1]=d
	  shp.dataframe$sumd[i+1]=sumd	  
}

# Create the specified point
point.n = seq(0,max(shp.dataframe$sumd),by=10)    # 10 m interval
x.n     = array(dim=(length(point.n)-1))
y.n     = array(dim=(length(point.n)-1))
for ( j in 2:length(point.n))
{
    for (i in 2:(length(x)))
	{
	   if (shp.dataframe$sumd[i] > point.n[j])
		{
		  x1=shp.dataframe$x[i-1]
		  x2=shp.dataframe$x[i]
		  y1=shp.dataframe$y[i-1]
		  y2=shp.dataframe$y[i]
		  dx = x2-x1 
		  dy = y2-y1
		  dD = (point.n[j]-shp.dataframe$sumd[i-1])/shp.dataframe$d[i]
		  x.n[j-1] = x1+dD*dx
		  y.n[j-1] = y1+dD*dy
          break
		}
	}
}
n.dataframe = as.data.frame(cbind(x=x.n,y=y.n,d=0,sumd=point.n[-1]))
new.point = merge(n.dataframe,shp.dataframe,all=T)
new.point = new.point[order(new.point[,4]),]

fout=file(".\\out_shp\\new_point.txt")
open(fout,open="w")
write(sprintf("%s,%s,%s","x","y","len"),file=fout,sep="\n",append=TRUE)
write(sprintf("%f,%f,%f",new.point$x,new.point$y,new.point$sumd),file=fout,sep="\n",append=TRUE)
close(fout)

# Create PolyLine
#dd &lt;- data.frame(Id=seq(1,length.out=length(x)),X=x,Y=y)
#ddTable &lt;- data.frame(Id=seq(1,length.out=length(x)))
#ddShapefile &lt;- convert.to.shapefile(dd, ddTable, &quot;Id&quot;, 3)
#write.shapefile(ddShapefile, &quot;test&quot;, arcgis=T)

Bivariate Flood Frequency Analysis by R Copula Function

Tags

, ,


The joints distribution frequency analysis by copula family is presented with R code  (example#code). The  data used in this analysis were  flood peak, flood volume and flood duration derived from the historical flood. In the example code show the required library package, how to estimate function’s parameters, how to determine the joint probability function from observed data, how to generate/ synthesis data from copula function and how to create the contour plot.
copula_floodNote : the input format has 4 colums : YEAR, Flood Volume, Flood Duration and Flood Peak