[GRASS5] ps.atlas

[Jachym Cepicky]

Hallo developers,
after some time, I puted together script for making something like
atlas of maps using ps.map.

The script splits the current region into number of columns and rows
and in each field it runs ps.map with defined input configuratin file.

The resulting atlas of maps can be stored in vector file. Desired
overlaping can be also set.

Example result can be seen at

    http://les-ejk.cz/~jachym/img/grass/krkonose_prehledka.jpg
    [375.04KB]

It would need some testing. It works with GNU bash, version 3.00.16.

Thanks

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
#% 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


## 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!"
    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);"| db.execute
fi

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

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

        # 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=$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) VALUES ($no_of_region,$no_of_region);"|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