-
Notifications
You must be signed in to change notification settings - Fork 57
/
DetectPotentialCloudShadow.m
114 lines (105 loc) · 4.8 KB
/
DetectPotentialCloudShadow.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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% A stratied random sampling approach was applied in this version.
% Speed up by combining the correction for NIR and SWIR bands together.
% by Shi.
% at 23,Feb.,2017
%
% No topo correction when no DEMs (Shi 9/13/2017)
% fixed the bug of the wrong selection of samples (Shi 2/24/2017)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function masker_shadow = DetectPotentialCloudShadow(data_meta, data_nir,data_swir,data_clear_land,masker_observation,...
slope_data,aspect_data)
% History:
% 1. A total of samples are revised as 40,000 from 50,000. (8. March, 2018)
percent_low=0.175;
sun_zenith_deg=data_meta.Zen;
sun_azimuth_deg=data_meta.Azi;
dim=data_meta.Dim;
% get potential layer
%% SCS+C for topo correct...
if ~isempty(slope_data)&&~isempty(aspect_data)
[data_nir,data_swir]=getDataTopoCorrected(data_nir,data_swir,data_clear_land,sun_zenith_deg,sun_azimuth_deg, slope_data,aspect_data,dim,data_meta.Resolution(1));
end
clear sun_zenith_deg sun_azimuth_deg slope_data aspect_data;
%% band 4 flood fill
% estimating background (land) Band 4 & 5 ref
backg_b4=prctile(data_nir(data_clear_land),100*percent_low);
backg_b5=prctile(data_swir(data_clear_land),100*percent_low);
data_nir(masker_observation==0)=backg_b4;
data_nir(isnan(data_nir(:)))=backg_b4;
clear backg_b4;
% fill in regional minimum Band 4 ref
data_nir_filled=imfill(data_nir);
data_nir_dif=data_nir_filled-data_nir;
clear data_nir data_nir_filled;
%% band 5 flood fill
% estimating background (land) Band 5 ref
data_swir(masker_observation==0)=backg_b5;
data_swir(isnan(data_swir(:)))=backg_b5;
clear backg_b5;
clear data_clear_land percent_low;
% fill in regional minimum Band 5 ref
data_swir_filled=imfill(data_swir);
data_swir_dif=data_swir_filled-data_swir;
clear data_swir data_swir_filled;
%% compute shadow probability
shadow_prob=min(data_nir_dif,data_swir_dif);
% get potential cloud shadow
masker_shadow=zeros(dim,'uint8'); % shadow mask
masker_shadow(shadow_prob>200)=1;
masker_shadow(masker_observation==0)=255;
end
% CSC+C stratied on cos i with 0.1 increasement. a total of 50,000 pixels.
function [data_nir,data_swir] = getDataTopoCorrected(data_nir_ori,data_swir_ori,index_exclude_cloud_water,sun_zenith_deg,sun_azimuth_deg, slope_data,aspect_data,dim,resl )
% History:
% 1. Create this function. (1. January, 2017)
% 2. A total of samples are revised as 40,000 from 50,000. (8. March, 2018)
% 3. When c is calculated as Nan, this function will not make the topo
% correction. (8. March, 2018)
sun_zenith_rad=deg2rad(sun_zenith_deg);
sun_zenith_cos=cos(sun_zenith_rad);
sun_zenith_sin=sin(sun_zenith_rad);
clear sun_zenith_deg sun_zenith_rad sun_zenith_rad;
cos_sita=sun_zenith_cos.*cos(deg2rad(slope_data))+sun_zenith_sin.*sin(deg2rad(slope_data)).*cos(deg2rad(sun_azimuth_deg-aspect_data));
clear aspect_data;
cos_sita_exclude_cloud=cos_sita(index_exclude_cloud_water);
% random stratied sampling
cos_sita_min=min(cos_sita_exclude_cloud);
cos_sita_max=max(cos_sita_exclude_cloud);
% total_sample=50000;
total_sample=40000;
cos_sita_interval=0.1;
samples_ids= stratiedSampleHanlder(cos_sita_exclude_cloud,cos_sita_min,cos_sita_max,dim,total_sample,cos_sita_interval,0,resl);
cos_sita_samples=cos_sita_exclude_cloud(samples_ids);
clear cos_sita_exclude_cloud cos_sita_min cos_sita_max total_sample cos_sita_interval;
% for NIR
data_nir_ori_tmp=data_nir_ori(index_exclude_cloud_water);
data_samples_nir=data_nir_ori_tmp(samples_ids);
clear data_nir_ori_tmp;
c_fitted=polyfit(double(cos_sita_samples),double(data_samples_nir),1);
% figure;plot(double(cos_sita_samples),double(data_samples_nir),'r.');
c=c_fitted(1,2)/c_fitted(1,1);
clear c_fitted;
if isnan(c)
data_nir=data_nir_ori;
clear data_nir_ori;
else
data_nir=double(data_nir_ori).*(cos(deg2rad(slope_data)).*sun_zenith_cos+c)./(cos_sita+c);
clear data_nir_ori;
end
% for SWIR
data_swir_ori_tmp=data_swir_ori(index_exclude_cloud_water);
data_samples_swir=data_swir_ori_tmp(samples_ids);
clear data_swir_ori_tmp;
c_fitted=polyfit(double(cos_sita_samples),double(data_samples_swir),1);
c=c_fitted(1,2)/c_fitted(1,1);
clear c_fitted samples_ids;
if isnan(c)
data_swir=data_swir_ori;
clear data_swir_ori;
else
data_swir=double(data_swir_ori).*(cos(deg2rad(slope_data)).*sun_zenith_cos+c)./(cos_sita+c);
clear data_swir_ori;
end
end