s.menu - creating circles
Lars Schylberg
larss at fmi.kth.se
Tue Feb 14 12:03:05 EST 1995
In an E-mail message, motte wrote:
>
> Erin O'Doherty writes:
> >
> > I need something that works like the "site characteristics report" in
> > s.menu but I want to include a radius around the site not a square.
> > I can't get away with just r.buffer because the buffers run together
> > and I'd have to sort out the sites to avoid overlap. I was hoping to
> > modify the code in s.menu to do a "circle" (I know it won't be a real
> > circle because it's a raster map) instead of a square, but I couldn't
> > figure out where it learns which cells to run stats on. I could
> > figure out a way to do it if I just had one site map to do this to,
> > but I have hundreds. What I want to do is describe the vegetation in
> > an error polygon around aerial telemetry locations. My veg map has a
> > resolution of about 3 m and the error polygon has a radius of 280 m.
> > If two sites fall in the same place I want to sample that area twice.
> > So if somebody could point me to the files in s.menu that I'd need to
> > change or figure out a slick way I could do this and repeat it for
> > many sets of sites, I'd be eternally grateful. Left to my own
> > devices I'll probably try something with r.mapcalc and r.stats.
> >
Hello !
I have an old script that you could modify a little to get it to do that.
I wrote this a couple a years ago. Ship the part about the ascii file
and read the site file directly. Then you just need to create circles instead
of the boxes that I created. A litle nice execise in geometry. Finally you
need to run r.stats to find your results.
Be happy and do some script-hacking.
Good luck
Lars
Lars Schylberg Email: lasc at celsiustech.se
CelsiusTech IT or: larss at fmi.kth.se
S-175 88 Jarfalla Tel. +46 8 580 847 03
Sweden Fax. +46 8 580 123 20
---------------- cut here ---------------------------------------------
#!/bin/sh
#
# input.arch.sh
#
# Author: Lars Schylberg
# email: larss at fmi.kth.se
#
# Date 930428
#
#------------------------------------------------------------------------
# Check if GRASS is running
#
test "$GISRC" || echo "GRASS is not running" || exit 2
#----------------------------------------------------------------------------
#
# Evaluate arguments
#
if [ $# != 2 ]
then
echo
echo Usage: `basename $0`
echo ' input=input file'
echo ' output=mapname '
echo
echo ' This program imports archelogical data from an ascii file,'
echo ' creates vector areas around the sites and '
echo ' converts these areas to raster areas. '
echo ' A site is also created with the original points.'
echo
echo ' Input file format is: '
echo ' X Y height width rotation(degres) category '
echo ' Values are separated by a space. '
echo
echo ' Cordinate system is geographic X northing, Y easting
echo ' rotation is positive clockwise
echo
exit 1
fi
# Parse input arguments
for i do
case $i in
input=*)
IN=`echo $i | sed s/input=//` ;;
in=*)
IN=`echo $i | sed s/in=//` ;;
i=*)
IN=`echo $i | sed s/i=//` ;;
output=*)
OUT=`echo $i | sed s/output=//` ;;
out=*)
OUT=`echo $i | sed s/out=//` ;;
o=*)
OUT=`echo $i | sed s/o=//` ;;
*)
echo ""
echo "Unrecognized option: $i"
echo 'Options: input=input file '
echo ' output=mapname '
exit 1
esac
done
#-------------------------------------------------------------------------
#
# Check the input arguments
#
# Test if site files is available
#
if test ! -f $IN
then
echo 'File ' $IN ' is not found'
exit 2
fi
eval `g.findfile element=dig file=$OUT`
if [ "$file" ]
then
echo "Vector file: $OUT - exists already"
exit 2
fi
eval `g.findfile element=cell file=$OUT`
if [ "$file" ]
then
echo "Raster file: $OUT - exists already"
exit 2
fi
eval `g.findfile element=site file=$OUT`
if [ "$file" ]
then
echo "Site file: $OUT - exists already"
exit 2
fi
#-----------------------------------------------------------
DIG_ASCII=$LOCATION/dig_ascii/$OUT
DIG_ATT=$LOCATION/dig_att/$OUT
SITE=$LOCATION/site_lists/$OUT
TMP=tmp.`basename $0`.$$
g.region -p | sed "s/\..*//" > $TMP
N=`grep north $TMP | sed "s/[^0-9]*//"`
S=`grep south $TMP | sed "s/[^0-9]*//"`
E=`grep east $TMP | sed "s/[^0-9]*//"`
W=`grep west $TMP | sed "s/[^0-9]*//"`
ZONE=`grep zone $TMP | sed "s/[^0-9]*//"`
/bin/rm $TMP
# Create header information
echo "ORGANIZATION: " > $DIG_ASCII
echo "DIGIT DATE: `date`" >> $DIG_ASCII
echo "DIGIT NAME: `whoami`" >> $DIG_ASCII
echo "MAP NAME: $IN" >> $DIG_ASCII
echo "MAP DATE: " >> $DIG_ASCII
echo "MAP SCALE: 50000" >> $DIG_ASCII
echo "OTHER INFO: Created by input.arch.sh" >> $DIG_ASCII
echo "ZONE: $ZONE" >> $DIG_ASCII
echo "WEST EDGE: $W" >> $DIG_ASCII
echo "EAST EDGE: $E" >> $DIG_ASCII
echo "SOUTH EDGE: $S" >> $DIG_ASCII
echo "NORTH EDGE: $N" >> $DIG_ASCII
echo "MAP THRESH: 0.00" >> $DIG_ASCII
echo "VERTI:" >> $DIG_ASCII
# Write vectors, orthogonal rectangles
# cat $IN | \
# awk '{printf "A 5\n %d %d\n %d %d\n %d %d\n %d %d\n %d %d\n", \
# $1+$3/2, $2-$4/2, $1+$3/2, $2+$4/2, \
# $1-$3/2, $2+$4/2, $1-$3/2, $2-$4/2, \
# $1+$3/2, $2-$4/2 }' >> $DIG_ASCII
# Write vectors, rectangles with rotation
# Helmert transformation
# H = height
# W = width
# X = X + H*cos(PHI) + W*sin(PHI)
# Y = Y - H*sin(PHI) + W*cos(PHI)
cat $IN | \
gawk '{ PHI=$5*3.14/180; H=$3/2; W=$4/2; \
printf "A 5\n %d %d\n %d %d\n %d %d\n %d %d\n %d %d\n", \
$1+H*(cos(PHI))-W*(-sin(PHI)), $2+H*(sin(PHI))-W*(cos(PHI)), \
$1+H*(cos(PHI))+W*(-sin(PHI)), $2+H*(sin(PHI))+W*(cos(PHI)), \
$1-H*(cos(PHI))+W*(-sin(PHI)), $2-H*(sin(PHI))+W*(cos(PHI)), \
$1-H*cos((PHI))-W*(-sin(PHI)), $2-H*(sin(PHI))-W*(cos(PHI)), \
$1+H*cos((PHI))-W*(-sin(PHI)), $2+H*(sin(PHI))-W*(cos(PHI)) \
}' >> $DIG_ASCII
cat $DIG_ASCII
# Make a vector attribute file
cat $IN | awk '{ printf" A %d %d %d\n", $2, $1, $6 }' > $DIG_ATT
# Make an site file with the same name
cat $IN | awk '{ printf"%d|%d|#%d\n", $2, $1, $6 }' > $SITE
# Make ascii to binary conversion
# of vector file
v.in.ascii input=$OUT output=$OUT
# Create topology
v.support map=$OUT option=build
# Do vector to raster conversion
v.to.rast input=$OUT output=$OUT
# Display raster data
d.rast $OUT
# Show vector file
# d.vect map=$OUT color=yellow
# Show site points
# d.sites sitefile=$OUT color=red
More information about the grass-user
mailing list