[GRASS-user] find rivers, river paths, basin area, cumulative basin area along the path

[Pankaj Kr Sharma]

Dear list,

I am interested in finding basin area limited to say 300sq km.

1. Easier approach: is As suggested by Mr. M Metz earlier. It's quite close.

Use SFD in r.watershed. Keep threshold in terms of number of cells equal
to  the area required. Use a raster dem 1 cells (r.grow) bigger than the
basin boundary to get positive accumulation values.

Then run r.stream.extract. Keep threshold equal to the number used above.
Find starting points . Use layer 2 option and type_code=0.

These points are having the catchment equal to the threshold. (at least, as
far as I have checked. I am still on job.)

The river starts forming when the given threshold is reached.


2. Not so easy approach:

--------------------------------------------------------------------------------------------------
echo "delete  from cat_temp_done;" > file.sql;
db.execute input=file.sql;

r.watershed -s --overwrite -m memory=12000
elevation=ybdem at PERMANENTthreshold=12345 accumulation=step1_acc
drainage=step1_fdir
basin=step1_basin stream=step1_stream

r.stream.extract --overwrite  memory=12000
elevation=ybdem at PERMANENTaccumulation=step1_acc threshold=12345
stream_length=2
stream_rast=step1_stream_rse stream_vect=step1_stream_rse
direction=step1_stream_fdir

v.db.addcolumn map=step1_stream_rse at PERMANENT layer=1 columns="basin_area
double precision,basin_area_point double precision"

v.db.addtable map=step1_stream_rse at PERMANENT layer=2


#find start points

v.extract -t --verbose input=step1_stream_rse at PERMANENT layer=2 type=point
cats=0 output=zero_points_delete
v.category --overwrite --verbose
input=zero_points_delete at PERMANENToutput=zero_points_no_cat_delete
option=del
v.category --overwrite --verbose
input=zero_points_no_cat_delete at PERMANENToutput=zero_points_with_cat_delete
option=add
v.select  --overwrite --verbose ainput=step1_stream_rse at PERMANENT alayer=1
atype=point binput=zero_points_with_cat_delete at PERMANENT blayer=1
btype=point output=zero_points at PERMANENT operator=equals

#find outlets

v.extract -t --verbose input=step1_stream_rse at PERMANENT layer=2 type=point
cats=2 output=outlet_points_delete
v.category --overwrite --verbose
input=outlet_points_delete at PERMANENToutput=outlet_points_no_cat_delete
option=del
v.category --overwrite --verbose
input=outlet_points_no_cat_delete at PERMANENToutput=outlet_points_with_cat_delete
option=add
v.select  --overwrite --verbose ainput=step1_stream_rse at PERMANENT alayer=1
atype=point binput=outlet_points_with_cat_delete at PERMANENT blayer=1
btype=point output=outlet_points at PERMANENT operator=equals



start_cats=`echo "select cat from zero_points_step1_stream_rse;" |
db.select -c;`
end_cat=`echo "select cat from outlet_points_step1_stream_rse;" | db.select
-c;`
for i in $start_cats
do
for j in $end_cat
do
echo "i:$i, j:$j"
v.net.salesman --overwrite
input=step1_stream_rse at PERMANENToutput=river_path_$i type=line
alayer=1 nlayer=1 ccats=$i,$j
v.db.addtable map=river_path_$i at PERMANENT columns="cat integer, basin_area
double precision, basin_area_point double precision"
path_line_cats=`echo "select cat from river_path_$i;" | db.select -c;`
prev_area_value=0;
for x in $path_line_cats
do
echo "i:$i, j:$j, x: $x"
r.stream.basins
dirs=step1_stream_fdir at PERMANENTstreams=step1_stream_rse at PERMANENTcats=$x
basins=river_path_basins_$i_$x
area_value=`r.stats -a -n --verbose
input=river_path_basins_$i_$x at PERMANENT| awk '{print $2}'`
echo "Update  river_path_$i set basin_area=$area_value where cat=$x;" >
file.sql;
db.execute input=file.sql;
area_value_point=$area_value-$prev_area_value
echo "Update  river_path_$i set basin_area_point=$area_value_point where
cat=$x;" > file.sql;
db.execute input=file.sql;
$prev_area_value=$area_value
done
line_cats=`echo "select cat from river_path_$i and cat not in (select cat
from cat_temp_done);" | db.select -c;`
for k in $line_cats
do
echo "i:$i, j:$j, x: $k"
r.stream.basins
dirs=step1_stream_fdir at PERMANENTstreams=step1_stream_rse at PERMANENTcats=$k
basins=river_path_basins_$i_$k
area_value=`r.stats -a -n --verbose
input=river_path_basins_$i_$k at PERMANENT| awk '{print $2}'`
echo "Update  step1_stream_rse set basin_area=$area_value where cat=$k;" >
file.sql;
db.execute input=file.sql;
area_value_point=$area_value-$prev_area_value
echo "Update  step1_stream_rse set basin_area_point=$area_value_point where
cat=$k;" > file.sql;
db.execute input=file.sql;
$prev_area_value=$area_value
echo "insert into gis_schema.cat_temp_done (cat) values($k);" > file.sql;
db.execute input=file.sql;
done
done
done

