[GRASS5] ps.atlas

Mon Jul 11 16:28:34 EDT 2005

hallo,

On Mon, Jul 11, 2005 at 09:03:48PM +1200, Hamish wrote:
> Hi Jachym,
> 
> You might be able to simplify your script by having it use v.mkgrid +
> 'v.label type=centroid' + ps.map's 'labels' instruction (use CVS).
> 
> Perhaps combine the row and coln columns that v.mkgrid makes to have
> "1A, 4D, ..." style labels.
> 
> regards,
> Hamish

eeee... you may be right.. v.mkgrid could be another solution. if you
don't might, I would let the script like it is. I implemented the "1A..
" style of lables and added some (hopefully) useful columns to database.

if there are no bigger reservations to the concept, I would let is be
as it is.

Jachym

-- 
Jachym Cepicky
e-mail: jachym.cepicky at centrum.cz
URL: http://les-ejk.cz
GPG: http://les-ejk.cz/gnupg_public_key/
-------------- next part --------------
#!/bin/sh

############################################################################
#
# script:	ps.atlas
# author:	Jachym Cepicky <jachym.cepicky at centrum.cz>
# description:	Makes "atlas of maps" of region for defined number of rows 
#               and columns
# date:		2005-07-11
# copyright:	This program is free software under the GNU General Public
#		License (>=v2). 
# needs:        - sed
#               - grep
#               - awk
#               - bc
#               - ps2pdf
#############################################################################

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

#%Module
#% description: Prints atlas of maps with ps.map
#%end
#%option
#%  key: input
#%  type: string
#%  description: Mapfile for ps.map
#%  required : yes
#%End
#%option
#% key: output
#% type: string
#% description: Prefix of output files
#% required : yes
#%End
#%option
#% key: rows
#% type: integer
#% description: Number of rows in region
#% required : yes
#% answer: 1
#%End
#%option
#% key: cols
#% type: integer
#% description: Number of columns in region, maximum is 676
#% required : yes
#% answer: 1
#%End
#%option
#% key: scale
#% type: string
#% description: Scale of the output map, e.g. 1:25000
#% answer: 1:1
#% required : no
#%End
#%option
#% key: overlap
#% type: integer
#% description: Overlap in map units
#% answer: 0
#% required : no
#%End
#%option
#% key: format
#% type: string
#% description: Format of the paper (a3, a4,...) (for ps2pdfwr).
#% answer: a4
#% required : no
#%End
#%flag
#% key: p
#% description: Create pdf
#%end
#%flag
#% key: v
#% description: Create vector file, which will contain the map boundaries
#%end
#%flag
#% key: r
#% description: Rotate plot
#%End

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

#  sed
if [ -z "`which sed`" ] ; then
   echo "ERROR: Script needs sed."
   exit 1
fi

#  grep
if [ -z "`which grep`" ] ; then
   echo "ERROR: Script needs grep."
   exit 1
fi

#  awk
if [ -z "`which awk`" ] ; then
   echo "ERROR: Script needs awk."
   exit 1
fi

#  bc
if [ -z "`which bc`" ] ; then
   echo "ERROR: Script needs bc."
   exit 1
fi

#  ps2pdfwr
if [ -z "`which ps2pdfwr`" ] && [ $GIS_FLAG_p -eq 1 ]; then
   echo "ERROR: Script needs ps2pdfwr."
   exit 1
fi

# input file
if [ ! -e "$GIS_OPT_input" ]; then
    echo "ERROR: Input file $GIS_OPT_input does not exist."
    exit 1;
fi

if [ "$GIS_OPT_cols" -gt "676" ]; then
    echo "WARNING: Maximum number of columns is 676."
    echo "WARNING: You desired $input of columns."
    echo "WARNING: I will use 676 columns for this region."
    GIS_OPT_cols=676;
fi

## vars
eval `g.gisenv`
: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
input="$GIS_OPT_input"
out="$GIS_OPT_output"
rows="$GIS_OPT_rows"
cols="$GIS_OPT_cols"
scale="$GIS_OPT_scale"
rotate=""
overlap="$GIS_OPT_overlap"
oldregion="region_$$"
newvector=$out"_atlas"
TMP1="`g.tempfile pid=$$`"

# removing vector, if exists
if [ $GIS_FLAG_v -eq 1 ] && [ "`g.mlist vect |grep $newvector`" != "" ]; then
    echo "Vector exists, removing..."
    g.remove vect=$newvector
fi
# rotate flag
if [ $GIS_FLAG_r -eq 1 ]; then
    rotate="-r";
fi

# informing
echo "ps.map $rotate in=$input out=$out.eps scale=$scale"

# use only hlaf of overlaping
overlap=`echo "$overlap/2"|bc -l`;
#######################################################################
cleanup() 
{
    g.region $oldregion 
    rm -f $GISDBASE/$LOCATION_NAME/$MAPSET/windows/$oldregion
    rm -f $TMP1 
}

# what to do in case of user break:
exitprocedure()
{
    echo "User break!"
    echo "drop table $newvector"|db.execute
    cleanup
    exit 1
}
# shell check for user break (signal list: trap -l)
trap "exitprocedure" 2 3 15

