[GRASSLIST:1837] Re: border lines

[aaime]

> Another question:
>
> I don't find a way to extract only the borderlines from rasterized
> polygons. I need this borders in a new raster map (for expample the borders
> between two different soil types).
>
> I found the script r.edge.dig. It's not exact what I need because this
> extract only between two given categories and not between all
> values in a given raster map.
>
> Another attemp was to use r.poly, but I don't find the way to make lines
> from the areas (the other way round there is a script available). After
> that I would use v.to.rast
>
>

Hi Dieter, (cc to Markus, cc to GRASSLIST)
I've looked at the problem at home, and I've prepared another version
of r.poly that has a flag to output border lines, and not area features. I've
sent it to Markus so you'll find it in the CVS next week (well, I hope so).
Now, I don't know if you can use CVS, so I've also found another solution:
here you'll find v.area2line, a script that does the opposite of v.line2area.
After running v.line2area you'll have to run v.llabel in order to attach a 
label to each line, and finally v.to.rast to get borders in raster format.
I also found a solution based on a fully raster approach: let's say that
you raster map is r_area. Then you can extract borders in the following way:

r.neighbors r_area o=r_area_diff m=diversity
r.mapcalc
r_area_border=if(r_area_diff>1,1,null())

You may also want to run r.thin on results to get one pixel wide borders.
Hope this helps. Regards
Andrea

PS: Markus, would you add v.area2line to GRASS5? Althought it's a new
command, it's only a slightly modified version of v.line2area...

-------------------------------------------------------------------------
Here's the source for v.area2line. You must put it in 
$GISBASE/scripts/v.area2line and make it executable using
chmod a+x v.area2line
-------------------------------------------------------------------------

#!/bin/sh
################################################################################# 
*** This is v.area2line . ***
#
# @(#) takes lines and translates them to areas by taking the file out to 
ASCII,# @(#) changing the 'A's to 'L's, then bringing the file back in
# @(#) dig_att file is not copied back since area labels cannot match
# @(#) lines correctly
#
# Written by Andrea Aime (aaime at libero.it),
# based on v.line2area by Rick Thompson (CAST)
#
################################################################################ 
if [ "$1" = "-help" -o "$1" = "help" ]
then
 echo v.area2line vectormap
 exit
fi
 
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
 
 
clear
echo ''
if [ $# = 1 ] ; then
  name="$1"
else
  g.ask type=old prompt="Which vector map to translate?" element=dig 
desc=vector unixfile=/tmp/$$
  . /tmp/$$
  rm -f /tmp/$$
  if [ -z $name ]; then
    echo "Vector map $name not found."
    exit
  fi
fi
 
################################################################################ 
v.out.ascii in=$name out=$name
cat $LOCATION/dig_ascii/$name | sed 's/A  /L  /' > $LOCATION/dig_ascii/$name.2
v.in.ascii in=$name.2 out=$name.2
v.support -r $name.2 op=build
rm $LOCATION/dig_ascii/$name $LOCATION/dig_ascii/$name.2
 
echo ''
$ECHON "Want to rename $name.2 to $name? (n) "
read ans
if [ "$ans" = "Y" -o "$ans" = "y" ] ; then
  g.rename vect=$name.2,$name
  echo ''
  echo "$name now contains lines."
else
  clear
  echo ''
  echo "$name.2 now contains lines."
fi