#!/bin/bash
#
# Dylan Beaudette 3/07/2008
#
#
# TODO: this module alters the region, alining to the grid two times-- this can result in a net shifting of the region extents over time. when a non-integer region resolution occurs export of the grid file will fail
#

#%Module
#% description: A simple example approach to plotting GRASS data with GMT; modify accordingly.
#%End
#%option
#% key: input
#% type: string
#% gisprompt: old,cell,raster
#% description: Name of input raster map, or basename of R,G,B triplet: .red .greed .blue
#% required : yes
#%END

#%option
#% key: output_file
#% type: string
#% description: Name of the outputfile (mapname.eps)
#% required : yes
#%END


#%option
#% key: paper_size
#% type: string
#% description: A common paper size for passing to gmtdefaults: i.e. 11x17, legal, letter, poster
#% required : yes
#%END

#%option
#% key: points_vector
#% type: string
#% gisprompt: vector
#% description: Name of a vector points file to export to simple GMT XY format
#% required : no
#%END

#%option
#% key: title
#% type: string
#% description: The title of the map
#% required : no
#%END
#%option
#% key: subtitle
#% type: string
#% description: Sub-title for the upper right hand corner of the map
#% required : no
#%END
#%option
#% key: elevation
#% type: string
#% gisprompt: old,cell,raster
#% description: Name of raster to be used for an intensity file, or contouring. Note that all maps must have the same resolution for intensity shading
#% required : no
#%END
#%option
#% key: elev_resolution
#% type: integer
#% description: Desired output resolution
#% required : no
#%END
#%option
#% key: resolution
#% type: integer
#% description: Desired output resolution
#% required : no
#%END


#%flag
#% key: c
#% description: Is raster name a basename for a triplet of R,G,B images?
#%end

#%flag
#% key: n
#% description: Do not export raster data, only re-make map.
#%end

#%flag
#% key: l
#% description: Countour the elevation data instead of an intensity file?
#%end

#%flag
#% key: p
#% description: Convert output PS|EPS files to PDF with no JPEG compression
#%end

# examples:
# ./template.sh -n -c -l -p in=m elevation=elev_4m res=1  title="Pygmy Location" points_vector=pedons@waypts
# ./template.sh  -c -l -p -n  in=m res=1 elevation=elev_4m  title="Pinnacles" points_vector=pedons@waypts  && kghostview map1.ps
# ./template.sh -c -l -n -p in=m res=1 elevation=elev_4m  title="South Trail" points_vector=pedons@waypts

# ./template.sh -n -c -p in=ikonos res=4 elevation=elev_4m  title="North Wilderness" points_vector=pedons@waypts  && kghostview map1.ps

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

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

# check if we have awk
AWK=`which awk`
if [ "$AWK" = "" ] ; then
    echo "$PROG: awk required, please install awk/gawk first" 1>&2
    exit 1
fi

# setting environment, so that awk works properly in all languages
unset LC_ALL
 LC_NUMERIC=C


#borrowed from Hamish's r.out.gmt
# what to do in case of user break:
exitprocedure()
{
    echo "User break!" 1>&2
#     cleanup
    exit 1
}
# shell check for user break (signal list: trap -l)
trap "exitprocedure" 2 3 15

# cleanup()
# {
#     if [ "$MAP_TYPE" = "DCELL" ] ; then
# 	g.remove tmp_gmt_$$ > /dev/null
#     fi
# }


output_raster="$GIS_OPT_input"
outfile="$GIS_OPT_output_file"

#setup the X,Y offset from the bottom left corner, in inches
#might need to play around with these so that none of the image is cropped by the printer
#  x_offset=".25i"     #smallest margin that my printer can use...
 x_offset="c"
 y_offset="c"


#setup the elevation contouring 
contour_interval=12	#meters
annotated_contour_interval=12	#meters
	
#setup the contour line colors
contour_color2="1/1/1"
contour_color1="230/230/230"


#setup the region based on some input:
g.region res="$GIS_OPT_resolution"


# raster based on maptype
 if [ $GIS_FLAG_c -eq 1 ] ; then #RGB triplet, just check the red map
	MAP_TYPE=`r.info -t "$output_raster".red | cut -f2 -d'='`
 else
	MAP_TYPE=`r.info -t "$output_raster" | cut -f2 -d'='`
