[GRASSLIST:6932] Creating multiple maps on subregions

[Jachym Cepicky]

hallo,
I wrote some script, which  ,,cuts`` current region to defined number of rows and columns and for each sub-region creates hardcopy map using ps.map. User can define size of region overlaping.

It would be good, if someone could test it, before I post it to GRASS Add-ons
site. 

Is there someone, how founds it useful?

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

############################################################################
#
# skript:	ps.map
# autor:	Jachym Cepicky <jachym.cepicky at centrum.cz>
# popis:	Makes "atlas" of region for defined number of rows and columns
# 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: string
#% description: Number of rows in region
#% required : yes
#% answer: 1
#%End
#%option
#% key: cols
#% type: string
#% 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


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

maxnorth="`g.region -p | grep north | sed -e s/.*:\ *//`";
maxsouth="`g.region -p | grep south | sed -e s/.*:\ *//`";
maxwest="` g.region -p | grep west | sed -e s/.*:\ *//`";
maxeast="` g.region -p | grep east | sed -e s/.*:\ *//`";


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

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

########################################################################
north=$maxnorth;
east=$maxeast;


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

south=`echo "($north - $nstep)" |bc -l`;
west=`echo  "($east - $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 $west $maxwest |awk '{printf("%d", $1 >= $2);}'` == 1 ];
    do
        echo "east->west"  # 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;
        #g.region n=$north s=$south w=$west e=$east;
        
        east=$west
        west=`echo "($east - $estep)" |bc -l`;
        
        no_of_region=$(( no_of_region=$no_of_region+1))
        echo "-------------- REGION NUMBER $no_of_region OF $(( $rows * $cols )) ---------------"

        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 east!
    east=$maxeast
    west=`echo "($east - $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