[GRASS-user] generating lines from points

[Paolo Craveri]

Martin,
I'm working on a script to create an area from a list of cat (here
attached), useful when we want areas from gps waypoints without direct
editing by users.

This is a very "beginner" script with not good code (work in
progress);  but handle the problem of connecting  point with category
'1' to point  with category '2' etc.


2009/1/18 Martin Landa <landa.martin at gmail.com>:
>
> If not, do you consider useful something like
>
> v.edit map tool=addline cats=1,2 layer=1
>
> or create multiple lines
>
> v.edit map tool=addline cats=1,2,3,4 layer=1 # line from 1 to 2 and
> line from 3 to 4
>

It would be great.

> Sure you can always write a simple script and use v.in.ascii approach
> (can be quite slow).

Not quite simple to me!

ciao

-- 
-- 
Paolo C.
Lat. 44° 39' 11.08'' N  Long. 7° 23' 25.26'' E
-------------- next part --------------
#!/bin/sh
############################################################################
#
# MODULE:       v.points2area.sh
# AUTHOR(S):    craveripaololivio
# PURPOSE:
# COPYRIGHT:    (C) 2008 GRASS Development Team/craveripaololivio
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#############################################################################/

#%Module
#%  description: make area from a list of selected points
#%End

#%option
#% key: input
#% type: string
#% gisprompt: old,vector,vector
#% description: Vector input map
#% key_desc: name
#% required : yes
#%end

#%option
#% key: layer
#% type: integer
#% description: Layer of input map
#% answer: 1
#% required : no
#%end

#%option
#% key: catlist
#% type: string
#% gisprompt: list of points categories
#% key_desc: list
#% description: list of points categories
#% required: yes
#%end

#%option
#% key: output
#% type: string
#% gisprompt: new,vector,vector
#% key_desc: name
#% description: Name for output vector map
#% required: yes
#%end


if [ -z "$GISBASE" ] ; then
    echo "You must be in GRASS GIS to run this program." 1>&2
    exit 1
fi

if [ "$1" != "@ARGS_PARSED@" ] ; then
    exec g.parser "$0" "$@"
fi

# # test for input vector map
eval `g.findfile element=vector file=$GIS_OPT_INPUT`
if [ ! "$file" ] ; then
     g.message -e message="Vector map <$GIS_OPT_INPUT> not found in mapset search path"
     exit 1
fi
# # does input vector has points ?
if [ "$(v.info -t $GIS_OPT_INPUT | grep points | cut -d"=" -f2)" -eq "0" ]; then
     g.message -e message="Input vector map contains no points"
     exit 1
fi
# # does layer of input vector map exists ? TODO:
# if [ -z "$(v.category  $GIS_OPT_INPUT opt=report | grep -A2 "Layer: $GIS_OPT_LAYER" | grep point)" ]; then
#      g.message -e message="Layer $GIS_OPT_LAYER not found in input vector map"
#      exit 1
# fi

# what to do in case of user break:
exitprocedure()
{
     g.message -e 'User break!'
     #delete  TMP files:
     g.remove vect=${TMPVECTPOINT},${TMPVECTAREA} > /dev/null
     rm ${TMPFILE}
     exit 1
}
# shell check for user break (signal list: trap -l)
trap "exitprocedure" 2 3 15

PROGRAM=`basename $0`
set -x
# # # setup temporary file
TMPFILE="$(g.tempfile pid=$$)"
if [ $? -ne 0 ] || [ -z "$TMPFILE" ] ; then
     g.message -e message="unable to create temporary files"
     exit 1
fi
# # # setup temporary vector
TMPVECTPOINT="vecttmpextractpoits_$$"
TMPVECTAREA="vecttmparea_$$"
# # does temporaries file exists already? (What rotten luck!)
if [ -n "$(g.mlist vect | sed -n  -e '/'"${TMPVECTAREA}"'/p' -e '/'"${TMPVECTPOINT}"'/p')" ]; then
	g.message -e message="What rotten luck!"
	exit 1
fi

v.extract -t input=${GIS_OPT_INPUT} output=${TMPVECTPOINT} type=point \
	layer=${GIS_OPT_LAYER} list=${GIS_OPT_CATLIST} new=-1

if [ ${GIS_OPT_LAYER} -ne "1" ]; then
	table="${TMPVECTPOINT}_${GIS_OPT_LAYER}"
# # # 	se non creo la tabella sul layer 1 ottengo un errore!!
	v.db.addtable map=${TMPVECTPOINT} layer=1
else
	table="${TMPVECTPOINT}"
fi
echo "CREATE TABLE $table (cat integer, x double precision, y double precision)" | db.execute
# v.db.addtable "${TMPVECTPOINT}"  col="x double precision, y double precision"
v.db.connect map=${TMPVECTPOINT} layer=${GIS_OPT_LAYER} table=$table key=cat
v.to.db "${TMPVECTPOINT}" layer=${GIS_OPT_LAYER} option=cat
v.to.db "${TMPVECTPOINT}" layer=${GIS_OPT_LAYER} option=coor col=x,y
# # # standard mode. First line
echo "VERTI:" > ${TMPFILE}
# # # get the first category from list
firstcat=`expr match "${GIS_OPT_CATLIST}" '\([0-9]*\)'`
# # # first and last list element are the same
# # # so boundaries can always close correctly
GIS_OPT_CATLIST="${GIS_OPT_CATLIST},${firstcat}"
# # # # # populate ascii temp file from points coordinates
# # # # # TODO: better way to do that?
# # # area vertex counter
count=0
for i in $(echo $GIS_OPT_CATLIST | sed -e 's/,/ /g'); do
	count=$(($count+1))
	v.db.select -c "${TMPVECTPOINT}" layer=${GIS_OPT_LAYER} col=x,y where="cat=$i" fs=" " >> ${TMPFILE}
done
# # # # # create area TODO: find somthing  better than sed -e '2i\B '$count'\'
cat ${TMPFILE}
cat ${TMPFILE} | sed -e '2i\B '"$count"'\' | v.in.ascii out="${TMPVECTAREA}" format=standard --overwrite

# # # # # add centroid to closed boundary
v.centroids in=${TMPVECTAREA} out=${GIS_OPT_OUTPUT}

# # # clenup
g.remove vect=${TMPVECTPOINT},${TMPVECTAREA} > /dev/null
rm ${TMPFILE}

if [ -z "$GRASS_VERBOSE" ] || [ "$GRASS_VERBOSE" -gt 0 ] ; then
   g.message message="Area ${GIS_OPT_OUTPUT} created from point of category ${GIS_OPT_CATLIST}"
fi