fi



############ new approach #########

# region extents: used with -R flag
# resolution: used with the -I flag
#
# elegant approach to extracting extent and 
# region resolution as suggested by
# Bruce Raup - National Snow and Ice Data Center
#
region_assignments=`g.region -g`
eval $region_assignments

# setup flags for GMT commands
region="-R$w/$e/$s/$n"
inc="-I$ewres/$nsres"

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





# output main raster file, but only if we need to
# note the -N-9999 sets that value to NULL in the resulting GRD files. this might break some FP maps...
if [ ! $GIS_FLAG_n -eq 1 ] ; then
	case "$MAP_TYPE" in
		CELL)
		if [ $GIS_FLAG_c -eq 1 ] ; then #RGB triplet
			r.out.bin input=$output_raster.red output=- null=-9999   | xyz2grd -G$output_raster.red.grd $region $inc -ZTLh -F -N-9999
			r.out.bin input=$output_raster.green output=- null=-9999 | xyz2grd -G$output_raster.green.grd $region $inc -ZTLh -F -N-9999
			r.out.bin input=$output_raster.blue output=-  null=-9999 | xyz2grd -G$output_raster.blue.grd $region $inc -ZTLh -F -N-9999
		else #just one raster
			r.out.bin input=$output_raster output=- null=-9999 | xyz2grd -G$output_raster.grd $region $inc -ZTLh -F -N-9999
		fi ;;
		FCELL)
			r.out.bin input=$output_raster output=- null=-9999 | xyz2grd -G$output_raster.grd $region $inc -ZTLf -F -N-9999 ;; 
		DCELL)
			r.out.bin input=$output_raster output=- null=-9999 | xyz2grd -G$output_raster.grd $region $inc -ZTLd -F -N-9999 ;; 
	esac
else
	echo "Not exporting raster images..."
fi



#optionally output intensity or contouring file
if [ ! -z "$GIS_OPT_elevation" ] ; then
	
	#get the map type
	INTENSITY_MAP_TYPE=`r.info -t "$GIS_OPT_elevation" | cut -f2 -d'='`
	
	#first set the region accordingly
	#fix this!
# 	 intensity_nsres=`r.info -s "$GIS_OPT_elevation" | grep nsres | awk -F"=" '{print $2}'`
# 	 intensity_ewres=`r.info -s "$GIS_OPT_elevation" | grep ewres | awk -F"=" '{print $2}'`
	
	
	intensity_nsres=$GIS_OPT_elev_resolution
	intensity_ewres=$GIS_OPT_elev_resolution
	
	g.region res=$intensity_nsres -a
	
	
	
	# the elev file for an intensity grid
	case "$INTENSITY_MAP_TYPE" in
		CELL)
			r.out.bin input=$GIS_OPT_elevation output=- | xyz2grd -Gelevation.grd $region -I$intensity_ewres/$intensity_nsres -ZTLh -F ;;
		FCELL)
			r.out.bin input=$GIS_OPT_elevation output=- | xyz2grd -Gelevation.grd $region -I$intensity_ewres/$intensity_nsres -ZTLf -F ;; 
		DCELL)
			r.out.bin input=$GIS_OPT_elevation output=- | xyz2grd -Gelevation.grd $region -I$intensity_ewres/$intensity_nsres -ZTLd -F ;; 
	esac
		
	
	if [ $GIS_FLAG_l -eq 1 ]; then # just contours
		
		#needs some work
		#leave PS file open for one more bit of code
		more_PS_flag="-K"
		countour_interval_text=${contour_interval}"m contour interval"
	else #intensity file
		#needs some work
		#do not leave PS file open for one more bit of code
# 		more_PS_flag="-K"
		
		#convert elev grid into intensity grid
# 		note that elevation grid must be the same res as raster grid..!
		grdgradient elevation.grd -Gelevation.int -A020 -Nt0.8 
# 		grdgradient elevation.grd -Gelevation.int -A020 -Ne0.6
		
		#make a variable for the intensity parameters:
		intensity_parameters="-Ielevation.int"
	fi
		
	
fi

