Processing/LIDARDEMs: lidar2dem.sh

# Shell script for generating an Azgcorr-format ascii DEM from LIDAR point-cloud data
#
# Author: Ben Taylor
# Date: 21/04/2008
#
# Arguments:
# $1: GRASS mapset directory
# $2: Intermediate file name (optional but recommended - can be reused later and GRASS map name is based on it)
# $3: --reuse to reuse intermediate file (assuming it already exists) or --overwrite to overwrite it if it already exists (optional)
#
# Usage/argument handling derived from http://rsalveti.wordpress.com/2007/04/03/bash-parsing-arguments-with-getopts/ under a
# Creative Commons Attribution-Share Alike license (see http://creativecommons.org/licenses/by-sa/2.5/br/deed.en_GB).
# Remainder of script copyright Plymouth Marine Laboratory 2008, available under the same license - you are free to use
# and redistribute as you wish, provided the author is credited and any derivative works are distributed only under the same license.
# No warranty of any kind is given.

# Suffixes - intermediate files are the given basename (or default) with these suffixes
data_suffix=".xyz" # Merged point cloud
hist_suffix=".hist" # Histogram
mask_suffix=".mask" # Mask points

temp=/tmp # Temp directory - must exist
detail_file=$temp"/lidar"$hist_suffix
data_file=$temp"/lidar"$data_suffix
mask_file=$temp"/lidar"$mask_suffix
dirfile=$temp"/dirlist.txt"

USEMASK=0
REUSE=0
PADDING=2000 # Amount in m to pad boundary points

usage()
{
cat << eof

usage: $0 arguments

Converts a set of Lidar point cloud files into an azgcorr-style DEM using GRASS

Arguments:
    -h:         Display usage
    -g <path>:  Full path to GRASS mapset directory to be used (required)
    -b <name>:  Base name for existing data files (data, histogram, point, mask)
                Should be set to identifying name for this data set.
                Optional but recommended (so they can be re-used later if desired)
    -m:         Create and use a GRASS mask (default)
    -n:         Do not use GRASS mask
    -r:         Reuse existing intermediate files (optional, default)
    -o:         Overwrite any existing mask/point files (optional)
eof
}

while getopts “hg:b:mnro” OPTION
do
    case $OPTION in
        h)
            usage
            exit 1
            ;;
        g)
            GRASSDIR=$OPTARG
            ;;
        b)
            BASEFILE=$OPTARG
            ;;
        r)
            REUSE=1
            echo "Reusing existing files"
            ;;
        o)
            REUSE=-1
            echo "Overwriting any existing files"
            ;;
        m)
            USEMASK=1
            ;;
        n)
            USEMASK=-1
            echo "Not using a GRASS mask"
            ;;
        ?)
            usage
            exit 1
            ;;
    esac
done

if [ -z $GRASSDIR ] ; then
    echo "Missing GRASS mapset path"
    usage
    exit 1
fi

# If a basefile name was supplied, need to check if we're reusing existing intermediate files ($REUSE=-1)
# If we're not (or they don't already exist) then create them from scratch
if [ $BASEFILE ] ; then
    
    detail_file=${BASEFILE}${hist_suffix}
    mask_file=${BASEFILE}${mask_suffix}
    data_file=${BASEFILE}${data_suffix}
    
    if [ $REUSE = -1 ] ; then
        if [ $USEMASK = -1 ] ; then
            echo "Merging raw data files (creating "$data_file")..."
            trim_lidar.sh > $data_file
            echo "Not generating GRASS mask"
        else
            echo "Merging raw data files (creating "$data_file") and generating mask ("$mask_file")..."
            \ls -1 *.all.bz2 > $dirfile
            generate_mask.sh $dirfile $mask_file $data_file
        fi
        
        echo "Generating histogram ("$detail_file")..."
        lidar_histogram.sh $data_file | tee $detail_file
    else
        if [ -e $data_file ] && [ -e $mask_file ] && [ $USEMASK != -1 ] ; then
            echo "Reading existing data file: "$data_file
            echo "Reading existing mask file: "$mask_file
        else
            if [ -e $data_file ] && [ $USEMASK = -1 ] ; then
                echo "Reading existing data file: "$data_file
                echo "Not generating GRASS mask"
            else
                if [ $USEMASK = -1 ] ; then
                    echo "Merging raw data files (creating "$data_file")..."
                    trim_lidar.sh > $data_file
                    echo "Not generating GRASS mask"
                else
                    echo "Merging raw data files (creating "$data_file") and generating mask ("$mask_file")..."
                    \ls -1 *.all.bz2 > $dirfile
                    generate_mask.sh $dirfile $mask_file $data_file
                fi
            fi
        fi
        
        if [ -e $detail_file ] ; then
            echo "Reading existing histogram file: "$detail_file
        else
            echo "Generating histogram ("$detail_file")..."
            lidar_histogram.sh $data_file | tee $detail_file
        fi
    fi
