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 |
---|
18 | data_suffix=".xyz" # Merged point cloud |
---|
19 | hist_suffix=".hist" # Histogram |
---|
20 | mask_suffix=".mask" # Mask points |
---|
21 | |
---|
22 | temp=/tmp # Temp directory - must exist |
---|
23 | detail_file=$temp"/lidar"$hist_suffix |
---|
24 | data_file=$temp"/lidar"$data_suffix |
---|
25 | mask_file=$temp"/lidar"$mask_suffix |
---|
26 | dirfile=$temp"/dirlist.txt" |
---|
27 | |
---|
28 | USEMASK=0 |
---|
29 | REUSE=0 |
---|
30 | PADDING=2000 # Amount in m to pad boundary points |
---|
31 | |
---|
32 | usage() |
---|
33 | { |
---|
34 | cat << eof |
---|
35 | |
---|
36 | usage: $0 arguments |
---|
37 | |
---|
38 | Converts a set of Lidar point cloud files into an azgcorr-style DEM using GRASS |
---|
39 | |
---|
40 | Arguments: |
---|
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) |
---|
50 | eof |
---|
51 | } |
---|
52 | |
---|
53 | while getopts âhg:b:mnroâ OPTION |
---|
54 | do |
---|
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 |
---|
86 | done |
---|
87 | |
---|
88 | if [ -z $GRASSDIR ] ; then |
---|
89 | echo "Missing GRASS mapset path" |
---|
90 | usage |
---|
91 | exit 1 |
---|
92 | fi |
---|
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 |
---|
96 | if [ $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 |
---|
143 | fi |
---|
144 | |
---|
145 | # Read bounds from histogram, add PADDING m around edges |
---|
146 | echo "Reading histogram..." |
---|
147 | North_Max=$((`cat $detail_file | grep '^n:' | cut -c 4-` + PADDING)) |
---|
148 | South_Max=$((`cat $detail_file | grep '^s:' | cut -c 4-` - PADDING)) |
---|
149 | East_Max=$((`cat $detail_file | grep '^e:' | cut -c 4-` + PADDING)) |
---|
150 | West_Max=$((`cat $detail_file | grep '^w:' | cut -c 4-` - PADDING)) |
---|
151 | Ubound=`cat $detail_file | grep '^Upper:' | cut -c 8-` |
---|
152 | Lbound=`cat $detail_file | grep '^Lower:' | cut -c 8-` |
---|
153 | |
---|
154 | echo "North="$North_Max |
---|
155 | echo "South="$South_Max |
---|
156 | echo "East="$East_Max |
---|
157 | echo "West="$West_Max |
---|
158 | echo "Upper="$Ubound |
---|
159 | echo "Lower="$Lbound |
---|
160 | |
---|
161 | Base_Name=`basename $data_file` |
---|
162 | Map_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 |
---|
165 | echo "Running GRASS..." |
---|
166 | if [ $USEMASK = -1 ] ; then |
---|
167 | grass62 -text $GRASSDIR << eof |
---|
168 | g.region n=$North_Max s=$South_Max e=$East_Max w=$West_Max res=5 |
---|
169 | r.in.xyz input=$data_file output=${Map_Base_Name}_raw x=2 y=3 z=4 fs=" " zrange=${Lbound},${Ubound} --overwrite |
---|
170 | r.surf.idw input=${Map_Base_Name}_raw output=${Map_Base_Name} --overwrite |
---|
171 | r.out.ascii input=${Map_Base_Name} output=${Map_Base_Name}_conv.dem null=0 |
---|
172 | eof |
---|
173 | else |
---|
174 | grass62 -text $GRASSDIR << eof |
---|
175 | g.region n=$North_Max s=$South_Max e=$East_Max w=$West_Max res=5 |
---|
176 | v.in.ascii input=$mask_file output=v_${Map_Base_Name}_points format=point fs=, x=1 y=2 --overwrite |
---|
177 | v.hull -a input=v_${Map_Base_Name}_points output=v_${Map_Base_Name}_poly --overwrite |
---|
178 | v.to.rast input=v_${Map_Base_Name}_poly output=${Map_Base_Name}_mask use=val value=1 --overwrite |
---|
179 | r.mask -o input=${Map_Base_Name}_mask maskcats=* |
---|
180 | r.in.xyz input=$data_file output=${Map_Base_Name}_raw x=2 y=3 z=4 fs=" " zrange=${Lbound},${Ubound} --overwrite |
---|
181 | r.surf.idw input=${Map_Base_Name}_raw output=${Map_Base_Name} --overwrite |
---|
182 | r.out.ascii input=${Map_Base_Name} output=${Map_Base_Name}_conv.dem null=0 |
---|
183 | r.mask -r input=${Map_Base_Name}_mask |
---|
184 | eof |
---|
185 | fi |
---|
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 | |
---|
191 | echo "Converting DEM format..." |
---|
192 | demheaderconvert.sh ${Map_Base_Name}_conv.dem ${Map_Base_Name}.dem |
---|
193 | |
---|
194 | echo "Deleting temporary files..." |
---|
195 | rm -f $dirfile |
---|
196 | if [ ! $BASEFILE ] ; then |
---|
197 | rm -f $data_file |
---|
198 | rm -f $detail_file |
---|
199 | rm -f $mask_file |
---|
200 | fi |
---|
201 | |
---|
202 | echo "Lidar conversion complete!" |
---|