[GRASS-SVN] r43645 - in grass-addons/raster: . r.csr
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Sep 23 06:06:38 EDT 2010
Author: neteler
Date: 2010-09-23 10:06:38 +0000 (Thu, 23 Sep 2010)
New Revision: 43645
Added:
grass-addons/raster/r.csr/
grass-addons/raster/r.csr/r.csr
Log:
Eric Patton: colored, shaded-relief raster
Added: grass-addons/raster/r.csr/r.csr
===================================================================
--- grass-addons/raster/r.csr/r.csr (rev 0)
+++ grass-addons/raster/r.csr/r.csr 2010-09-23 10:06:38 UTC (rev 43645)
@@ -0,0 +1,505 @@
+#! /bin/bash
+#
+############################################################################
+#
+# MODULE: r.csr for Grass 6.*/7.*
+#
+# AUTHOR: Eric Patton, Geological Survey of Canada (Atlantic)
+# <epatton at nrcan dot gc dot ca>
+#
+# PURPOSE: To allow the batch creation of coloured, shaded-relief
+# Grass rasters, optionally with tiff export. Uses
+# (modifiable) shading defaults of 45alt, 315az, 20x.
+#
+# COPYRIGHT: (C) 2006-2010 by Eric Patton
+#
+# This program is free software under the GNU General Public
+# License (>=v3). Read the file COPYING that comes with GRASS
+# for details.
+#
+# Last Modified: March 16, 2010
+# by Author
+#############################################################################
+#
+#%Module
+#% description: Single or batch creation of coloured, shaded-relief rasters, with optional tiff export.
+#%END
+
+#%flag
+#% key: a
+#% description: Export ESRI Arc ascii grid of input DEM
+#%END
+
+#%flag
+#% key: m
+#% description: Median filter raster using 3x3 window.
+#%END
+
+#%flag
+#% key: b
+#% description: Bypass shaded-map creation (still expects a shade map named MAP_shade).
+#%END
+
+#%flag
+#% key: r
+#% description: Export region is set by g.region, not input raster extents.
+#%END
+
+#%flag
+#% key: s
+#% description: Export tiff image and worldfile of shaded-relief raster.
+#%END
+
+#%flag
+#% key: t
+#% description: Export tiff image and worldfile of coloured, shaded-relief raster.
+#%END
+
+#%option
+#% key: map
+#% type: string
+#% gisprompt: old,cell,raster
+#% required: yes
+#% description: Input raster filename or wildcard seach pattern
+#%END
+
+#%option
+#% key: resolution
+#% type: integer
+#% required: no
+#% description: Output resolution of colored, shaded relief raster
+#%END
+
+#%option
+#% key: shademap
+#% type:string
+#% required: no
+#% gisprompt: old,cell,raster
+#% description: Name of raster map to use for shaded-relief (overrides -b flag)
+#%END
+
+#%option
+#% key: passes
+#% type: integer
+#% required: no
+#% answer: 1
+#% description: With -m flag enabled, number of times to repeat median filter
+#%END
+
+#%option
+#% key: altitude
+#% type: integer
+#% required : no
+#% answer: 45
+#% description: Altitude of the sun in degrees above the horizon (must be 1-89)
+#%END
+
+#%option
+#% key: azimuth
+#% type: integer
+#% description: Azimuth of the sun in degrees to the east of north (must be 0-360)
+#% required : no
+#% answer: 315
+#%END
+
+#%option
+#% key: zmult
+#% type: double
+#% description: Factor for exaggerating relief (default=10)
+#% required : no
+#% answer: 10
+#%END
+
+#%option
+#% key: units
+#% type: string
+#% description: Set scaling factor (applies to lat./long. locations only, none: scale=1)
+#% required : no
+#% options: none,meters,feet
+#% answer: none
+#%END
+
+#%option
+#% key: colormap
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of another raster in current mapset to copy color table from
+#% required: no
+#%END
+
+#%option
+#% key: rules
+#% type: string
+#% description: Name of rules file in current directory containing rules for color table
+#% required: no
+#%END
+
+if [ -z "$GISBASE" ] ; then
+ echo "You must be in GRASS GIS to run this program."
+ exit 1
+fi
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+
+SCRIPT=`basename $0`
+
+# Setup clean exit for Ctrl-C or similar breaks.
+trap 'echo -e "\n\nUser break or similar caught; Exiting.\n" ; exit 1' 2 3 15
+
+# Capture environment variables.
+eval `g.gisenv`
+: ${GISDBASE?} ${GISBASE?} ${LOCATION_NAME?} ${MAPSET?}
+
+ALTITUDE=$GIS_OPT_altitude
+AZIMUTH=$GIS_OPT_azimuth
+ZMULT=$GIS_OPT_zmult
+UNITS=$GIS_OPT_units
+PATTERN=$GIS_OPT_map
+RULES=$GIS_OPT_rules
+COLOR_RAST=$GIS_OPT_colormap
+SHADEMAP=$GIS_OPT_shademap
+PASSES=$GIS_OPT_passes
+COLOR_DIR=colr
+GRASS_VERSION=`g.version | awk '{print $2}'`
+RESOLUTION=$GIS_OPT_resolution
+
+# Not sure if there's a better way to check if any rasters
+# match the search pattern, but here goes:
+MATCHES=`g.mlist pattern="$PATTERN" | wc -l`
+
+# Abort if no matches are found.
+if [ "$MATCHES" -eq 0 ] ; then
+ echo "$SCRIPT: No rasters matched the search pattern!"
+ exit 1
+fi
+
+# Now iterate through each match of PATTERN and colour raster,shade, and combine.
+for MAP in `g.mlist pattern="$PATTERN"` ; do
+
+ ORIGINAL_MAP=$MAP
+ # Because we'll need the unmodified input file name in the r.color section
+
+ # Use specified resolution if given, otherwise set resolution to that of input raster.
+ if [ -z "$RESOLUTION" ] ; then
+
+ eval `r.info -s $MAP`
+ : ${ewres?} ${nsres?}
+
+ RESOLUTION=$ewres
+
+ g.region res=$RESOLUTION
+ else
+
+ g.region -a res=$RESOLUTION
+ fi
+
+ # Check for the presence of -r flag to determine if the export region should determined by g.region
+ # or the maximum extent of each input raster. The resolution set by
+ # RESOLUTION parameter is still respected, if -r flag is set.
+ # If the -r flag is NOT passed, the default behavior will be
+ # to set the region to that of each raster in turn.
+
+ if [ "$GIS_FLAG_r" -eq 0 ] ; then
+
+ g.region -a rast=$MAP res=$RESOLUTION
+ fi
+
+
+ echo -e "\n\n\n======================================"
+ echo -e "\n$SCRIPT"
+ echo -e "\nProcessing input map $MAP..."
+ echo -e "\n\nFormatting check begun..."
+
+ # Check filename of raster to make sure no period exists in filename.
+ # r.composite doesn't accept periods in the input names.
+
+ NAME_CHECK=`echo "$MAP" | grep "\."`
+
+ if [ -n "$NAME_CHECK" ] ; then
+
+ echo -e "\nA period (".") exists in the input raster filename."
+ echo -e "Replacing period with an underscore so r.composite will execute.\n"
+
+ NEW_MAP=`echo "$MAP" | tr "." "_"`
+
+ # The syntax of the overwrite flags varies slightly from different versions of Grass.
+ if [ "$GRASS_VERSION" != "6.5.svn" ] ; then
+
+ g.rename rast=${MAP},${NEW_MAP} 2>&1 > /dev/null
+
+ else
+ g.rename --o --q rast=${MAP},${NEW_MAP} 2>&1 > /dev/null
+
+ fi
+
+ if [ $? -eq 0 ] ; then
+ echo -e "\nSuccessfully renamed $MAP to $NEW_MAP..."
+
+ MAP=$NEW_MAP
+ ORIGINAL_MAP=$NEW_MAP
+ else
+ echo -e "\n$SCRIPT: Error: Unable to rename $MAP to $NEW_MAP. Exiting."
+ exit 1
+ fi
+
+ fi
+
+
+
+ ###########################################
+ # PROCEDURES FOR MEDIAN FILTERING (-m flag)
+ ###########################################
+
+
+ if [ "$GIS_FLAG_m" -eq 1 ] ; then
+
+ echo -e "\n\n--------------------------------------"
+ echo -e "\nFilling Null Values"
+ echo -e "\nFilling holes with 3x3 median filter in ${MAP}..."
+
+ # Procedure for filling nulls with multiple median filter passes
+ if [ "$PASSES" -gt 1 ] ; then
+
+ COUNT=1
+ while [ "$COUNT" -le "$PASSES" ] ; do
+
+ echo ""
+ r.neighbors input=${MAP} output=fill.tmp method=median size=3
+
+ r.mapcalc "${MAP}_fill = if(isnull(${MAP}), fill.tmp, ${MAP})"
+
+ MAP=${MAP}_fill
+
+ echo -e "Finished median filter pass $COUNT..."
+
+ COUNT=$(($COUNT + 1))
+ g.remove rast=fill.tmp 1>&2 > /dev/null
+
+ done
+
+ else
+ echo ""
+ r.neighbors input=${MAP} output=fill.tmp method=median size=3
+ r.mapcalc "${MAP}_fill = if(isnull(${MAP}), fill.tmp, ${MAP})"
+ MAP=${MAP}_fill
+ g.remove rast=fill.tmp 1>&2 > /dev/null
+
+ fi
+
+ if [ $? -eq 0 ] ; then
+ echo "Finished median filtering."
+
+ else
+ echo -e "\n$SCRIPT: Error: r.neighbors failed during median filtering. Exiting."
+ exit 1
+ fi
+
+ fi
+
+ RED=${MAP}_r
+ GREEN=${MAP}_g
+ BLUE=${MAP}_b
+ SHADE=${MAP}_shade
+
+
+ ##################################
+ # PROCEDURES FOR SHADING (-b flag)
+ ##################################
+
+ # If $SHADEMAP is given, then use it as the shaded-relief map and combine it with
+ # $MAP. If $SHADEMAP is not given and -b is given, bypass creating a shaded-relief
+ # from the input raster, and use a map named $MAP_shade as the shaded-relief map.
+ # If neither $SHADEMAP nor the -b flag are given, generate a shaded-relief map from scratch
+ # using the input raster.
+
+
+ echo -e "\n\n--------------------------------------"
+
+ # If -b is not given, use r.shaded.relief to create a shade map.
+ if [ "$GIS_FLAG_b" -eq 0 ] ; then
+ echo -e "\nr.shaded.relief\n"
+
+ if [ "$GRASS_VERSION" != "7.0.svn" ] ; then
+ r.shaded.relief --o map=${MAP} shadedmap=${SHADE} altitude=${ALTITUDE} azimuth=${AZIMUTH} zmult=${ZMULT} units=${UNITS}
+ else
+ r.shaded.relief --o input=${MAP}
+ oBonavista_AllBathy_Feb16_2009_10m_fill_shadeutput=${SHADE} altitude=${ALTITUDE} azimuth=${AZIMUTH} zmult=${ZMULT} units=${UNITS}
+ fi
+
+ if [ $? -ne 0 ] ; then
+ echo -e "\n$SCRIPT: r.shaded.relief was unsuccessful."
+
+ # Cleanup if r.shaded.relief failed
+ eval `g.findfile element=cell mapset=$MAPSET file=${SHADE}`
+
+ if [ -n "$file" ] ; then
+ g.remove rast=${SHADE} > /dev/null
+ exit 1
+ else
+ exit 1
+ fi
+ fi
+
+ else # Use the $SHADEMAP if it was given.
+
+ if [ -n "$SHADEMAP" ] ; then
+ SHADE=${SHADEMAP}
+
+ else
+ # Check if $MAP_shade already exists, since the bypass shade flag (-b) was passed,
+ # and no $SHADEMAP parameter was given.
+
+ eval `g.findfile element=cell mapset=$MAPSET file=${SHADE}`
+ SHADE_CHECK=$file
+
+ echo -e "\nThe shaded-relief option has been skipped (-b flag). Note that"
+ echo -e "r.his still expects a shaded-relief file named ${SHADE} to"
+ echo -e "combine with the colour map ${MAP}."
+ echo -e "\nChecking if shade map already exists..."
+
+ if [ -n "$SHADE_CHECK" ] ; then
+ echo -e "Shade map $SHADE found."
+
+ else
+ echo -e "\n$SCRIPT: Error: No shade map found!"
+ echo "r.csr -b expects an existing shade map named ${SHADE} in the current MAPSET."
+ echo -e "Disable -b flag to continue."
+ exit 1
+ fi
+ fi
+ fi
+
+
+ echo -e "\n\n--------------------------------------"
+ echo -e "\nCombining colours with shade - r.his and r.composite\n"
+ echo -e "\nCreating new csr maps...\n"
+
+
+ ####################################################
+ # PROCEDURES FOR ASSIGNING COLOR TABLES TO INPUT MAP
+ ####################################################
+
+ # Check if color rules file or existing raster colormap is specified.
+
+ if [ -z "$RULES" -a -z "$COLOR_RAST" ] ; then
+
+ # Use the original input raster color table if it exists.
+ if [ -f ${GISDBASE}/${LOCATION_NAME}/${MAPSET}/${COLOR_DIR}/${ORIGINAL_MAP} ] ; then
+
+ r.colors map=${MAP} rast=${ORIGINAL_MAP}
+
+ # If the raster map was renamed from $ORIGINAL_MAP to $MAP, and a color
+ # table exists for $MAP, use it.
+ elif [ -f ${GISDBASE}/${LOCATION_NAME}/${MAPSET}/${COLOR_DIR}/${MAP} ] ; then
+
+ r.colors map=${MAP} rast=${MAP}
+
+ else
+ echo "Pre-existing color map for ${MAP} not found"
+ echo "and neither rules nor colormap parameters given."
+ echo -e "Defaulting to rainbow colormap...\n"
+
+ r.colors map=${MAP} color=bcyr
+
+ fi
+
+ # Copy the color table of another raster if COLOR_RAST parameter is given.
+ elif [ -n "$COLOR_RAST" ] ; then
+
+ r.colors map=${MAP} rast=${COLOR_RAST}
+
+ elif [ -n "$RULES" ] ; then
+
+ #Use the color table specified in a rules text file in the current directory.
+ cat $RULES | r.colors map=${MAP} color=rules
+ fi
+
+ ####################################
+ # COMBINATION OF COLOR AND SHADE MAPS
+ ####################################
+
+ # Color map and shade map get combined in this step.
+
+ r.his -n --o h_map=$MAP i_map=$SHADE r_map=$RED g_map=$GREEN b_map=$BLUE
+
+ # Confirmation that r.his has actually created 3 rgb maps.
+ eval `g.findfile element=cell mapset=$MAPSET file=${RED}`
+ RED_CHECK=$file
+
+ eval `g.findfile element=cell mapset=$MAPSET file=${GREEN}`
+ GREEN_CHECK=$file
+
+ eval `g.findfile element=cell mapset=$MAPSET file=${BLUE}`
+ BLUE_CHECK=$file
+
+ if [ -z "$RED_CHECK" -o -z "$GREEN_CHECK" -o -z "$BLUE_CHECK" ] ; then
+ echo -e "\n$SCRIPT: r.his failed to create one or more RGB maps. Exiting."
+ echo -e "\nCleaning any RGB maps..."
+
+ g.remove rast=$RED,$GREEN,$BLUE > /dev/null
+ exit 1
+ fi
+
+ COMBINED=${MAP}_shade_comb
+
+ # Merging RGB maps.
+
+ r.composite -d --o red=${RED} green=${GREEN} blue=${BLUE} output=${COMBINED}
+
+ if [ $? -ne 0 ] ; then
+ echo -e "\n$SCRIPT: Error occurred during r.composite export. Exiting."
+ exit 1
+ fi
+
+ ##############################################
+ # EXPORTING OUTPUT - TIFFS AND ARC ASCII GRIDS
+ ##############################################
+
+ # Procedures for exporting tiffs
+
+ echo -e "\n--------------------------------------"
+ echo -e "\nExport Section\n"
+
+ # Procedure for exporting Arc ascii grid for -a flag.
+ if [ "$GIS_FLAG_a" -eq 1 ] ; then
+ echo -n -e "\nCreating Arc ascii grid...standby."
+
+ r.out.arc input=${MAP} output=${MAP}.asc dp=3
+ echo -e "\nCreated Arc ascii grid.\n"
+
+ fi
+
+ if [ "$GIS_FLAG_t" -eq 1 ] ; then
+ r.out.tiff -t input=${COMBINED} output=${COMBINED}
+ gdal_translate -a_srs "`g.proj -wf`" ${COMBINED}.tif tmp.tif && mv tmp.tif ${COMBINED}.tif
+
+ if [ "$?" -eq 0 ] ; then
+ rm ${COMBINED}.tfw
+ echo -e "\nCreated tiff ${COMBINED}.tif."
+ fi
+
+ # Check to see if the -s flag was passed to export a tiff image of the shaded-relief map too.
+ if [ "$GIS_FLAG_s" -eq 1 ] ; then
+ r.out.tiff -t input=${SHADE} output=${SHADE}
+ gdal_translate -a_srs "`g.proj -wf`" ${SHADE}.tif tmp.tif && mv tmp.tif ${SHADE}.tif
+
+ if [ "$?" -eq 0 ] ; then
+ rm ${SHADE}.tfw
+ echo -e "\nCreated tiff ${SHADE}.tif."
+ fi
+ fi
+ fi
+
+ # CLEANUP
+
+ echo -e "\n\n--------------------------------------"
+ echo -e "\nCleaning up rgb files...\n"
+
+ g.remove rast=$RED,$GREEN,$BLUE > /dev/null
+ echo "Done."
+
+done
+
+exit 0
Property changes on: grass-addons/raster/r.csr/r.csr
___________________________________________________________________
Added: svn:executable
+ *
More information about the grass-commit
mailing list