___________________________________________________________________________________________________--

The above codes will neatly arrange the required values in the columns.
Which I shall later use on as label on points on river path.

3. Not so easy -- Not so easy approach

Replace the one line code of v.net.salesman used above with this:

--------------------------------------------------------------------------------------------------------------------
elcat=1972;
echo "delete  from  work.cat_temp;" > file.sql;
db.execute input=file.sql;
echo "delete  from  work.cat_temp_done;" > file.sql;
db.execute input=file.sql;
echo "delete  from  work.yrse100sqkm_stream_path1;" > file.sql;
db.execute input=file.sql;
while [  "$elcat" ]
do
cats_to_write=`echo "select cat from work.yrse100sqkm_stream where
end_long=(select start_long from work.yrse100sqkm_stream where cat=$elcat)
and end_lat=(select start_lat from work.yrse100sqkm_stream where
cat=$elcat);" | db.select -c;`
echo "delete  from  work.cat_temp where cat= $elcat;" > file.sql;
db.execute input=file.sql;
echo "cats_to_write: $cats_to_write";
if [ ! -n "$cats_to_write" ]
then
echo "Update  work.yrse100sqkm_stream set typecode=0 where cat=$elcat;" >
file.sql;
db.execute input=file.sql;
else
echo "Update  work.yrse100sqkm_stream set typecode=1 where cat=$elcat;" >
file.sql;
db.execute input=file.sql;
for k in $cats_to_write
do
echo "insert into work.cat_temp (cat) values($k);" > file.sql;
db.execute input=file.sql;
prev_path=9999999;
prev_path=`echo "select cat from work.yrse100sqkm_stream_path1 where
path_cat=$elcat; " | db.select  -c;`
for p in $prev_path
do
echo "p:k: $p,$k";
echo "insert into work.yrse100sqkm_stream_path1 (cat,path_cat)
values($p,$k);" > file.sql;
db.execute input=file.sql;
done
echo "insert into work.yrse100sqkm_stream_path1 (cat,path_cat)
values($elcat,$k);" > file.sql;
db.execute input=file.sql;
done
fi
elcat=`echo "select  cat from work.cat_temp limit 1;" | db.select -c;`
done
------------------------------------------------------------------------------------------------------------------------------------------------------------------------

or with this


--------------------------------------------------------------------------------------------------------------------------------------------------------------------------

elcat=1972;
echo "delete  from  work.cat_temp;" > file.sql;
db.execute input=file.sql;
echo "delete  from  work.cat_temp_done;" > file.sql;
db.execute input=file.sql;
while [ $elcat ]
do
cats=`echo "select cat from work.yrse100sqkm_stream where end_long=(select
start_long from work.yrse100sqkm_stream where cat=$elcat) and
end_lat=(select start_lat from work.yrse100sqkm_stream where cat=$elcat)
and cat <> $elcat; " | db.select -c;`
if [ ! -n "$cats" ]
     then
    echo "insert into work.river_startline(cat) values($elcat);" > file.sql;
    db.execute input=file.sql;
fi
echo "insert into work.cat_temp_done (cat) values($elcat);" > file.sql;
db.execute input=file.sql;
cats_to_write=`echo "select cat from work.yrse100sqkm_stream where
end_long=(select start_long from work.yrse100sqkm_stream where cat=$elcat)
and end_lat=(select start_lat from work.yrse100sqkm_stream where
cat=$elcat)  and cat not in (select cat from work.cat_temp_done); " |
db.select -c;`
echo "delete  from  work.cat_temp where cat= $elcat;" > file.sql;
db.execute input=file.sql;
for k in $cats_to_write
do
echo "insert into work.cat_temp (cat) values($k);" > file.sql;
db.execute input=file.sql;
done
elcat=`echo "select  cat from work.cat_temp limit 1;" | db.select -c;`
done

and then

echo "delete  from  yrse100sqkm_stream_path;" > file.sql;
db.execute input=file.sql;
elcats=`echo "select cat from work.river_startline;" | db.select -c;`
for i in $elcats
do
cats=$i
declare -i order;
order=1;
while [ $cats ]
do
echo "insert into work.yrse100sqkm_stream_path (cat,path_cat,orderno)
values($cats,$i,$order);" > file.sql;
db.execute input=file.sql;
echo "cat: $cats,path: $i,order: $order";
cats=`echo "select cat  from work.yrse100sqkm_stream  where
start_long=(select end_long from work.yrse100sqkm_stream where cat=$cats)
and start_lat=(select end_lat from work.yrse100sqkm_stream where cat=$cats)
limit 1; " | db.select -c;`
order=order+1;
done
done


_________________________________________________________________________________

Definitely, knowledge of  python would have made things easier.

Thanks to all for the help.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20111203/c66b5237/attachment.html