s.menu

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.
>
>
The problem you're mentioning has been one of my concerns as well over
the last two years (basically, I wanted to obtain information about the
environment of archaeological sites within a certain radius). There is no
possibility to do this inside s.menu, unfortunately, so you should revert
to using shell-scripting. My solution was to zoom in on each site (using
the site_list as input for x,y-coordinates and then determining the
maximum extent I needed from there), and then create a 'circle' around
the site running r.los on a map called 'empty' (just zero's) with a
maximum distance of the desired radius. This produces a circular raster
map that could be used as a MASK that can be used to provide data from
r.stats. I could send you this script if you're interested, but I assume
the procedure will be clear to you. If you're really thinking of hacking
into s.menu, I'd be very much interested in any solution you might find.

Philip Verhagen
--

I thought I might add that you could produce circles of any desired distance
around points with the contributed program v.circle. The resultant vector
layer can be rasterized. This becomes the MASK that you can use as discussed
above. The advantage is that you can use a site file as input and the program
runs fairly quickly. This procedure for producing what you want is still a
multiple step procedure though. I've included the syntax for v.circle in
case you haven't seen it yet.

Usage:
v.circle [-s] [radius=value] [radius_uom=name] [area=value]
   [area_uom=name] sitefile=name output=name

Flags:
  -s Automatically run "v.support" on newly created vector file.

Parameters:
      radius Radius of circle(s) with "site_lists" point(s) as center(s). If ra
dius selected then area values are not used for computations. If both radius and
area selected, then radius has precedence over area.
               default: 0.0
  radius_uom Radius unit of measure, ie. (m)meters, ft(feet), (mi)miles.
               default: m
        area Area of circle(s) with "site_lists" point(s) as center(s). If area
selected then radius values are not used for computations.
               default: 0.0
    area_uom Area unit of measure, ie. sqm(square meters), ac(acres), sqmi(squar
e miles), hec(hectares).
               default: sqm
    sitefile GRASS site_lists file (input).
               default:
      output Vector file to be created (output).
               default:

--
-------------------------

Paul Loechl p-loechl@cecer.army.mil
USA CERL (217)352-6511x7668

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@celsiustech.se
CelsiusTech IT or: larss@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@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

All the talk of circles reminded me put out a feeler on the following:
does anyone deal with data *on* circles? In other words, does anyone
have basic data that are in the form of directions? (An example of
this type of data may include directions that animals leave a
nest/hole in search of food/mates/etc.) Just for fun-zees, I'm
contemplating writing a program to statistically test for randomness
of directions (testing the hypothesis that the n points are randomly
distributed about the circumference of the circle). I have all
the internal nasties already written (e.g., Kuiper's V and Watson's
U^2 stats) so writing a GRASS wrapper would be pretty easy.

Let me know.

--Darrell
James Darrell McCauley, PhD http://soils.ecn.purdue.edu/~mccauley/
Dept of Agricultural Engineering mccauley@ecn.purdue.edu
Purdue University tel: 317.494.9772 fax: 317.496.1115