[GRASS-user] shell script to break a vector line by a vector point

Michaël Rabotin michaelrabotin at gmail.com
Mon Nov 6 09:24:03 EST 2006


Hi everybody
A few weeks ago, a discussion abour how to break a vector line by a vector
point was beginned.
I create a script in shell to solve this issue in grass 6.0.1. You need the
add on v.append and v.out.ascii.db to launch this script (add on avalaible
on http://grass.gdf-hannover.de/wiki/GRASS_AddOns)
You choose a vector line and a vector point (the points need to be on the
line to cut the lines) and the script create a new intersected vector line
Any comments or suggestions are welcome.

Cheers

Mick

 #!/bin/bash
#break a vector line by a vector point and create a new vector line
#need the Add On v.out.ascii.db and v.append
#the different vectors do not have the columns x,y,along,lcat and type
#because this script create these columns and there are no securities for
#checking their presence
#in line 179 of this script, we delete an esri column because we worked with
esri
#vector imported in grass format

if  [ -z "$GISBASE" ]
then
    echo ""
    echo "You must be in GRASS to launch this script"
    echo ""
    exit 1
fi
eval `g.gisenv`
: ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
#verification of the presence of v.out.ascii.db
if [ ! -e "/usr/lib/grass/scripts/v.out.ascii.db" ] ; then
    echo "You need the Add on v.out.ascii.db"
    echo " Segmentation process exit"
    exit
fi


#verification of the presence of v.append
if [ ! -e "/usr/lib/grass/scripts/v.append" ] ; then
    echo "You need the Add on v.append  "
    echo " Segmentation process exit"
    exit
fi

# ask the mapset where the user want to work
while test ! $Seg_mapset
     do
     echo  -n " Name of the mapset for the work : "
     read Seg_mapset
     done

while test ! -d "$GISDBASE/$LOCATION_NAME/$Seg_mapset"
do
    echo ""
    echo " the mapset $Seg_mapset doesn't exist!"
    echo ""
    echo " Choose an another name"
    read Seg_mapset
done

g.gisenv set="MAPSET=$Seg_mapset"
echo " We work now in the mapset $Seg_mapset"
echo ""


echo "enter the name of the vector line: "
read rh
#we verify the presence of this vector
g.findfile element=vector file="${rh}" > /dev/null
    if [ $? -ne 0 ] ; then
        echo "The object $rh doesn't exist !!"
        echo "Segmentation process exit"
        exit
    fi
# we verify the type (line) of this vector
v.info map=$rh | grep "lines:" | cut -c30 > /tmp/line
lines=`cat /tmp/line`
    if [ "$lines" = "0" ] ; then
        echo "The Object $rh doesn't have any lines !!"
        echo "Segmentation process exit"
        exit
    fi


echo " "
echo "The vector line choosen is $rh"
echo ""
echo "enter the name of the points vector : "
read exut
#we verify the presence of this vector
g.findfile element=vector file="${exut}" > /dev/null
    if [ $? -ne 0 ] ; then
        echo "The Object $exut doesn't exist !!"
        echo "Segmentation process exit"
        exit
    fi
# we verify the type (point) of this vector
v.info map=$exut | grep "points:" | cut -c30 > /tmp/pt
pts=`cat /tmp/pt`
    if [ "$pts" = "0" ] ; then
        echo "The Object $exut doesn't have any points!!"
        echo "Segmentation process exit"
        exit
    fi
rm /tmp/line
rm /tmp/pt

echo "The vector point choosen is $exut"
echo ""
echo "enter a name for the new vector line: "

while test ! $rhfin
     do
     echo  -n "Enter the name: "
     read rhfin
     done

while test -d "$GISDBASE/$LOCATION_NAME/$MAPSET/vector/$rhfin"
do
    echo ""
    echo " The Object $rhfin already exist"
    echo ""
    echo " Choose an another name"
    read rhfin
done



echo "the name for the new vector line is $rhfin"
echo " "
echo "***************"
echo "add the columns lcat,along,x and y in the vector points"
echo "and calcul of the x,y coordinate of these points"
echo "ALTER TABLE $exut ADD COLUMN lcat INTEGER"|db.execute
echo "ALTER TABLE $exut ADD COLUMN along DOUBLE"|db.execute
echo "ALTER TABLE $exut ADD COLUMN x DOUBLE"|db.execute
echo "ALTER TABLE $exut ADD COLUMN y DOUBLE"|db.execute

v.to.db map=$exut type=point option=coor units=me column=x,y

echo " "
echo "***************"
echo "calcul of the distance between vector line and vector point"
#allow to calculate the cat of the line where the point is
v.distance from=$exut to=$rh from_type=point to_type=line upload=cat
column=lcat
#allow to find the distance entre the beginning of the line and the point
v.distance from=$exut to=$rh from_type=point to_type=line upload=to_along
column=along

echo " "
echo "***************"
echo "export of the cat of the selected lines"

v.out.ascii.db input=$exut output=/tmp/segcat columns=lcat

cut -d"|" -f4 /tmp/segcat>/tmp/segcat2

echo " "
echo "***************"
echo "extract of the selected lines"

v.extract input=rh output=rh2 type=line file=/tmp/segcat2

echo " "
echo "***************"
echo "export of the points of the selected lines"
echo "and calcul of the distance between the beginning of the line and the
point, result in layer 2"
v.to.points input=rh2 type=line output=rh3

echo " "
echo "***************"
echo "calcul of the x,y coordinate of the points"
echo "ALTER TABLE rh3_2 ADD COLUMN x DOUBLE"|db.execute
echo "ALTER TABLE rh3_2 ADD COLUMN y DOUBLE"|db.execute
v.to.db map=rh3 type=point layer=2 option=coor units=me column=x,y

echo " "
echo "***************"
echo "creation of the column type"
echo "ALTER TABLE rh3_2 ADD COLUMN type VARCHAR(10)"|db.execute
echo "ALTER TABLE $exut ADD COLUMN type VARCHAR(10)"|db.execute
echo "UPDATE rh3_2 SET type='rh'"|db.execute
echo "UPDATE $exut SET type='exut'"|db.execute

echo " "
echo "***************"
echo "export the attributes of the points of  lines and of the points"
#segcat3: fichier contenant points et attributs des segments rh3
db.select -c table=rh3_2>/tmp/segcat3
db.select -c table=$exut>/tmp/segcat4
echo "we suppress an esri column "
cut -d"|" -f1,3,4,5,6,7 /tmp/segcat4>/tmp/segcat5


echo " "
echo "***************"
echo "sum of all the points"

sort -k2n -t'|' /tmp/segcat3 >/tmp/segcat30

sort -k2n -t'|' /tmp/segcat5 >/tmp/segcat50

sort -m -k2n -k3g -t'|' /tmp/segcat30 /tmp/segcat50>/tmp/segcat6

cut -d"|" -f2 /tmp/segcat6>/tmp/segcat7
uniq /tmp/segcat7>/tmp/segcat8


cut -d"|" -f2,3,4,5,6 /tmp/segcat6>/tmp/segcat9


echo " "
echo "***************"
echo "creation of the ascii input file for the creation of the new vector
line"

 for k in `cat /tmp/segcat8`;do
    grep "^$k" /tmp/segcat9 > /tmp/segcat10_"$k"
    grep -n exut$ /tmp/segcat10_"$k" |cut -d":" -f1 > /tmp/segcat11_"$k"
    typeset -i nb
    nb=`cat /tmp/segcat11_"$k"`
    head -"$nb" /tmp/segcat10_"$k" > /tmp/segcat20_"$k"
    echo "`cat /tmp/segcat10_"$k"`" | wc -l > /tmp/segcat13_"$k"
    typeset -i nbexut
    nbexut=`cat /tmp/segcat13_"$k"`
    typeset -i nb2
    nb2=`cat /tmp/segcat13_"$k"`-`cat /tmp/segcat11_"$k"`+1

    tail -"$nb2" /tmp/segcat10_"$k" > /tmp/segcat21_"$k"

    touch /tmp/import_"$k"
    echo "`cat /tmp/segcat20_"$k"`" | wc -l > /tmp/segcat30_"$k"
    nb20=`cat /tmp/segcat30_"$k"`
    echo "`cat /tmp/segcat21_"$k"`" | wc -l > /tmp/segcat31_"$k"
    nb21=`cat /tmp/segcat31_"$k"`

    cut -d'|' -f3,4 /tmp/segcat20_"$k" | awk -F'|' '{OFS=" ";$1=$1;print "
"$0}' > /tmp/segcat40_"$k"
    cut -d'|' -f3,4 /tmp/segcat21_"$k" | awk -F'|' '{OFS=" ";$1=$1;print "
"$0}' > /tmp/segcat41_"$k"
    echo "L  "$nb20" 1" >>/tmp/import_"$k"
    echo "`cat /tmp/segcat40_"$k"`" >>/tmp/import_"$k"
    echo " 1 1">>/tmp/import_"$k"
    echo "L  "$nb21" 1" >>/tmp/import_"$k"
    echo "`cat /tmp/segcat41_"$k"`" >>/tmp/import_"$k"
    echo " 1 2">>/tmp/import_"$k"

 done
#the file import will be the fileascii input
touch /tmp/import
    echo "ORGANIZATION:" >>/tmp/import
    echo "DIGIT DATE:" >>/tmp/import
    echo "DIGIT NAME:" >>/tmp/import
    echo "MAP NAME:" >>/tmp/import
    echo "MAP DATE:" >>/tmp/import
    echo "MAP SCALE:" >>/tmp/import
    echo "OTHER INFO:" >>/tmp/import
    echo "ZONE:" >>/tmp/import
    echo "MAP THRESH:" >>/tmp/import
    echo "VERTI:" >>/tmp/import
    for k in `cat /tmp/segcat8`;do
        echo "`cat /tmp/import_"$k"`" >>/tmp/import
    done
 rm /tmp/import_*


echo " "
echo "***************"
echo "creation of the vector line $rhfin "

v.in.ascii input=/tmp/import output="$rh"_tmp format=standard
v.category input=$rh option=print > /tmp/cattot
sort -k1n /tmp/segcat8 > /tmp/segcat8s
diff /tmp/cattot /tmp/segcat8s > /tmp/cattot2
grep -v "," /tmp/cattot2 > /tmp/cattot3
cut -d" " -f2 /tmp/cattot3 > /tmp/cattot4
v.extract input=$rh output=ext_$rh type=line file=/tmp/cattot4

rm /tmp/import
rm /tmp/cattot*



echo "`cat /tmp/segcat8s`"| awk -F'|' '{ORS=") or (cat=";$1=$1;print $0}' >
/tmp/segwhere
awk '{print NF}' /tmp/segwhere > /tmp/segwherenb
typeset -i nbwhere
nbwhere=`cat /tmp/segwherenb`-2
cut -d" " -f1-"$nbwhere" /tmp/segwhere > /tmp/segwhere2

db.copy from_driver=dbf
from_database="$GISDBASE/$LOCATION_NAME/$MAPSET/dbf/" from_table="$rh"
to_driver=dbf to_database="$GISDBASE/$LOCATION_NAME/$MAPSET/dbf/"
to_table="$rh"_tmp where="(cat=`cat /tmp/segwhere2`"
v.db.connect map="$rh"_tmp driver=dbf
database="$GISDBASE/$LOCATION_NAME/$MAPSET/dbf/" table="$rh"_tmp key=cat
layer=1

v.category input="$rh"_tmp output="$rh"_tmp2 type=line option=del
v.category input="$rh"_tmp2 output="$rh"_tmp3 type=line option=add
#utilisation du module v.append pour merger les deux rÃ(c)seaux
v.append vect1="$rh"_tmp3 vect2=ext_$rh vmerged=$rhfin vtype=line vlayer=1

rm /tmp/seg*
g.remove vect=rh3,rh2,ext_$rh,"$rh"_tmp,"$rh"_tmp2,"$rh"_tmp3
echo " "
echo "***************"
echo "Segmentation process end"
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20061106/5d23079e/attachment.html


More information about the grass-user mailing list