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 <- data.frame(Id=seq(1,length.out=length(x)),X=x,Y=y)
#ddTable <- data.frame(Id=seq(1,length.out=length(x)))
#ddShapefile <- convert.to.shapefile(dd, ddTable, "Id", 3)
#write.shapefile(ddShapefile, "test", arcgis=T)
13.727896
100.524123
