-
Notifications
You must be signed in to change notification settings - Fork 2
/
Lidar2HemiEval_Prep.m
243 lines (183 loc) · 8.95 KB
/
Lidar2HemiEval_Prep.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
function Lidar2HemiEval_Prep(sfile)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % GENERAL DESCRIPTION
% % Preparatory script that generates canopy height models and gathers
% % input and runs L2HE_Prep_R.R
% %
% % VERSION
% % v1.0 updated 06.07.2020
% %
% % AUTHOR:
% % Clare Webster (1,2)
% % (1) WSL Institute for Snow and Avalanche Research SLF, Davos, CH
% % (2) University of Edinburgh, School of GeoSciences, Edinburgh, UK
% %
% % CONTRIBUTING MATERIAL
% % Calculation of canopy height models follows
% % 'generating_a_pit_free_chm.bat' from LAStools, based on
% % http://www.riegl.com/uploads/tx_pxpriegldownloads/khosravipour_SilviLaser2013.pdf
% %
% % USAGE
% % > Lidar2HemiEval_Prep('L2HEPrep_settings.m')
% %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOAD SETTINGS
if nargin == 0
error('No settings file specified')
end
[setp,setf,~] = fileparts(sfile);
cd(setp);
try
if exist(sfile,'file')
prepset = feval(setf);
else
error('Specified settings folder is not in directory')
end
catch
error('Error while reading L2HEPrep settings file')
end
% Check input compatability
if prepset.in.gengrid ~= 0 && prepset.in.gengrid ~= 1
error('Unknown setting for grid generation. Options are 0 or 1 to disable or enable functionality')
end
if ~isempty(prepset.out.chmfmt) && ~strcmp(prepset.out.chmfmt,'.tif') && ~strcmp(prepset.out.chmfmt,'.asc')
error('Unknown CHM output format. Options are .tif or .asc or empty')
end
if (strcmp(prepset.out.chmfmt,'.asc') || isempty(prepset.out.chmfmt)) &&...
(prepset.in.trunks || prepset.in.branches)
disp('Warning: CHM output format changed from to .tif to enable calculation of trunks and branches')
prepset.in.chmfmt = '.tif';
end
if prepset.in.biome < 0 || prepset.in.biome > 24
error('Unknown input biome for calculating diameter at breast height')
end
% Species settings
switch prepset.in.species
case 'NS'
species = 1;
case 'DF'
species = 2;
case 'SP'
species = 3;
case 'SB'
species = 4;
case 'NA'
species = 0;
otherwise
error('Input species not recognised');
end
% Check spacing input
if isempty(prepset.in.spacing) && prepset.in.gengrid == 1
disp('Grid spacing left empty. Using default of 1 metre')
prepset.in.spacing = 1;
elseif isempty(prepset.in.spacing) && prepset.in.gengrid == 0
prepset.in.spacing = 1; % will be unused but variable needs value for R script to run
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Clip large lidar to analysis area + buffer
tmpdir = fullfile(prepset.in.basefolder,'temp');
if ~exist(tmpdir,'dir')
mkdir(tmpdir)
end
tx1 = prepset.in.xlimits(1);
tx2 = prepset.in.xlimits(2);
ty1 = prepset.in.ylimits(1);
ty2 = prepset.in.ylimits(2);
[~,lasname] = fileparts(prepset.in.lasfile);
if prepset.in.cliplas
new(1,:) = [tx1-prepset.in.buffer,ty1-prepset.in.buffer];
new(2,:) = [tx2+prepset.in.buffer,ty1-prepset.in.buffer];
new(3,:) = [tx2+prepset.in.buffer,ty2+prepset.in.buffer];
new(4,:) = [tx1-prepset.in.buffer,ty2+prepset.in.buffer];
new(5,:) = [tx1-prepset.in.buffer,ty1-prepset.in.buffer];
clipfilename = fullfile(tmpdir,strcat(lasname,'_clipdim.txt'));
dlmwrite(clipfilename,new,'Precision',12,'Delimiter','\t');
fprintf('\n Clipping Lidar.... \n \n')
lasname = strcat(lasname,'_clipped');
lasfilename = fullfile(prepset.in.basefolder,'Data_Surface','DSM',strcat(lasname,'.las'));
system([prepset.in.lastoolpath 'lasclip.exe -i ' prepset.in.lasfile ' -poly '...
clipfilename ' -o ' lasfilename ]);
else
lasfilename = prepset.in.lasfile;
end
if prepset.in.gengrid
outdir = fullfile(prepset.in.basefolder,'Output_Prep',lasname);
if ~exist(outdir,'dir'); mkdir(outdir); end
dlmwrite(fullfile(outdir,strcat(lasname,'_analysisarea.txt')),[tx1,tx2,ty1,ty2],'Precision',12,'Delimiter','\t')
end
[chmfname, normlidar] = create_CHM(lasfilename,prepset,tmpdir);
if prepset.in.trunks || prepset.in.branches
if prepset.in.trunks; tdir = fullfile(prepset.in.basefolder,'Data_Surface','DBH'); if ~exist(tdir,'dir'); mkdir(tdir); end; end
if prepset.in.branches; bdir = fullfile(prepset.in.basefolder,'Data_Surface','LTC'); if ~exist(bdir,'dir'); mkdir(bdir); end; end
fprintf('\n Starting segmentation in R.... \n \n')
try
system([prepset.in.rpath 'Rscript.exe ' prepset.in.preppath ' ' prepset.in.basefolder ' '...
chmfname ' ' normlidar ' ' lasname ' ' num2str(prepset.in.trunks) ' ' num2str(prepset.in.branches) ' ' ...
num2str(prepset.in.gengrid) ' ' num2str(prepset.in.spacing) ' ' ...
prepset.in.epsgstr ' ' num2str(prepset.in.biome) ' ' num2str(species)]);
catch
error('Failed to execute R segmentation script. Check R file path / version number')
end
r_finfile = fullfile(tmpdir,'R_finish_confirmation.txt');
if exist(r_finfile,'file')
fprintf('\n Segmentation in R finished successfully.... \n')
delete(r_finfile) % then delete it because it's not needed
else
error('Failure in R segmentation script')
end
delete(normlidar) % deleted here and not in make CHM function because the segmentation algorithm needs it
end
% clear tmpdir
if exist(tmpdir,'dir')
rmdir(tmpdir,'s')
end
if ~prepset.out.saveCHM
delete(chmfname)
if ~prepset.in.gengrid
delete(fullfile(prepset.in.basefolder,'Output_Prep',lasname))
end
end
fprintf('\n Finished. \n \n')
end
function [chmfname, normlidar] = create_CHM(lasfilename,prepset,tmpdir)
if ~exist(tmpdir,'dir'); mkdir(tmpdir); end
[~,lasname] = fileparts(lasfilename);
outdir = fullfile(prepset.in.basefolder,'Output_Prep',lasname);
if ~exist(outdir,'dir'); mkdir(outdir); end
chmfname = fullfile(prepset.in.basefolder,'Output_Prep',lasname,strcat(lasname,'_chm',prepset.out.chmfmt));
grndlidar = fullfile(tmpdir,strcat(lasname,'_ground.laz'));
normlidar = fullfile(tmpdir,strcat(lasname,'_norm.laz'));
fprintf('\n Calculating ground points... \n')
system([prepset.in.lastoolpath 'lasground_new.exe -i ' lasfilename ' -wilderness -o ' grndlidar]);
fprintf('\n Making CHM.... \n \n')
%%% The following steps were translated into Matlab from
%%% 'generating_a_pit_free_chm.bat' from LAStools:
%%% :: a batch script for generating a pit-free CHM as outlined
%%% :: in the Silvilaser 2013 poster by A. Khosravipour et al.
% CHM settings
step = prepset.in.chmstep;
kill = 1;
% create normalised point cloud
system([prepset.in.lastoolpath 'lasheight.exe -i ' grndlidar ' -replace_z -o ' normlidar]);
% make the chm
system([prepset.in.lastoolpath 'blast2dem.exe -i ' normlidar ' -keep_first -step ' num2str(step) ...
' -odir ' tmpdir ' -odix _00 -obil']);
system([prepset.in.lastoolpath 'blast2dem.exe -i ' normlidar ' -keep_first -drop_z_below 2 -step ' ...
num2str(step) ' -kill ' num2str(kill) ' -odir ' tmpdir ' -odix _02 -obil']);
system([prepset.in.lastoolpath 'blast2dem.exe -i ' normlidar ' -keep_first -drop_z_below 5 -step ' ...
num2str(step) ' -kill ' num2str(kill) ' -odir ' tmpdir ' -odix _05 -obil']);
system([prepset.in.lastoolpath 'blast2dem.exe -i ' normlidar ' -keep_first -drop_z_below 10 -step ' ...
num2str(step) ' -kill ' num2str(kill) ' -odir ' tmpdir ' -odix _10 -obil']);
system([prepset.in.lastoolpath 'blast2dem.exe -i ' normlidar ' -keep_first -drop_z_below 15 -step ' ...
num2str(step) ' -kill ' num2str(kill) ' -odir ' tmpdir ' -odix _15 -obil']);
system([prepset.in.lastoolpath 'blast2dem.exe -i ' normlidar ' -keep_first -drop_z_below 20 -step ' ...
num2str(step) ' -kill ' num2str(kill) ' -odir ' tmpdir ' -odix _20 -obil']);
system([prepset.in.lastoolpath 'blast2dem.exe -i ' normlidar ' -keep_first -drop_z_below 25 -step ' ...
num2str(step) ' -kill ' num2str(kill) ' -odir ' tmpdir ' -odix _25 -obil']);
system([prepset.in.lastoolpath 'blast2dem.exe -i ' normlidar ' -keep_first -drop_z_below 30 -step ' ...
num2str(step) ' -kill ' num2str(kill) ' -odir ' tmpdir ' -odix _30 -obil']);
system([prepset.in.lastoolpath 'lasgrid -i ' fullfile(tmpdir,'*.bil') ' -merged -step ' ...
num2str(step) ' -highest -o ' chmfname]);
delete(grndlidar)
fprintf('\n ....done \n \n')
end