|
@@ -1,437 +0,0 @@
|
|
|
-#!/bin/sh
|
|
|
-
|
|
|
-############################################################################
|
|
|
-#
|
|
|
-# MODULE: r.in.gdalwarp for GRASS 6
|
|
|
-# AUTHOR(S): Cedric Shock
|
|
|
-# PURPOSE: To warp and import data
|
|
|
-# COPYRIGHT: (C) 2006 by Cedric Shock
|
|
|
-#
|
|
|
-# This program is free software under the GNU General Public
|
|
|
-# License (>=v2). Read the file COPYING that comes with GRASS
|
|
|
-# for details.
|
|
|
-#
|
|
|
-#############################################################################
|
|
|
-
|
|
|
-# Requires:
|
|
|
-# gdalwarp
|
|
|
-
|
|
|
-#%Module
|
|
|
-#% description: Warps and imports GDAL supported raster file complete with correct NULL values
|
|
|
-#% keywords: raster, rotate, reproject
|
|
|
-#%End
|
|
|
-#%flag
|
|
|
-#% key: p
|
|
|
-#% description: Don't reproject the data, just patch it.
|
|
|
-#%end
|
|
|
-#%flag
|
|
|
-#% key: e
|
|
|
-#% description: Extend location extents based on new dataset
|
|
|
-#%end
|
|
|
-#%flag
|
|
|
-#% key: c
|
|
|
-#% description: Make color composite image if possible
|
|
|
-#%end
|
|
|
-#%flag
|
|
|
-#% key: k
|
|
|
-#% description: Keep band numbers instead of using band color names
|
|
|
-#%end
|
|
|
-#%option
|
|
|
-#% key: input
|
|
|
-#% type: string
|
|
|
-#% description: Raster file or files to be imported. If multiple files are specified they will be patched together.
|
|
|
-#% multiple: yes
|
|
|
-#% required : yes
|
|
|
-#%end
|
|
|
-#%option
|
|
|
-#% key: output
|
|
|
-#% type: string
|
|
|
-#% description: Name for resultant raster map. Each band will be name output.bandname
|
|
|
-#% required : yes
|
|
|
-#%end
|
|
|
-#%option
|
|
|
-#% key: s_srs
|
|
|
-#% type: string
|
|
|
-#% description: Source projection in gdalwarp format
|
|
|
-#% required : yes
|
|
|
-#%end
|
|
|
-#%option
|
|
|
-#% key: method
|
|
|
-#% type: string
|
|
|
-#% description: Reprojection method to use
|
|
|
-#% options:nearest,bilinear,cubic,cubicspline
|
|
|
-#% answer:nearest
|
|
|
-#% required: yes
|
|
|
-#%end
|
|
|
-#%option
|
|
|
-#% key: warpoptions
|
|
|
-#% type: string
|
|
|
-#% description: Additional options for gdalwarp
|
|
|
-#% required : no
|
|
|
-#%end
|
|
|
-
|
|
|
-
|
|
|
-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
|
|
|
-
|
|
|
-g.message -d "[r.in.gdalwarp]"
|
|
|
-
|
|
|
-#### setup temporary file
|
|
|
-TMP="`g.tempfile pid=$$`"
|
|
|
-if [ $? -ne 0 ] || [ -z "$TMP" ] ; then
|
|
|
- g.message -e "Unable to create temporary files"
|
|
|
- exit 1
|
|
|
-fi
|
|
|
-
|
|
|
-#######################################################################
|
|
|
-# name: exitprocedure
|
|
|
-# purpose: removes all temporary files
|
|
|
-#
|
|
|
-exitprocedure()
|
|
|
-{
|
|
|
- g.message -e 'User break!'
|
|
|
- rm -f "${TMP}"*
|
|
|
- exit 1
|
|
|
-}
|
|
|
-trap "exitprocedure" 2 3 15
|
|
|
-
|
|
|
-# check if we have gdalwarp
|
|
|
-if [ $GIS_FLAG_P -ne 1 ] ; then
|
|
|
- if [ ! -x "`which gdalwarp`" ] ; then
|
|
|
- g.message -e "gdalwarp is required, please install it first"
|
|
|
- exit 1
|
|
|
- fi
|
|
|
-fi
|
|
|
-
|
|
|
-defaultIFS="$IFS"
|
|
|
-
|
|
|
-# warpimport file map
|
|
|
-# Warps and imports a single file as a named map
|
|
|
-# This is NOT a function
|
|
|
-# It is a subroutine that makes extensive use of globals
|
|
|
-warpimport () {
|
|
|
- FILE="$1"
|
|
|
- MAP=$2
|
|
|
- WARPFILE="${TMP}warped.geotiff"
|
|
|
- TMPMAPNAME=${2}_tmp
|
|
|
-
|
|
|
- # Warp it and convert it to geotiff with alpha band to be sure we get nulls:
|
|
|
- g.message -d message="gdalwarp -s_srs \"$GIS_OPT_S_SRS\" -t_srs \"`g.proj -wf`\" \"$FILE\" \"$WARPFILE\" $GIS_OPT_WARPOPTIONS $METHOD"
|
|
|
-
|
|
|
- gdalwarp -s_srs "$GIS_OPT_S_SRS" -t_srs "`g.proj -wf`" \
|
|
|
- "$FILE" "$WARPFILE" $GIS_OPT_WARPOPTIONS $METHOD
|
|
|
-
|
|
|
- if [ $? -ne 0 ] ; then
|
|
|
- g.message -e message="`basename $0`: gdalwarp failure."
|
|
|
- rm -f "${TMP}"*
|
|
|
- exit 1
|
|
|
- fi
|
|
|
-
|
|
|
- #Import it into a temporary map:
|
|
|
- r.in.gdal $FLAGS input="$WARPFILE" output="$TMPMAPNAME" --quiet
|
|
|
- if [ $? -ne 0 ] ; then
|
|
|
- g.message -e message="`basename $0`: r.in.gdal failure."
|
|
|
- rm -f "${TMP}"*
|
|
|
- exit 1
|
|
|
- fi
|
|
|
-
|
|
|
- # Remove temporary file
|
|
|
- rm -f "$WARPFILE"
|
|
|
-
|
|
|
- #Get a list of channels:
|
|
|
- PATTERN="$TMPMAPNAME*"
|
|
|
- g.message -d message="Pattern: [$PATTERN]"
|
|
|
- CHANNEL_LIST=`g.mlist type=rast pattern=$PATTERN`
|
|
|
- g.message -d message="Channel list: [$CHANNEL_LIST]"
|
|
|
- # Get a list of suffixes:
|
|
|
- CHANNEL_SUFFIXES=`echo "$CHANNEL_LIST" | sed "s/^${TMPMAPNAME}/sfx=/"`
|
|
|
- # test for single band data
|
|
|
- if [ "$CHANNEL_SUFFIXES" = 'sfx=' ] ; then
|
|
|
- CHANNEL_SUFFIXES=""
|
|
|
- fi
|
|
|
- g.message -d message="Channel suffixes: [$CHANNEL_SUFFIXES]"
|
|
|
-
|
|
|
- # Add to the list of all suffixes:
|
|
|
- SUFFIXES=`echo "$SUFFIXES
|
|
|
-$CHANNEL_SUFFIXES" | sort -u`
|
|
|
-
|
|
|
- IFS=$defaultIFS
|
|
|
- for SUFFIX in $CHANNEL_SUFFIXES ; do
|
|
|
- eval "$SUFFIX"
|
|
|
- LAST_SUFFIX=$sfx
|
|
|
- done
|
|
|
-
|
|
|
- # Find the alpha layer
|
|
|
- if [ $GIS_FLAG_K -eq 1 ] ; then
|
|
|
- ALPHALAYER=${TMPMAPNAME}${LAST_SUFFIX}
|
|
|
- else
|
|
|
- ALPHALAYER=${TMPMAPNAME}.alpha
|
|
|
- fi
|
|
|
- # test to see if the alpha map exists
|
|
|
- g.findfile element=cell file="$ALPHALAYER" > /dev/null
|
|
|
- if [ $? -ne 0 ] ; then
|
|
|
- ALPHALAYER=""
|
|
|
- fi
|
|
|
-
|
|
|
- # Calculate the new maps:
|
|
|
- for SUFFIX in $CHANNEL_SUFFIXES ; do
|
|
|
- eval "$SUFFIX"
|
|
|
-
|
|
|
- # if no suffix, processing is simple (e.g. elevation has only 1 band)
|
|
|
- # --never reached?
|
|
|
- if [ -z "$sfx" ] ; then
|
|
|
- # run r.mapcalc to crop to region
|
|
|
- GRASS_VERBOSE=1 \
|
|
|
- r.mapcalc "${MAP}${sfx} = ${TMPMAPNAME}${sfx}"
|
|
|
- continue
|
|
|
- fi
|
|
|
-
|
|
|
- g.message -d message="alpha=[$ALPHALAYER] MAPsfx=[${MAP}${sfx}] tmpname=[${TMPMAPNAME}${sfx}]"
|
|
|
- if [ -n "$ALPHALAYER" ] ; then
|
|
|
- # Use alpha channel for nulls:
|
|
|
-# problem: I've seen a map where alpha was 1-255; 1 being transparent. what to do?
|
|
|
-# (Geosci Australia Gold layer, format=tiff)
|
|
|
- GRASS_VERBOSE=1 \
|
|
|
- r.mapcalc "${MAP}${sfx} = if( $ALPHALAYER, ${TMPMAPNAME}${sfx}, null() )"
|
|
|
- else
|
|
|
- g.copy rast="${TMPMAPNAME}${sfx}","${MAP}${sfx}" --quiet
|
|
|
- fi
|
|
|
-
|
|
|
- # Copy the color tables:
|
|
|
- r.colors map=${MAP}${sfx} rast=${TMPMAPNAME}${sfx} --quiet
|
|
|
-
|
|
|
- # Make patch lists:
|
|
|
- sfx2=`echo $sfx | sed "s/\./_/"`
|
|
|
- # This is a hack to make the patch lists empty:
|
|
|
- if [ $TITER -eq 0 ] ; then
|
|
|
- eval "PATCHES_$sfx2=\"\""
|
|
|
- fi
|
|
|
- eval "PATCHES_$sfx2=\${PATCHES_$sfx2},${MAP}${sfx}"
|
|
|
- done
|
|
|
-
|
|
|
- # if no suffix, processing is simple (e.g. elevation has only 1 band)
|
|
|
- if [ -z "$sfx" ] ; then
|
|
|
- # run r.mapcalc to crop to region
|
|
|
- GRASS_VERBOSE=1 \
|
|
|
- r.mapcalc "${MAP}${sfx} = ${TMPMAPNAME}${sfx}"
|
|
|
- r.colors map=${MAP}${sfx} rast=${TMPMAPNAME}${sfx} --quiet
|
|
|
- fi
|
|
|
-
|
|
|
- # Remove the old channels:
|
|
|
- CHANNEL_LIST_COMMA=`echo "$CHANNEL_LIST" | tr '\n' ',' | sed -e 's/ /,/g' -e 's/,$//'`
|
|
|
- g.message -d message="CHANNEL_LIST_COMMA=[$CHANNEL_LIST_COMMA]"
|
|
|
- g.remove rast="$CHANNEL_LIST_COMMA" --quiet
|
|
|
-
|
|
|
-}
|
|
|
-
|
|
|
-nowarpimport () {
|
|
|
- FILE=$1
|
|
|
- MAP=$2
|
|
|
-
|
|
|
- r.in.gdal -o $FLAGS input="$FILE" output="$MAP" --quiet
|
|
|
- if [ $? -ne 0 ] ; then
|
|
|
- g.message -e message="`basename $0`: r.in.gdal failure."
|
|
|
- rm -f "${TMP}"*
|
|
|
- exit 1
|
|
|
- fi
|
|
|
-
|
|
|
- #Get a list of channels:
|
|
|
- PATTERN="$MAP*"
|
|
|
- g.message -d message="pattern: [$PATTERN]"
|
|
|
- CHANNEL_LIST=`g.mlist type=rast pattern=$PATTERN`
|
|
|
- g.message -d message="channel list: [$CHANNEL_LIST]"
|
|
|
-
|
|
|
- # Get a list of suffixes:
|
|
|
- echo $CHANNEL_LIST
|
|
|
- CHANNEL_SUFFIXES=`echo "$CHANNEL_LIST" | sed "s/^${MAP}/sfx=/"`
|
|
|
- echo $CHANNEL_SUFFIXES
|
|
|
- # test for single band data
|
|
|
- if [ "$CHANNEL_SUFFIXES" = 'sfx=' ] ; then
|
|
|
- CHANNEL_SUFFIXES=""
|
|
|
- fi
|
|
|
- g.message -d message="channel suffixes: [$CHANNEL_SUFFIXES]"
|
|
|
-
|
|
|
- # Add to the list of all suffixes:
|
|
|
- echo $SUFFIXES
|
|
|
- SUFFIXES=`echo "$SUFFIXES
|
|
|
-$CHANNEL_SUFFIXES" | sort -u`
|
|
|
-
|
|
|
- IFS="$defaultIFS"
|
|
|
- echo $CHANNEL_SUFFIXES
|
|
|
- for SUFFIX in $CHANNEL_SUFFIXES ; do
|
|
|
- echo 'x', $sfx
|
|
|
- eval "$SUFFIX"
|
|
|
-
|
|
|
- # skip if only one band ($sfx is empty)
|
|
|
- if [ -z "$sfx" ] ; then
|
|
|
- g.message -d "Single band data"
|
|
|
- continue
|
|
|
- fi
|
|
|
-
|
|
|
- # Make patch lists:
|
|
|
- sfx2=`echo "$sfx" | sed "s/\./_/"`
|
|
|
- # This is a hack to make the patch lists empty:
|
|
|
- if [ $TITER -eq 0 ] ; then
|
|
|
- eval "PATCHES_$sfx2=\"\""
|
|
|
- fi
|
|
|
- eval "PATCHES_$sfx2=\${PATCHES_$sfx2},${MAP}${sfx}"
|
|
|
- done
|
|
|
-}
|
|
|
-
|
|
|
-
|
|
|
-# Compute flags for r.in.gdal
|
|
|
-FLAGS=""
|
|
|
-if [ $GIS_FLAG_E -eq 1 ] ; then
|
|
|
- FLAGS="${FLAGS} -e"
|
|
|
-fi
|
|
|
-if [ $GIS_FLAG_K -eq 1 ] ; then
|
|
|
- FLAGS="${FLAGS} -k"
|
|
|
-fi
|
|
|
-
|
|
|
-# And options for gdalwarp:
|
|
|
-
|
|
|
-case "${GIS_OPT_METHOD}" in
|
|
|
- "nearest") METHOD="-rn"
|
|
|
- ;;
|
|
|
- "bilinear") METHOD="-rb"
|
|
|
- ;;
|
|
|
- "cubic") METHOD="-rc"
|
|
|
- ;;
|
|
|
- "cubicspline") METHOD="-rcs"
|
|
|
- ;;
|
|
|
-esac
|
|
|
-
|
|
|
-IFS=,
|
|
|
-SUFFIXES=""
|
|
|
-# We need a way to make sure patches are intialized correctly
|
|
|
-TITER=0
|
|
|
-# map list is for single band patching
|
|
|
-MAPLIST=""
|
|
|
-
|
|
|
-# Import all the tiles:
|
|
|
-for input in $GIS_OPT_INPUT ; do
|
|
|
- IFS="$defaultIFS"
|
|
|
- TMPTILENAME="${GIS_OPT_OUTPUT}_tile_$TITER"
|
|
|
- if [ -f "$input" ] ; then
|
|
|
- if [ $GIS_FLAG_P -eq 1 ] ; then
|
|
|
- nowarpimport "$input" "$TMPTILENAME"
|
|
|
- else
|
|
|
- warpimport "$input" "$TMPTILENAME"
|
|
|
- fi
|
|
|
-
|
|
|
- if [ "$TITER" -eq 0 ] ; then
|
|
|
- MAPLIST="$TMPTILENAME"
|
|
|
- else
|
|
|
- MAPLIST="$MAPLIST,$TMPTILENAME"
|
|
|
- fi
|
|
|
-
|
|
|
- TITER="`expr $TITER + 1`"
|
|
|
- else
|
|
|
- g.message -w message="Missing input $input"
|
|
|
- fi
|
|
|
-done
|
|
|
-NUM_TILES="$TITER"
|
|
|
-
|
|
|
-# If there's more than one tile patch them together, otherwise just rename that tile.
|
|
|
-if [ $NUM_TILES -eq 1 ] ; then
|
|
|
- if [ -n "$CHANNEL_SUFFIXES" ] ; then
|
|
|
- for SUFFIX in $CHANNEL_SUFFIXES ; do
|
|
|
- eval "$SUFFIX"
|
|
|
- #Rename tile 0 to be the output
|
|
|
- GRASS_VERBOSE=1 \
|
|
|
- g.rename rast="${GIS_OPT_OUTPUT}_tile_0${sfx}","${GIS_OPT_OUTPUT}${sfx}"
|
|
|
- done
|
|
|
- else # single-band, single-tile
|
|
|
- echo 'x', $sfx
|
|
|
- #Rename tile 0 to be the output
|
|
|
- GRASS_VERBOSE=1 \
|
|
|
- g.rename rast="${GIS_OPT_OUTPUT}_tile_0${sfx}","${GIS_OPT_OUTPUT}${sfx}" --quiet
|
|
|
- fi
|
|
|
-else
|
|
|
- # Patch together each channel:
|
|
|
- g.message -d message="SUFFIXES: [$SUFFIXES]"
|
|
|
-
|
|
|
- if [ -n "$SUFFIXES" ] ; then
|
|
|
- # multi-band data
|
|
|
- IFS="$defaultIFS"
|
|
|
- for SUFFIX in $SUFFIXES ; do
|
|
|
- eval "$SUFFIX"
|
|
|
- sfx2=`echo "$sfx" | sed "s/\./_/"`
|
|
|
- eval "PATCHES=\${PATCHES_$sfx2}"
|
|
|
- # Patch these together (using nulls only):
|
|
|
- g.message "Patching [$sfx] channel"
|
|
|
- r.patch input="$PATCHES" output="${GIS_OPT_OUTPUT}${sfx}"
|
|
|
- # Remove the temporary patches we made
|
|
|
- g.remove rast="$PATCHES" --quiet
|
|
|
- done
|
|
|
- else
|
|
|
- # single-band data
|
|
|
- g.message "Patching tiles (this may take some time)"
|
|
|
- g.message -d message="patch list = [$MAPLIST]"
|
|
|
-
|
|
|
- # HACK: for 8bit PNG, GIF all the different tiles can have
|
|
|
- # different color tables so patches output end up all freaky.
|
|
|
- # r.mapcalc r#,g#,b# manual patching + r.composite?
|
|
|
- # or d.out.file + r.in.png + r.region?
|
|
|
- # *** there is a GDAL utility to do this: gdal_merge.py
|
|
|
- # http://www.gdal.org/gdal_merge.html
|
|
|
-
|
|
|
- #r.patch input="$MAPLIST" output="${GIS_OPT_OUTPUT}"
|
|
|
- IFS=" "
|
|
|
- for COLOR in r g b ; do
|
|
|
- for MAP in `echo "$MAPLIST" | tr ',' ' '` ; do
|
|
|
- GRASS_VERBOSE=0 \
|
|
|
- r.mapcalc "${MAP}_${COLOR} = ${COLOR}#$MAP"
|
|
|
- r.colors map="${MAP}_${COLOR}" color="grey255" --quiet
|
|
|
- r.null "${MAP}_${COLOR}" setnull=255 --quiet # so patching works
|
|
|
- done
|
|
|
-
|
|
|
- GRASS_VERBOSE=1\
|
|
|
- r.patch input=`echo "$MAPLIST" | sed -e "s/,/_${COLOR},/g" -e "s/$/_${COLOR}/"` \
|
|
|
- output="${MAP}_${COLOR}_all"
|
|
|
- done
|
|
|
- r.composite red="${MAP}_r_all" green="${MAP}_g_all" blue="${MAP}_b_all" \
|
|
|
- output="${GIS_OPT_OUTPUT}"
|
|
|
-
|
|
|
- if [ $? -eq 0 ] ; then
|
|
|
- IFS=','
|
|
|
- for MAP in $MAPLIST ; do
|
|
|
- g.mremove -f rast="${MAP}*" --quiet
|
|
|
- done
|
|
|
- fi
|
|
|
- fi
|
|
|
-fi
|
|
|
-
|
|
|
-# Set up color counter
|
|
|
-COLORS=0
|
|
|
-
|
|
|
-IFS="$defaultIFS"
|
|
|
-for SUFFIX in $SUFFIXES ; do
|
|
|
- eval "$SUFFIX"
|
|
|
- # Keep track of colors
|
|
|
- if [ -z "$sfx" ] ; then
|
|
|
- # There's already a no suffix, can't make colors
|
|
|
- COLORS=4 # Can only go up from here ;)
|
|
|
- fi
|
|
|
- if [ "$sfx" = ".red" ] || [ "$sfx" = ".green" ] || [ "$sfx" = ".blue" ]; then
|
|
|
- COLORS="`expr $COLORS + 1`"
|
|
|
- fi
|
|
|
-done
|
|
|
-
|
|
|
-echo $COLORS
|
|
|
-# Make a composite image if asked for and colors exist
|
|
|
-if [ $COLORS -eq 3 ] && [ $GIS_FLAG_C -eq 1 ] ; then
|
|
|
- g.message "Building Color Image"
|
|
|
- r.composite red="${GIS_OPT_OUTPUT}.red" green="${GIS_OPT_OUTPUT}.green" \
|
|
|
- blue="${GIS_OPT_OUTPUT}.blue" output="${GIS_OPT_OUTPUT}"
|
|
|
- g.message message="Written: ${GIS_OPT_OUTPUT}"
|
|
|
-fi
|
|
|
-
|
|
|
-# Clean up:
|
|
|
-rm -f "${TMP}"*
|