#optionally export a points file to XY format
if [ ! -z "$GIS_OPT_points_vector" ] ; then

	#export in a format that pstext and psxy can use
	v.out.ascii in=$GIS_OPT_points_vector |	awk -F"|" '{print $1,$2,10,0,4,"BL",$3}' > $GIS_OPT_points_vector.xy

fi


#setup some projection, and image size variables

# extract extent values for calculating map aspect ratio
extent_assignments=`g.region -ge`
eval $extent_assignments

# preserve aspect ratio from UTM E,N coordinates:
aspect_ratio=`echo $ns_extent $ew_extent | awk '{printf("%f", $1 / $2)}'`



#setup the map width in inches
#this is going to affect how much information can be placed in the margins of the paper.
#commonly a printer will enforce a .25" margin...
#on 11x17 paper, a 1:24,000 scale image will leave .21" on each side...
# map_width=10.58    #1:24,000 scale on 11x17 when using DOQQ bounds
# map_width=10.50    #1:24,184 scale on 11x17 -- compromise for ken's printer's min margins...
# map_width=15


# compute default grid ticks
# default units are meters
if [ `echo $aspect_ratio | cut -f1 -d'.'` -ge 1 ] ; then
    MAX_EXTENT=$ns_extent
    #paper should be oriented in Portrait mode
    echo Using Portrait mode
    paper_orientation_flag="-P"
    
    if [ $GIS_OPT_paper_size = "poster" ] ; then 
    	paper_size="Custom_2592x3456"
    	map_width=36
    fi
    
    if [ $GIS_OPT_paper_size = "11x17" ] ; then 
    	paper_size="11x17"
    	map_width=9
    fi
    
    if [ $GIS_OPT_paper_size = "letter" ] ; then 
    	paper_size="letter"
    	map_width=6
    fi
    
    if [ $GIS_OPT_paper_size = "legal" ] ; then 
		paper_size="legal"
    	map_width=7
    fi
    
else
    MAX_EXTENT=$ew_extent
    
    #paper should be oriented in Landscape mode (the default)
    echo "Using Landscape mode"
    paper_orientation_flag=""
    
    if [ $GIS_OPT_paper_size = "11x17" ] ; then 
    	paper_size="11x17"
    	map_width=15
    fi
    
    if [ $GIS_OPT_paper_size = "letter" ] ; then 
		paper_size="letter"	
    	map_width=8
    fi
    
    if [ $GIS_OPT_paper_size = "legal" ] ; then 
    	paper_size="legal"
    	map_width=8.5
    fi
fi


# calculate the map length based on the original aspect ratio
map_length=`echo $map_width $aspect_ratio | awk '{printf("%f", $1 * $2)}'`

# calculate the y location (in paper units) for the field sheet title
field_sheet_title_y=`echo $map_length | awk '{printf("%f", $1 + 0.33)}'`

#compile the projection string, with width/length variables
#linear projection, width inches / length inches
projection_string="X${map_width}i/${map_length}i"



# thanks Hamish for the better formatting
# compute the scale in meters
# 39.369 inch/meter
#  (round to a whole number of some sort?)
scale=`echo $ew_extent $map_width | awk '{printf("%d", ($1 * 39.369) / $2)}'`



#setup the x,y axis
#default units are meters
# 
# cruddy awk method
# 
# annotated_tic_interval=`echo $MAX_EXTENT | awk '{printf("%d", (int($1)/10) )}'`
# tic_interval=`expr $annotated_tic_interval '/' 10`

# 
# use a simple algorithm in perl, as posted at:
# http://www.perlmonks.org/?node_id=599865
# 
# 
annotated_tic_interval=`perl -w round.pl $MAX_EXTENT`
tic_interval=`expr $annotated_tic_interval '/' 10`


# debugging
echo "Annotated Tick Interval = $annotated_tic_interval"
echo "Tick Interval = $tic_interval"

# setup the comment at the bottom of the map
map_comment="1:$scale $countour_interval_text"



