Processing/LIDARDEMs: lidar2dem.sh

File lidar2dem.sh, 7.2 KB (added by benj, 11 years ago)
Line 
1# Shell script for generating an Azgcorr-format ascii DEM from LIDAR point-cloud data
2#
3# Author: Ben Taylor
4# Date: 21/04/2008
5#
6# Arguments:
7# $1: GRASS mapset directory
8# $2: Intermediate file name (optional but recommended - can be reused later and GRASS map name is based on it)
9# $3: --reuse to reuse intermediate file (assuming it already exists) or --overwrite to overwrite it if it already exists (optional)
10#
11# Usage/argument handling derived from http://rsalveti.wordpress.com/2007/04/03/bash-parsing-arguments-with-getopts/ under a
12# Creative Commons Attribution-Share Alike license (see http://creativecommons.org/licenses/by-sa/2.5/br/deed.en_GB).
13# Remainder of script copyright Plymouth Marine Laboratory 2008, available under the same license - you are free to use
14# and redistribute as you wish, provided the author is credited and any derivative works are distributed only under the same license.
15# No warranty of any kind is given.
16
17# Suffixes - intermediate files are the given basename (or default) with these suffixes
18data_suffix=".xyz" # Merged point cloud
19hist_suffix=".hist" # Histogram
20mask_suffix=".mask" # Mask points
21
22temp=/tmp # Temp directory - must exist
23detail_file=$temp"/lidar"$hist_suffix
24data_file=$temp"/lidar"$data_suffix
25mask_file=$temp"/lidar"$mask_suffix
26dirfile=$temp"/dirlist.txt"
27
28USEMASK=0
29REUSE=0
30PADDING=2000 # Amount in m to pad boundary points
31
32usage()
33{
34cat << eof
35
36usage: $0 arguments
37
38Converts a set of Lidar point cloud files into an azgcorr-style DEM using GRASS
39
40Arguments:
41    -h:         Display usage
42    -g <path>:  Full path to GRASS mapset directory to be used (required)
43    -b <name>:  Base name for existing data files (data, histogram, point, mask)
44                Should be set to identifying name for this data set.
45                Optional but recommended (so they can be re-used later if desired)
46    -m:         Create and use a GRASS mask (default)
47    -n:         Do not use GRASS mask
48    -r:         Reuse existing intermediate files (optional, default)
49    -o:         Overwrite any existing mask/point files (optional)
50eof
51}
52
53while getopts “hg:b:mnro” OPTION
54do
55    case $OPTION in
56        h)
57            usage
58            exit 1
59            ;;
60        g)
61            GRASSDIR=$OPTARG
62            ;;
63        b)
64            BASEFILE=$OPTARG
65            ;;
66        r)
67            REUSE=1
68            echo "Reusing existing files"
69            ;;
70        o)
71            REUSE=-1
72            echo "Overwriting any existing files"
73            ;;
74        m)
75            USEMASK=1
76            ;;
77        n)
78            USEMASK=-1
79            echo "Not using a GRASS mask"
80            ;;
81        ?)
82            usage
83            exit 1
84            ;;
85    esac
86done
87
88if [ -z $GRASSDIR ] ; then
89    echo "Missing GRASS mapset path"
90    usage
91    exit 1
92fi
93
94# If a basefile name was supplied, need to check if we're reusing existing intermediate files ($REUSE=-1)
95# If we're not (or they don't already exist) then create them from scratch
96if [ $BASEFILE ] ; then
97   
98    detail_file=${BASEFILE}${hist_suffix}
99    mask_file=${BASEFILE}${mask_suffix}
100    data_file=${BASEFILE}${data_suffix}
101   
102    if [ $REUSE = -1 ] ; then
103        if [ $USEMASK = -1 ] ; then
104            echo "Merging raw data files (creating "$data_file")..."
105            trim_lidar.sh > $data_file
106            echo "Not generating GRASS mask"
107        else
108            echo "Merging raw data files (creating "$data_file") and generating mask ("$mask_file")..."
109            \ls -1 *.all.bz2 > $dirfile
110            generate_mask.sh $dirfile $mask_file $data_file
111        fi
112       
113        echo "Generating histogram ("$detail_file")..."
114        lidar_histogram.sh $data_file | tee $detail_file
115    else
116        if [ -e $data_file ] && [ -e $mask_file ] && [ $USEMASK != -1 ] ; then
117            echo "Reading existing data file: "$data_file
118            echo "Reading existing mask file: "$mask_file
119        else
120            if [ -e $data_file ] && [ $USEMASK = -1 ] ; then
121                echo "Reading existing data file: "$data_file
122                echo "Not generating GRASS mask"
123            else
124                if [ $USEMASK = -1 ] ; then
125                    echo "Merging raw data files (creating "$data_file")..."
126                    trim_lidar.sh > $data_file
127                    echo "Not generating GRASS mask"
128                else
129                    echo "Merging raw data files (creating "$data_file") and generating mask ("$mask_file")..."
130                    \ls -1 *.all.bz2 > $dirfile
131                    generate_mask.sh $dirfile $mask_file $data_file
132                fi
133            fi
134        fi
135       
136        if [ -e $detail_file ] ; then
137            echo "Reading existing histogram file: "$detail_file
138        else
139            echo "Generating histogram ("$detail_file")..."
140            lidar_histogram.sh $data_file | tee $detail_file
141        fi
142    fi
143fi
144
145# Read bounds from histogram, add PADDING m around edges
146echo "Reading histogram..."
147North_Max=$((`cat $detail_file | grep '^n:' | cut -c 4-` + PADDING))
148South_Max=$((`cat $detail_file | grep '^s:' | cut -c 4-` - PADDING))
149East_Max=$((`cat $detail_file | grep '^e:' | cut -c 4-` + PADDING))
150West_Max=$((`cat $detail_file | grep '^w:' | cut -c 4-` - PADDING))
151Ubound=`cat $detail_file | grep '^Upper:' | cut -c 8-`
152Lbound=`cat $detail_file | grep '^Lower:' | cut -c 8-`
153
154echo "North="$North_Max
155echo "South="$South_Max
156echo "East="$East_Max
157echo "West="$West_Max
158echo "Upper="$Ubound
159echo "Lower="$Lbound
160
161Base_Name=`basename $data_file`
162Map_Base_Name=`echo $Base_Name | sed 's/\(.\+\)\(\.[^.]\+\)/\1/'` # Chop the .txt (or whatever) off the end of the intermediate file
163
164# Grass commands are different depending on whether you're using a mask or not
165echo "Running GRASS..."
166if [ $USEMASK = -1 ] ; then
167grass62 -text $GRASSDIR << eof
168g.region n=$North_Max s=$South_Max e=$East_Max w=$West_Max res=5
169r.in.xyz input=$data_file output=${Map_Base_Name}_raw x=2 y=3 z=4 fs=" " zrange=${Lbound},${Ubound} --overwrite
170r.surf.idw input=${Map_Base_Name}_raw output=${Map_Base_Name} --overwrite
171r.out.ascii input=${Map_Base_Name} output=${Map_Base_Name}_conv.dem null=0
172eof
173else
174grass62 -text $GRASSDIR << eof
175g.region n=$North_Max s=$South_Max e=$East_Max w=$West_Max res=5
176v.in.ascii input=$mask_file output=v_${Map_Base_Name}_points format=point fs=, x=1 y=2 --overwrite
177v.hull -a input=v_${Map_Base_Name}_points output=v_${Map_Base_Name}_poly --overwrite
178v.to.rast input=v_${Map_Base_Name}_poly output=${Map_Base_Name}_mask use=val value=1 --overwrite
179r.mask -o input=${Map_Base_Name}_mask maskcats=*
180r.in.xyz input=$data_file output=${Map_Base_Name}_raw x=2 y=3 z=4 fs=" " zrange=${Lbound},${Ubound} --overwrite
181r.surf.idw input=${Map_Base_Name}_raw output=${Map_Base_Name} --overwrite
182r.out.ascii input=${Map_Base_Name} output=${Map_Base_Name}_conv.dem null=0
183r.mask -r input=${Map_Base_Name}_mask
184eof
185fi
186
187# Alternative interpolation method (uses vectors - slower, uses more memory) - to use, use line below in place
188# of the r.surf.idw lines above
189#r.fillnulls input=${Map_Base_Name}_raw output=${Map_Base_Name} --overwrite
190
191echo "Converting DEM format..."
192demheaderconvert.sh ${Map_Base_Name}_conv.dem ${Map_Base_Name}.dem
193
194echo "Deleting temporary files..."
195rm -f $dirfile
196if [ ! $BASEFILE ] ; then
197    rm -f $data_file
198    rm -f $detail_file
199    rm -f $mask_file
200fi
201
202echo "Lidar conversion complete!"