-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathautoCmaskDetectCirrus.m
112 lines (87 loc) · 4.19 KB
/
autoCmaskDetectCirrus.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
function autoCmaskDetectCirrus(varargin)
%AUTOCMASKDETECTCIRRUS Detect cirrus clouds by differing predicted and
%observed cirrus band TOA reflectance.
% Specific parameters
% ------------------------
% 'DirIn' Directory of input stacked data.
% 'DirBackground' Directory of backup coeffs.
% 'DirOut' Directory of output data.
% 'ThrdNormCirrus' Default value: 0.5. Usually this can be adjusted to decrease commisison or omission errors.
% 'ThrdMinCirrus' Default value: 31. See Qiu et al., RSE, 2020.
% History:
% 1. add new labels for no enough observations for building time series
% modle (254) and filled pixel for individual image (255). 16 Aug., 2020 Shi
% autoCmaskDetectCirrus('DirIn', 'C:\Users\qsly0\Desktop\Example_Data_Cmask\CmaskInput','DirBackground', 'C:\Users\qsly0\Desktop\Example_Data_Cmask\CmaskOutput\Background', 'DirOut', 'C:\Users\qsly0\Desktop\Example_Data_Cmask\CmaskOutput\CirrusMask')
thrd_normdiff = 0.5;
thrd_mincir = 31; % 10000
save_pcd = 0;
dir_work = pwd;
ds=DefaultSetting();
dir_cmask_bgd = fullfile(dir_work, ...
ds.foldername_output, ds.foldername_bg);
dir_cmask_inpt = fullfile(dir_work,ds.foldername_input);
dir_cmask_cirrus_mask = fullfile(dir_work, ...
ds.foldername_output, ds.foldername_mask);
dir_cmask_pred = fullfile(dir_work, ...
ds.foldername_output, ds.foldername_pred);
p = inputParser;
p.FunctionName = 'detectParas';
% optional
% default values.
addParameter(p,'DirIn',dir_cmask_inpt);
addParameter(p,'DirBackground',dir_cmask_bgd);
addParameter(p,'DirOut',dir_cmask_cirrus_mask);
addParameter(p,'ThrdNormCirrus',thrd_normdiff);
addParameter(p,'ThrdMinCirrus',thrd_mincir);
% request user's input
parse(p,varargin{:});
dir_cmask_inpt = p.Results.DirIn;
dir_cmask_bgd = p.Results.DirBackground;
dir_cmask_cirrus_mask = p.Results.DirOut;
thrd_normdiff = p.Results.ThrdNormCirrus;
thrd_mincir = p.Results.ThrdMinCirrus;
warning('off','all');
% load Cmask background image
path_cmask_bgd = dir(fullfile(dir_cmask_bgd,'L*_BG'));
[im_bgd,~,jiUL,resolu,ZC]=enviread(fullfile(dir_cmask_bgd, path_cmask_bgd.name));
% get image parameters automatically
imf = dir(fullfile(dir_cmask_inpt,'L*')); % folder names
for i =1:length(imf)
name_cmask_img = imf(i).name;
% load water vapor
path_cmask_input = fullfile(dir_cmask_inpt,name_cmask_img,[name_cmask_img,'_MTLstack']);
im = enviread(path_cmask_input);
img_obserd = single(im(:,:,1));
wvs = single(im(:,:,2))./100;% back to 100
clear im;
% convert to X
yr = str2num(name_cmask_img(10:13));
doy = str2num(name_cmask_img(14:16));
Xs = datenum(yr, 1, 0) + doy;
w = 2*pi/365.25;
img_pred = im_bgd(:,:,1) + ...
im_bgd(:,:,2).*cos(Xs.*w)+im_bgd(:,:,3).*sin(Xs.*w)+...
im_bgd(:,:,4)./exp(wvs);
delta_obserd_pred = img_obserd - img_pred;
mask_cirrus = delta_obserd_pred > thrd_mincir & delta_obserd_pred./img_obserd>thrd_normdiff; % cirrus: 1
clear delta_obserd_pred
mask_cirrus( im_bgd(:,:,end)==-9999) = 254; % 254: lack of enough obervations for building Cmask time series model,
mask_cirrus(img_obserd==0)=255; % filled: 255
mask_cirrus = uint8(mask_cirrus);
output_folder = fullfile(dir_cmask_cirrus_mask,name_cmask_img);
if length(dir(output_folder))==0
mkdir(output_folder);
end
enviwrite(fullfile(output_folder,[name_cmask_img,'_Cmask']),mask_cirrus,'uint8',resolu,jiUL,'bsq',ZC);
clear output_folder mask_cirrus;
if save_pcd
output_folder_pcd = fullfile(dir_cmask_pred,name_cmask_img);
if ~isfolder(output_folder_pcd)
mkdir(output_folder_pcd);
end
enviwrite(fullfile(output_folder_pcd,[name_cmask_img,'_PredCirrusBand']),img_pred,'int16',resolu,jiUL,'bsq',ZC); % pcb: predicted cirrus band
end
clear output_folder_pcd img_pred;
fprintf('Processed %dth image\n', i);
end
end