----------
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!"