#######################################################################

g.region save=$oldregion

# must be because of awk and non-english locale
export LANG=C

# setign maxeast, maxnorth, maxwest and maxsouth variables with help of g.region
eval `g.region -p|grep '^north\|^east\|^west\|^south'| sed -e s/:\ */=/ | sed -e s/^/max/`;

# be simple ps.map
if [ "$cols" == "" ]; then
    cols=1;
fi

if [ "$rows" == "" ]; then
    rows=1;
fi

########################################################################
north=$maxnorth;
west=$maxwest;

# size of steps
estep=`echo "($maxeast  - $maxwest) / $cols"|  bc -l`;
nstep=`echo "($maxnorth - $maxsouth) / $rows"| bc -l`;

# setting first region
south=`echo "($north - $nstep)" |bc -l`;
east=`echo  "($west + $estep)" |bc -l`;

# number of region
no_of_region=0;

### ascii file
if [ $GIS_FLAG_v -eq 1 ]; then
    echo "Creating ascii file + vector table"
    echo -n "ORGANIZATION: GRASS Development Team
DIGIT DATE:   
DIGIT NAME:   -
MAP NAME:     $newvector
MAP DATE:     2005
MAP SCALE:    10000
OTHER INFO:   Atlas
ZONE:         0
MAP THRESH:   0.500000
VERTI:
" > $TMP1

    ### creating database
    echo "CREATE TABLE $newvector (cat int, region int, label varchar(10), row int, col varchar(5), east double, north double, west double, south double);"| db.execute
fi

currow=0;

### rows N-S
while [ `echo $south $maxsouth |awk '{printf("%d", $1 >= $2);}'` == 1 ];
do 
    echo "north->south"    # test

    # number of current row
    currow=$(( $currow+1 ))

    curcol=0;
    # columns W-E
    while [ `echo $east $maxeast |awk '{printf("%d", $1 <= $2);}'` == 1 ];
    do
        echo "west->east"  # test

        # number of current column
        curcol=$(( $curcol+1 ))

        # charakter of current column
        if [ "$curcol" -gt "26" ]; then
            pre_alpha=`echo $curcol/26| bc`;
            curcol_alpha=`echo $pre_alpha $curcol|awk '{printf("%c%c",$1+64, $2-26+64);}'`;
        else 
            curcol_alpha=`echo $curcol|awk '{printf("%c",$1+64);}'`;
        fi

        # overlapping 10% = 5% on both sides
        over_north=`echo $north $overlap|awk '{printf("%f", $1+$2);}'`;
        over_west=`echo $west $overlap|awk '{printf("%f", $1-$2);}'`;
        over_east=`echo $east $overlap|awk '{printf("%f", $1+$2);}'`;
        over_south=`echo $south $overlap|awk '{printf("%f", $1-$2);}'`;
        east_center=`echo "$over_east+(($over_west)-($over_east))/2" |bc `
        north_center=`echo "$over_south+(($over_north)-($over_south))/2" |bc `

        g.region n=$over_north s=$over_south w=$over_west e=$over_east;

        west=$east
        east=`echo "($west + $estep)" |bc -l`;

        no_of_region=$(( $no_of_region+1 ))
        echo "-------------- REGION NUMBER $no_of_region OF $(( $rows * $cols )) ---------------"

        # storing coordinaes to ascii file
        if [ $GIS_FLAG_v -eq 1 ]; then
            echo "B 5" >> $TMP1
            echo "$over_east $over_north" >> $TMP1
            echo "$over_east $over_south" >> $TMP1
            echo "$over_west $over_south" >> $TMP1
            echo "$over_west $over_north" >> $TMP1
            echo "$over_east $over_north" >> $TMP1
            echo -e "C 1 1\n$east_center $north_center\n1 $no_of_region" >> $TMP1
            echo "INSERT INTO $newvector (cat, region, col, row, label, north, east, south, west) VALUES ($no_of_region,$no_of_region, '$curcol_alpha', $currow, '$curcol_alpha$currow', $over_north, $over_east, $over_south, $over_west);"|db.execute
        fi

        #################################################################
        # create ps & pdf
        #################################################################
        ps.map -e $rotate in=$input out=$out$no_of_region.eps scale=$scale
        if [ "$GIS_FLAG_p" = "1" ]; then
            ps2pdfwr -sPAPERSIZE=$GIS_OPT_format $out$no_of_region.eps
        fi
        #################################################################
        # 
        #################################################################
    done

    north=$south
    south=`echo "($north - $nstep)" |bc -l`;

    # go back to west!
    west=$maxwest
    east=`echo  "($west + $estep)" |bc -l`;
done;

if [ $GIS_FLAG_v -eq 1 ]; then
    # import the vector of regions
    outtmp="$newvector"_tmp
    v.in.ascii in=$TMP1 out=$outtmp format=standard
    v.clean tool=break,rmdupl in=$outtmp out=$newvector
    g.remove vect=$outtmp
    v.db.connect map=$newvector table=$newvector
fi

### cleaning
cleanup