#set some GMT defaults here:
#note that LABEL_OFFSET and Y_AXIS_TYPE were adjusted to allow for the lack of space along horizontal margins
#we may need to adjust ANNOT_FONT_SIZE_PRIMARY to 8, so that the entire map will fit on the page
gmtset DOTS_PR_INCH 300 ANNOT_FONT_PRIMARY Times-Roman ANNOT_FONT_SIZE_PRIMARY 10 ANNOT_OFFSET_PRIMARY 0.1c HEADER_FONT Times-Roman LABEL_FONT Times-Roman LABEL_FONT_SIZE 12 LABEL_OFFSET 0.1125i HEADER_FONT_SIZE 24 PAPER_MEDIA $paper_size UNIX_TIME_POS 0i/-.5i PS_COLOR RGB Y_AXIS_TYPE ver_text COLOR_NAN 255/255/255 COLOR_BACKGROUND 255/255/255 D_FORMAT %10.0f


#A. make a basemap
	#note syntax for -B option
	psbasemap -J$projection_string $region -Bpf${tic_interval}a${annotated_tic_interval}:"Easting (m)":/f${tic_interval}a${annotated_tic_interval}:"Northing (m)"::."$GIS_OPT_title":WeSn -X$x_offset -Y$y_offset -U"$map_comment" $paper_orientation_flag -K -V > $outfile
	
	
	if [ ! -z "$GIS_OPT_subtitle" ] ; then
		#plot a sub title for the map:
		pslegend -J -R -Dx${map_width}i/${field_sheet_title_y}i/1.75i/0.5i/TR -K -O <<-EndOfLegend >>$outfile
			L 10 Times-Roman L $GIS_OPT_subtitle
		EndOfLegend
	fi	
	
	#B. make a suitible color pallette for the DOQQ-- the grid values range from 0-255
	# somehow the GMT environmental vars not accessible in GRASS session... and makecpt breaks
# 	makecpt -Cgray -T0/255/1 > pinn_work/sampling-density/field_maps/gray.cpt

# 	makecpt -Cgray -T0/255/1 -V > gray.cpt
 	OUTPUT_CPT='gray.cpt'
 
	#C. render the grid (DOQQ image) with a greyscale color pallette
	#note that we append a =2 after the grd filename- this instructs GMT to read it as an INT GRD
	#-J inherit the projection information from previous declaration
	#-R inherit the region information from previous declaration (important when regions of grd files do not match!!)
	#-K allows for further PS additions later
	#-O uses preset plotting grid
 	#grdimage /home/dylan/gis-data/PINN/field_sheet_creation/raster/doqq.grd=2 -Cdoqq.cpt -J -R -K -O -V >> $outfile
 	
 	if [ $GIS_FLAG_c -eq 1 ] ; then #RGB triplet
 		grdimage $output_raster.red.grd  $output_raster.green.grd $output_raster.blue.grd $intensity_parameters -J -R -K -O -V >> $outfile
 	else #simgle image:
		grdimage $output_raster.grd -C$OUTPUT_CPT $intensity_parameters -J -R -O -K -V >> $outfile
	fi
	
	if [ $GIS_FLAG_l -eq 1 ] ; then #make contours
# 		grdcontour elevation.grd -J -R -S4 -Q10 -C${contour_interval}  -A${annotated_contour_interval}+k$contour_color2+s6 -Wc2/$contour_color1 -Wa2/$contour_color2 -O -K >> $outfile
		
# 		only plot contours
		grdcontour elevation.grd -J -R -S4 -Q10 -C${contour_interval} -Wc2/$contour_color1 -O -K >> $outfile
	fi
	
	
	#plot exported points
	psxy $GIS_OPT_points_vector.xy  -J -R  -Sc0.2c -W1/1/1/1  -G255/1/1 -O -K >> $outfile
	
	#label them, with an outline of width 3
	# note that they are labeled with CAT
	pstext $GIS_OPT_points_vector.xy  -J -R -Dj0.25c/0.1cv255/255/255 -G1/1/1 -S3/255/255/255 -O >> $outfile
	
	
# 	#



# optionally convert EPS files to high quality PDF 
#
# note that does not always use the proper paper orientation
# 
if [ $GIS_FLAG_p -eq 1 ] ; then
	echo "Converting EPS files to PDF..."
	export GS_OPTIONS=-dPDFSETTINGS=/prepress
# 	for x in `ls *.eps`; do epstopdf $x ; done
	epstopdf $outfile && rm $outfile
fi