GRASS Comand Sequence

extract_morphometric_parameters

# Watershed SFDo 
r.watershed --o -a elevation=srtm.clip accumulation=accumulation drainage=drainage convergence=5 threshold=123
#Without SFDo
#r.watershed elevation=srtm.clip threshold=555 accumulation=accumulation1 drainage=drainage1 basin=ba500

#r.watershed elevation=filled_clip_srtm_gramame threshold=123 accumulation=accum drainage= dirs basin=ba1000 stream=streams

#Stream extraction:
r.stream.extract --o elevation=srtm.clip threshold=123 d8cut=infinity mexp=0 stream_rast=streams direction=drain_dir stream_vect=v_streams

#Delineation of basin:
r.stream.basins --o coors=280666.416304,9190588.658405 dir=drain_dir basins=basin

#Mask and cropping
r.mapcalc expression='mde_basin=clip_srtm_gramame*basin'
r.mapcalc expression='accum_basin=accumulation*basin'
r.mapcalc expression='drain_basin=drainage*basin'
r.mapcalc expression='streams_basin=streams*basin'
r.mapcalc expression='drain_dir_basin=drain_dir*basin'

#Extract Network River
r.thin --o input=streams_basin output=streams_basin_thin

r.to.vect --o input=streams_basin_thin output=v_network@permanent type=line

#Creation of slope and aspect maps
r.slope.aspect --o elevation=mde_basin slope=slope_basin aspect=asp_basin

#Basin mask (vector)
#Raster to vector
r.to.vect --o -sv input=basin output=basin feature=area

# Add two columns to the table: area and perimeter                               
v.db.addcol map=basin columns='area double precision'
v.db.addcol map=basin columns='perimeter double precision'

#Populate perimeter and area columns                                 
v.to.db --o map=basin type='line,boundary' layer=1 qlayer=1 option='perimeter' units='kilometers' columns='perimeter'
v.to.db --o map=basin type='line,boundary' layer=1 qlayer=1 option='area' units='kilometers' columns='area'

perimeter=`echo "SELECT perimeter FROM basin" | db.select -c`
area=`echo "SELECT area FROM basin" | db.select -c`

#Creation of order maps: strahler, horton, hack, shreeve
r.stream.order --o stream=streams_basin dir=drain_dir_basin strahler=r_strahler shreve=r_shreve horton=r_horton hack=r_hack

# Distance to outlet
touch outlet.tmp
echo "280802.291361|9190692.838850" > outlet.tmp

v.in.ascii --o input=outlet.tmp output=v_outlet x=1 y=2

v.to.rast --o input=v_outlet output=r_outlet use=cat type=point layer=1 value=1 rows=4096

r.stream.distance --o -o stream=r_outlet dir=drain_dir_basin distance=r_distance

#Ipsographic curve
r.ipso -ab map=mde_basin image=/home/marcello/Desktop/ipso

#Width Function
r.wf map=r_distance image=/home/marcello/Desktop/wf

#Creation of map of hillslope distance to river network
r.stream.distance --o stream=streams_basin dir=drain_dir_basin dem=mde_basin distance=r_hillslope_distance

#Mean elevation
r.average --o base=mde_basin cover=mde_basin output=r_height_average
r.info -r map=r_height_average

eval `r.info -r map=r_height_average`

# dúvida...
mean_elev=$min

r.mapcalc expression='r_aspect_mod=if(asp_basin > 90, 450 - asp_basin, 90 - asp_basin)'

#Centroid and mean slope
r.volume --o data=slope_basin clump=basin centroids=v_centroid1 > volume.txt

#em seguida...
 # baricenter_slope_baricenter = baricenter_slope_baricenter.split()
 # mean_slope = baricenter_slope_baricenter[28]

mean_slope=`head -n 6 volume.txt | tail -1 | awk '{print $2}'`

#Rectangle containing basin
basin_east=`head -n 6 volume.txt | tail -1 | awk '{print $5}'`
basin_north=`head -n 6 volume.txt | tail -1 | awk '{print $6}'`

eval `g.region -m vect=basin@user`

basin_resolution=$nsres

	# x_massimo = float(dict_region_basin['n']) + (basin_resolution * 10)
	# x_minimo = float(dict_region_basin['w']) - (basin_resolution * 10)
	# y_massimo = float(dict_region_basin['e']) + (basin_resolution * 10)
	# y_minimo = float(dict_region_basin['s']) - (basin_resolution * 10)
	# nw = dict_region_basin['w'], dict_region_basin['n'] 
	# se = dict_region_basin['e'], dict_region_basin['s'] 

	# #Directing vector 
	# delta_x = abs(float(basin_east) - east)
	# delta_y = abs(float(basin_north) - north)
	# L_orienting_vect = math.sqrt((delta_x**2)+(delta_y**2)) / 1000

	# #Prevalent orientation 
	# prevalent_orientation = math.atan(delta_y/delta_x)

	# #Compactness coefficient 
	# C_comp = perimeter_basin / ( 2 * math.sqrt( area_basin / math.pi))

	# #Circularity ratio 
	# R_c = ( 4 * math.pi * area_basin ) / (perimeter_basin **2)

# verificar depois uma forma de como realizar os cálculos... (shell, BD, PHP)

# Mainchannel
r.mapcalc 'r_mainchannel=if($r_hack==1,1,null())'

r.thin --o input=r_mainchannel output=r_mainchannel_thin

r.to.vec --o input=r_mainchannel_thin output=v_mainchannel feature=line

v.db.addcol map=v_mainchannel columns='length double precision'

#OBS: v.to.db -> pode ser utilizado para extrair outros parâmetros (ver o help dele depois)
v.to.db --o map=v_mainchannel type='line,boundary' layer=1 qlayer=1 option='length' units='kilometers' columns='length'

length_km=`echo "SELECT length FROM v_mainchannel" | db.select -c`

# Topological Diameter
# ...


# Mean slope of mainchannel
# ...