----------
X-Sun-Data-Type: text
X-Sun-Data-Description: text
X-Sun-Data-Name: text
X-Sun-Content-Lines: 48

Hi,

There is a contributed program called r.length which should do what you
want, given you use some kind of a mask.  You would have to mask because 
it will measure the _total_ length of each class of line.  It also will only work
on raster data. I have included the program as an attachment in case you
wanted to take a look.  Just do a "chmod +x r.length" first to get it to execute.

-Eric Lambert

=========================================================================
Eric Lambert					
ENE Division	         	OR		GIS Labratory/Geography
P.O. Box 9005					318 Davenport Hall
Champaign, IL 61826-9005			Urbana, IL 61801

E-mail:lambert@gis.uiuc.edu			Lab#: (217) 333-5919
=========================================================================



<> From grass-lists-owner@max.cecer.army.mil Thu Aug  4 10:35:52 1994
<> Date: Thu, 4 Aug 94 16:29 BST-1
<> From: nra@cix.compulink.co.uk (N R A  )
<> Subject: Length of lines
<> Sender: grass-lists-owner@max.cecer.army.mil
<> Reply-To: grassu-list@max.cecer.army.mil
<> To: grassu-list@max.cecer.army.mil
<> Reply-To: nra@cix.compulink.co.uk
<> Content-Length: 480
<> 
<> 
<> Can anyone help on this query - I've loaded our rivers as a vector data 
<> set, but the users want to determine the length of river between two 
<> points. This is the length including all the bends, not just the 'as the 
<> crow flys' one.
<> 
<> Is there any GRASS function that can do this ?
<> 
<> Thanks for any help,
<> 
<> Steve Culshaw
<> NRA North West
<> e-mail : nra@cix.compulink.co.uk (use this for general usage)
<> or         sculshaw@cix.compulink.co.uk (Private - only checked at best 
<> biweekly)
<> 
<> 
----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: r.length
X-Sun-Content-Lines: 216

:
# *** This is r.length . *** 
#
# Idea from C.Dana Tomlin's Geographic Information Systems and 
# Cartographic Modeling; written by Rick Thompson (Center for Advanced 
#                        Spatial Technologies)
#
# @(#) r.length : creates a map of raster line lengths.
# Last revision: 6/21/94

if [ `uname -m` = mips ] ; then
 ECHON="/usr/bsd43/bin/echo -n"
elif [ `uname -s` = SunOS -a `uname -r | sed 's/\...*$//'` = 5 ] ; then
 ECHON="/usr/ucb/echo -n"
else
 ECHON="echo -n"
fi

if [ -z "$GISRC" ] ; then
 clear
 echo "This command must be run from GRASS!"
 echo ''
 exit 1
fi

if [ -f $LOCATION/cell/MASK ] ; then
 g.remove rast=MASK >> /dev/null
fi