fi

# Read bounds from histogram, add PADDING m around edges
echo "Reading histogram..."
North_Max=$((`cat $detail_file | grep '^n:' | cut -c 4-` + PADDING))
South_Max=$((`cat $detail_file | grep '^s:' | cut -c 4-` - PADDING))
East_Max=$((`cat $detail_file | grep '^e:' | cut -c 4-` + PADDING))
West_Max=$((`cat $detail_file | grep '^w:' | cut -c 4-` - PADDING))
Ubound=`cat $detail_file | grep '^Upper:' | cut -c 8-`
Lbound=`cat $detail_file | grep '^Lower:' | cut -c 8-`

echo "North="$North_Max
echo "South="$South_Max
echo "East="$East_Max
echo "West="$West_Max
echo "Upper="$Ubound
echo "Lower="$Lbound

Base_Name=`basename $data_file`
Map_Base_Name=`echo $Base_Name | sed 's/\(.\+\)\(\.[^.]\+\)/\1/'` # Chop the .txt (or whatever) off the end of the intermediate file

# Grass commands are different depending on whether you're using a mask or not
echo "Running GRASS..."
if [ $USEMASK = -1 ] ; then
grass62 -text $GRASSDIR << eof
g.region n=$North_Max s=$South_Max e=$East_Max w=$West_Max res=5
r.in.xyz input=$data_file output=${Map_Base_Name}_raw x=2 y=3 z=4 fs=" " zrange=${Lbound},${Ubound} --overwrite
r.surf.idw input=${Map_Base_Name}_raw output=${Map_Base_Name} --overwrite
r.out.ascii input=${Map_Base_Name} output=${Map_Base_Name}_conv.dem null=0
eof
else
grass62 -text $GRASSDIR << eof
g.region n=$North_Max s=$South_Max e=$East_Max w=$West_Max res=5
v.in.ascii input=$mask_file output=v_${Map_Base_Name}_points format=point fs=, x=1 y=2 --overwrite
v.hull -a input=v_${Map_Base_Name}_points output=v_${Map_Base_Name}_poly --overwrite
v.to.rast input=v_${Map_Base_Name}_poly output=${Map_Base_Name}_mask use=val value=1 --overwrite
r.mask -o input=${Map_Base_Name}_mask maskcats=*
r.in.xyz input=$data_file output=${Map_Base_Name}_raw x=2 y=3 z=4 fs=" " zrange=${Lbound},${Ubound} --overwrite
r.surf.idw input=${Map_Base_Name}_raw output=${Map_Base_Name} --overwrite
r.out.ascii input=${Map_Base_Name} output=${Map_Base_Name}_conv.dem null=0
r.mask -r input=${Map_Base_Name}_mask
eof
fi

# Alternative interpolation method (uses vectors - slower, uses more memory) - to use, use line below in place
# of the r.surf.idw lines above
#r.fillnulls input=${Map_Base_Name}_raw output=${Map_Base_Name} --overwrite

echo "Converting DEM format..."
demheaderconvert.sh ${Map_Base_Name}_conv.dem ${Map_Base_Name}.dem

echo "Deleting temporary files..."
rm -f $dirfile
if [ ! $BASEFILE ] ; then
    rm -f $data_file
    rm -f $detail_file
    rm -f $mask_file
fi

echo "Lidar conversion complete!"