Creating road network files for input into DHSVM

/**** Start arc/info

plane% setup owclient arcinfo
plane% arc
arc: &workspace /nfs/meter/usr1/lxb/arcinfo/data

/**** Import ascii dem file to grid:

Add arc/info header to ascii dem file:
ncols 110
nrows 115
xllcorner 546700
y11corner 5173900
cellsize 30
nodata_value 0

arc: &station 9999

usage: ASCIIGRID <in_ascii_file> <out_grid> {INT|FLOAT}
arc: asciigrid dem_filename dem01 INT

/**** Define projection

arc: projectdefine GRID dem01
project: projection utm 
project: units meters
project: zone 10
project: xshift 0
project: yshift 0
project: parameters
project: end
arc: describe dem01

/**** Display grid

arc: grid
grid: shadeset rainbow
grid: mapextent dem01
grid: grids dem01 # linear
 
/**** Locate sinks and fill dem (in grid)
 
grid: flowdirgrid = flowdirection(dem01)
grid: sinkgrid = SINK(flowdirgrid)
grid: grids sinkgrid
grid: FILL dem01 filleddem SINK

/**** Prepare road coverage

Start with a road coverage with individual arcs as long as possible.  Nodes
are only located at arc (and class) intersections.  Each segment should have
a class associated with it.  The class designates the following
charactersistics: ditch width, ditch depth, manning's n, cutbank slope, and
crown type.  (At this point the crown type is not used by DHSVM).

/**** Add class category to the road network:

usage:  ADDITEM <in_info_file> <out_info_file> <item_name> <item_width>
		<output_width> <item_type> {decimal_places} {start_item}

arc: additem rdcov01.aat rdcov01.aat CLASS 1 1 i


/**** IF YOU DO NOT HAVE AN ACTUAL CULVERT LOCATION FILE

Use ARCINTERSECT to prepare a point coverage containing points where streams
and roads intersect.

arc: &run ARCINTERSECT strcov01 rdcov01 culvcov01 


/**** IF YOU DO HAVE AN ACTUAL CULVERT LOCATION FILE

Plot roads, streams and culvert locations to be sure there is a culvert at
each stream crossing.  Edit streams and culverts accordingly.  It is better 
to move culverts to match stream crossing locations than to move streams because
you need to keep the stream in the topographic hollow of the DEM.  Otherwise, 
the imposed stream channel will not capture all of the simulated subsurface flow.

arc: arctools

select <Edit Tools>
select <File> -> <Coverage: Open>


/**** Add a type attribute to the culvert file.  Set to 'S' for sink.

arc: additem culvcov01.pat culvcov01.pat TYPE 1 1 c
arc: ap
arcplot: aselect culvcov01 point
arcplot: calculate culvcov01 point type = 'S'

/**** Split road coverage at known culvert locations.

arc: copy rdcov01 rdcov02
arc: &run NODESFROMPOINTS rdcov02 culvcov01 5.0

Use pointnode to transfer 'type' from the culvert coverage to the road
network nodes.

usage: POINTNODE <point_cover> <node_cover> {search_radius}

arc: pointnode culvcov01 rdcov02 5.0


/**** Run first three scripts

arc: &run roadelevation rdcov02 filleddem

arc: &run roadslope rdcov02
arc: &run roaddivide rdcov02
arc: list rdcov02.div

The listed points indicate the number of 'divides' and 'sinks' you have
in your stream network.  Create a point coverage to compare against known 
culvert locations:

/**** Create point coverage from point event database

usage: EVENTSOURCE add point <source_name> <table_name> {database} 
{relate_type} {route_key_item} {event_key_item} {measure_item}

arc: eventsource add point rdsource rdcov02.div info linear elevation#
elevation# measure

** NOTE:  source name can be up to 8 characters long

usage: EVENTPOINT <in_cover> <in_route_system> <event_source> 
<out_cover>

arc: eventpoint rdcov02 elevation rdsource rddivs01


/**** Compare sink locations with known culvert point coverage.

arc: arctools

select <Edit Tools>
select <File> -> <Open: Coverage> -> <rddivs01>
select <Display> -> <Background env: General...> -> culvcov01		
						    rdcov01
						    strcov01
Delete sinks which repeat culvert locations

/**** This is sort of the sensitive part in derivation of the road network 
and you can do a variety of things depending on your knowledge of the actual 
road drainage.  If you do not have actual culvert locations, break the road 
arcs at divide and sink locations.  This is equivalent to assuming that a 
culvert exists at all areas of convergence:

arc: copy rdcov02 rdcov03

arc: &run nodesfrompoints rdcov03 divpoints 5.0

/**** If you have some idea of the actual drainage patterns of your roads, 
you will want to check the road drainage direction derived from the DEM.  Plot 
road drainage direction as follows, the drainage direction can be changed by 
switching the 'sign' attribute:

arc: &run roaddirection rdcov02 direct01
arc: arcplot
arcplot: linesize 0.0
arcplot: linecolor blue
arcplot: linepattern 0
arcplot: lineput 1
arcplot: linesymbol 1
arcplot: arcs rdcov02
arcplot: arrowsize 0.07
arcplot: arrowtype double solid
arcplot: arcarrows direct01 sign


Use pointnode to transfer 'type' from the culvert coverage to the road
network nodes.

usage: POINTNODE <point_cover> <node_cover> {search_radius}

arc: pointnode rddivs01 rdcov03 5.0


/**** Create network input file

arc: &run roadorder rdcov03
arc: &run roadnetwork rdcov03 rd_net.dat class


/**** Create row/column map
(This is the section which is affected by the grid error.  Create a copy of
the stream coverage so you can return to this point when arc/info is fixed.)

arc: copy rdcov03 rdcov04
arc: grid
grid: slope01 = slope(filleddem)
grid: aspect01 = aspect(filleddem)
grid: q

arc: &run rowcolmap aspect01 slope01 rowcol01

/**** Create channel map input file

&run roadmap rdcov04 rowcol01 rdmap01

&run roadaspect rdmap01 filleddem


/**** Create road class look-up table:

arcedit: create roads.dat info
	 class	2, 2, i
	 width  4, 4, f, 2
	 slope	4, 4, f, 2
	 ddepth 8, 8, f, 4
	 dwidth  8, 8, f, 4
	 rdslope 12, 12, c 

arcedit: q
arc: info
input user name: arc
info: select roads.dat
info: add from /nfs/meter/usr1/lxb/arcinfo/test/rd_class.csv
info: q stop

/** Note: Must use complete address with add from command.

Example rd_class.csv:

1, 3.1,  0.5, 0.089, 0.219
2, 3.1,  0.5, 0.321, 1.224 
3, 3.1,  23.3, 0.089, 0.219 
4, 3.1,  23.3, 0.089, 0.910
5, 3.1,  7.5, 0.089, 0.910
6, 3.1,  7.5, 0.321, 1.224 
7, 3.1,  23.3, 0.200, 1.224 
8, 3.1,  23.3, 0.200, 0.910


(Columns correspond to: road class, road width, cutbank slope,ditch depth,
ditch width.)

arc: relate add; rd_class; roads.dat; info; class; class; ordered; rw;;

/**** Estimate road segment cut width and depth

arc: &run roadcut rdmap01 rd_class//width rd_class//dwidth rd_class//ddepth
rd_class//slope

/**** Create stream map output file:

arc: &run roadmapfile rdmap01 rdcov04 rd_map.dat