clear
echo ''
if [ $# = 1 ] ; then
 name="$1"
else
 g.ask type=mapset prompt="Measure distance for which raster map?" element=cell desc=raster unixfile=/tmp/$$
 . /tmp/$$
 rm -f /tmp/$$
 if [ ! -f $LOCATION/cell/$name ]; then
 echo "Raster map $name not found in mapset $MAPSET."
 echo ''
 exit
 fi
fi
clear
g.region rast=$name

# Setting variables with the g.tempfile command.
tmpcat=`g.tempfile $$`
tmpcat2=`g.tempfile $$`
tmpcat3=`g.tempfile $$`
tmpcat4=`g.tempfile $$`
categ=`g.tempfile $$`
body=`g.tempfile $$`

# Selecting only categories that apply to our map with the
# r.describe -1 command. r.describe -1 will produce a single column of 
# the categories present in $name. We do not want category 0 included
# in this list, therefore, we will check for it with the variable
# $deschk. $deschk will produce the single column category list for
# $name, then keep only the first line (sed -n '1p') without any
# extra spaces (tr -d ' ')

deschk=`r.describe -1 $name | sed -n '1p' | tr -d ' '`
if [ $deschk = 0 ] ; then
 r.describe -1 $name | sed '1d' > $categ
else
 r.describe -1 $name > $categ
fi

# Producing a map of each category and calculating the length.

echo ''
echo 'Calculating lengths for each class.'

# Using the for command to individually go through each category.
for i in `cat $categ` ; do 
 r.mapcalc MASK = "if($name == $i)"
# $name must be resampled to have only the one category.
 r.resample -q in=$name out=$name.cells$i         
 g.remove rast=MASK >> /dev/null

# The next 2 r.mapcalc commands will construct 2 maps. The first
# r.mapcalc statement looks at the neighborhood around each cell
# of $name and place a value in only cells where the center cell is
# also in proximity with a cell directly over, under, or left or 
# right of it. These cells are uniquely either horizontally or
# vertically oriented to the cell being checked. These horiz./vertical
# cells are labelled with a "1". The 2nd r.mapcalc statement does
# the same thing with diagonally oriented cells and gives them a 
# value of "2".

 r.mapcalc $name.cells$i.1="if($name.cells$i && $name.cells$i[1,0],1,if($name.cells$i && $name.cells$i[0,1],1,if($name.cells$i && $name.cells$i[0,-1],1,if($name.cells$i && $name.cells$i[-1,0],1))))"

 r.mapcalc $name.cells$i.2="if($name.cells$i && $name.cells$i[-1,-1],2,if($name.cells$i && $name.cells$i[1,1],2,if($name.cells$i && $name.cells$i[1,-1],2,if($name.cells$i && $name.cells$i[-1,1],2))))"

# The next 2 variables ensure that complete fields will exist
# for all categories for each tmpcat file to be pasted together
# in $tmpcat4. Trial and error showed that no results from a
# layer could ruin the tmpcat files. $cat1chk cats the
# horiz./vertical map and tries to grab a category number
# (grep '#'), then cuts it out. The same principle applies to cat2chk.
  cat1chk=`cat $LOCATION/cats/$name.cells$i.1 | grep '#' | cut -d' ' -f2`
  cat2chk=`cat $LOCATION/cats/$name.cells$i.2 | grep '#' | cut -d' ' -f2`

# If any combination of $cat1chk and/or $cat2chk are 0, then an
# entry is sent to the tmpcat file to ensure an entry of 0 for
# the proper category is made.
# The r.stats commands count (-c) the cells quietly (q) without
# taking into account zero values (z). The first line counts the 
# total cells in the resampled map. The 2nd line counts only the 
# horiz./vertical "1"s, while the 3rd counts the diagonally 
# oriented cells. All are written to the $tmpcat files.

  if [ $cat1chk = 0 -a $cat2chk != 0 ] ; then
    r.stats -cqz $name.cells$i fs=: >> $tmpcat
    echo "1:0" >> $tmpcat2
    r.stats -cqz $name.cells$i.2 fs=: >> $tmpcat3
  elif [ $cat1chk != 0 -a $cat2chk = 0 ] ; then
    r.stats -cqz $name.cells$i fs=: >> $tmpcat
    r.stats -cqz $name.cells$i.1 fs=: >> $tmpcat2
    echo "2:0" >> $tmpcat3
  elif [ $cat1chk != 0 -a $cat2chk != 0 ] ; then
    r.stats -cqz $name.cells$i fs=: >> $tmpcat
    r.stats -cqz $name.cells$i.1 fs=: >> $tmpcat2
    r.stats -cqz $name.cells$i.2 fs=: >> $tmpcat3
  fi
done

# The $tmpcat1, etc. files are pasted together and separated by
# the : delimiter so that all the necessary information to 
# calculate the final perimeter is in one variable: $tmpcat4. The :
# delimiter though at first glance may be confusing will allow us
# to easily grab the proper field for the later awk statement.

paste -d':' $tmpcat $tmpcat2 $tmpcat3 > $tmpcat4

# Let's walk through what is being done by seeing one line of a
# sample $tmpcat4 file.
# An example is    2:125:1:116:2:49                

# 2:125   Here, 2 is the category and 125 are the total cells in the 
# resampled map that are pasted with 1:116 noting the horiz./vertical
# category 1 cells that are pasted to 2:49 noting the category 2 
# diagonal cells. 

# Producing a copy of the original $name map for the $name.length map.
# We will produce a new cats file for $name.length and substitute 
# it for the original.

g.copy rast=$name,$name.length >> /dev/null

# Setting a variable that will enable us to place the proper unit
# measurements in the final cats file. $unitchk checks for the 
# projection used in $name's cellhd file. proj:  1 means UTM meters
# and proj:   2 means feet in State Plane.

unitchk=`cat $LOCATION/cellhd/$name | grep proj | cut -d: -f2 | tr -d ' '`

# The awk command can be used to calculate the values that are needed.
# The formula says to begin with the field separator (FS) of :. Next,
# print $1 which is the category, then print a value that is the result
# of the addition of 2 different equations. First, fields $4 and $6 are 
# added together in both equations. These are the numbers of horiz./vertical
# and diagonal cells. Then subtract the total cells ($2) from them and divide 
# by 2 to find how many of the cells will be horiz./vertical or diagonal.
# Since the result may be a real number, take the integer (int) of the
# number, then subtract it from the respective horiz./vertical
# or diagonal cells to determine the actual number of cells for both 
# category 1 and category 2. The final step in the case of 
# horiz./vertical cells, is to multiply the result by the resolution.
# For diagonal cells, multiply by the square root of 2. The calculated
# categories will be sent to the $body variable.

if [ $unitchk = 1 ] ; then
  cat $tmpcat4 | awk 'BEGIN {FS = ":"} {print $1 ":" (nsres * ($4 - int((($4 + $6)- $2)/2))) + ((nsres * sqrt(2)) * ($6 - int((($4 + $6) - $2)/2))) " "units}' nsres=`cat $LOCATION/cellhd/$name | grep 'n-s' | cut -d: -f2 |tr -d ' '` units="meters" > $body
elif [ $unitchk = 2 ] ; then
  cat $tmpcat4 | awk 'BEGIN {FS = ":"} {print $1 ":" (nsres * ($4 - int((($4 + $6)- $2)/2))) + ((nsres * sqrt(2)) * ($6 - int((($4 + $6) - $2)/2))) " "units}' nsres=`cat $LOCATION/cellhd/$name | grep 'n-s' | cut -d: -f2 |tr -d ' '` units="feet" > $body

# If something other than 1 or 2 are the proj values, then use no
# units.                                           
else
  cat $tmpcat4 | awk 'BEGIN {FS = ":"} {print $1 ":" (nsres * ($4 - int((($4 + $6)- $2)/2))) + ((nsres * sqrt(2)) * ($6 - int((($4 + $6) - $2)/2)))}' nsres=`cat $LOCATION/cellhd/$name | grep 'n-s' | cut -d: -f2 |tr -d ' '` > $body
fi

# Make a header for the new raster cat file. Use g.tempfile for another
# variable named head.
head=`g.tempfile $$`

# Find the number of categories by cat-ing the $name's cats file, 
# grabbing the "#" symbol with the grep command, then cutting the 2nd
# field that is separated from the other fields by a ' ' (space).
headcat=`cat $LOCATION/cats/$name | grep '#' | cut -d' ' -f2`

# Print the following items to construct the $name.length cats file. 
echo "# $headcat categories" >> $head
echo "r.length of $name" >> $head
echo '' >> $head
echo '0.00 0.00 0.00 0.00' >> $head

# Assemble the head and body variables to make the final cats file.
cat $head $body > $LOCATION/cats/$name.length

# Removing the individual cell maps produced through the for loop.
for i in `cat $categ` ; do
 g.remove rast=$name.cells$i,$name.cells$i.1,$name.cells$i.2 >> /dev/null
done

# If the program is exited before completion, the trap command will
# remove any temp files generated by the program.
trap "rm $categ ; rm $tmpcat ; rm $tmpcat2 ; rm $tmpcat3 ; rm $tmpcat4 ; rm $head ; rm $body ; exit 0" 2

# Remove the tempfiles.
rm $categ $head $body $tmpcat $tmpcat2 $tmpcat3 $tmpcat4
echo ''
echo "$name.length is complete!"


