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