r.length script

[Rick Thompson]

r.length provides raster analysis ability on raster line segments. The script 
is pretty well commented. Any feedback is appreciated. Remember, better 
accuracy is found with finer resolutions. Remove my address trailer. Rick
-------------------------------------------------------------------------------
:
# *** 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 perimeter map for each category and calculating the
# perimeter. In order to associate each perimeter with its "mother"
# polygon, each polygon will be individually masked, then examined.

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         

# Remove the MASK so that the border can be grown around the polygon.

# If the MASK is not removed, then no border can be generated.
 g.remove rast=MASK >> /dev/null

# The next 2 r.mapcalc commands will construct 2 borders. 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!"
------------------------------------------------------------------------
Rick Thompson-   Research Assistant           E-mail: rick at cast.uark.edu
Center for Advanced Spatial Technologies      Telephone: (501) 575-5736
Ozark Hall Rm. 12                             Fax: 575-3846 
University of Arkansas              
Fayetteville, AR 72701