diff --git a/BRDFAdjust.m b/BRDFAdjust.m new file mode 100644 index 0000000..71d1fc4 --- /dev/null +++ b/BRDFAdjust.m @@ -0,0 +1,141 @@ +function ref_norm = BRDFAdjust( ref_ori,band_name,... + solar_zenith_angle,view_zenith_angle,solar_azimuth_angle,view_azimuth_angle,... + centre_lat, solar_zenith_angle_norm) +% BRDFADJUST adjusts Prepare U.S. Landsat Analysis Ready Data (ARD) into CCDC +%format. +% OLI 5,110,31,72 +% MSI 11,106,29,68 +%-------------------------------------------------------------------------- +% Summary: +% This function is to adjust the BRDFs of Landsat and Sentinel images to +% niar observations. +%-------------------------------------------------------------------------- +% History of revisions: +% modify this function. by Rong Shang 6/27/2018 +% create this function. by Zhanmang Liao 5/13/2017 +% lamda=3;%Ñband number +%========================================================================== + solar_zenith_angle = double(solar_zenith_angle); + view_zenith_angle = double(view_zenith_angle); + solar_azimuth_angle = double(solar_azimuth_angle); + view_azimuth_angle = double(view_azimuth_angle); + + solar_zenith_angle = deg2rad(solar_zenith_angle); + view_zenith_angle = deg2rad(view_zenith_angle); + relative_azimuth_angle = deg2rad(view_azimuth_angle-solar_azimuth_angle); + + % if have centre's latitude + if exist('centre_lat','var')&&~isempty(centre_lat) + solar_zenith_angle_norm_centre = 31.0076 - 0.1272*centre_lat + 0.01187*centre_lat^2 + ... + 2.4*10^(-5)*centre_lat^3 - 9.48*10^(-7)*centre_lat^4 - ... + 1.95*10^(-9)*centre_lat^5 + 6.15*10^(-11)*centre_lat^6; + solar_zenith_angle_norm_centre = deg2rad(solar_zenith_angle_norm_centre); + end + + % if denfine solar angle + if exist('solar_zenith_angle_norm','var')&&~isempty(solar_zenith_angle_norm) + solar_zenith_angle_norm_custom = deg2rad(solar_zenith_angle_norm); + end + + if exist('solar_zenith_angle_norm_custom','var')&&~isempty(solar_zenith_angle_norm_custom) + % do not change the solar zenith angle if no any inputs because the + % Landsat solar zenith angles are very similar. + solar_zenith_angle_norm = solar_zenith_angle_norm_custom; + else + if exist('solar_zenith_angle_norm_centre','var')&&~isempty(solar_zenith_angle_norm_centre) + solar_zenith_angle_norm = solar_zenith_angle_norm_centre; + else + % do not change solar angle if no self-defined one + solar_zenith_angle_norm = solar_zenith_angle; + end + end + + view_zenith_angle_norm = deg2rad(0); + relative_azimuth_angle_norm = deg2rad(180); + + % band name should be converted to band 1 2 3 4 5 and 7 for Landsat 7 data. + switch band_name + case 'Blue' + lamda = 1; % band 1 + case 'Green' + lamda = 2; % band 2 + case 'Red' + lamda = 3; % band 3 + case 'NIR' + lamda = 4; % band 4 + case 'SWIR1' + lamda = 5; % band 5 + case 'SWIR2' + lamda = 6; % band7 + otherwise + ref_norm = []; + warning('Please input right band name. (See BRDFAdjust)'); + return; + end + + %%======== BRDF parameters ======== + % band 1 2 3 4 5 7 + f_iso= [0.0774, 0.1306, 0.1690, 0.3093, 0.3430, 0.2658]; + f_vol= [0.0372, 0.0580, 0.0574, 0.1535, 0.1154, 0.0639]; + f_geo= [0.0079, 0.0178, 0.0227, 0.0330, 0.0453, 0.0387]; + tmp = f_iso(lamda); clear f_iso; + f_iso = tmp; + tmp = f_geo(lamda); clear f_geo; + f_geo = tmp; + tmp = f_vol(lamda); clear f_vol; + f_vol = tmp; + %clear tmp; + + %======== calculate the kernel ======== + [k_vol_norm,k_geo_norm]=kernel(solar_zenith_angle_norm, view_zenith_angle_norm, relative_azimuth_angle_norm); + % clear solar_zenith_angle_norm view_zenith_angle_norm relative_azimuth_angle_norm; + [k_vol_sensor,k_geo_sensor]=kernel(solar_zenith_angle, view_zenith_angle, relative_azimuth_angle); + % clear solar_zenith_angle view_zenith_angle relative_azimuth_angle; + + %======== calculate correcting parameter ======== + P1=f_iso+f_geo*k_geo_norm+f_vol*k_vol_norm; +% clear k_vol_norm; + P2=f_iso+f_geo*k_geo_sensor+f_vol*k_vol_sensor; +% clear k_vol_sensor; + C_lamda=P1./P2; + + ref_norm = int16(C_lamda.*double(ref_ori)); %corrected reflectance + +end + +function [k_vol,k_geo] = kernel(theta_i,theta_v,azimuth_angle) +b=1;r=1;h=2; + +%================ calculate k_vol ================ +cos_g = cos(theta_i).*cos(theta_v)+sin(theta_i).*sin(theta_v).*cos(azimuth_angle); +g = acos(max(-1,min(cos_g,1))); +clear cos_g; + +k_vol = ((0.5*pi-g).*cos(g)+sin(g))./(cos(theta_i)+cos(theta_v))-pi/4; +clear g; + + +%================ calculate k_geo ================ +theta_i1 = atan(max(b/r*tan(theta_i),0)); +theta_v1 = atan(max(b/r*tan(theta_v),0)); +clear theta_i theta_v r; + +g_1 = cos(theta_i1).*cos(theta_v1) + sin(theta_i1).*sin(theta_v1).*cos(azimuth_angle); +g_1 = acos( max(-1,min(1,g_1)) ); + +D = tan(theta_i1).^2+tan(theta_v1).^2-2*tan(theta_i1).*tan(theta_v1).*cos(azimuth_angle); +D = sqrt(max(D,0)); + +cos_t = h/b*( sqrt(D.^2+(tan(theta_i1).*tan(theta_v1).*sin(azimuth_angle)).^2) )./(sec(theta_i1)+sec(theta_v1)); +cos_t = max(-1,min(1,cos_t)); +clear D h b azimuth_angle; + +t = acos(cos_t); +clear cos_t; + +O = 1/pi*(t-sin(t).*cos(t)).*(sec(theta_i1)+sec(theta_v1)); +O = max(0,O); + +k_geo = O - sec(theta_i1)-sec(theta_v1) + 1/2*(1+cos(g_1)).*sec(theta_i1).*sec(theta_v1); +end + diff --git a/TopoCorrectIC.m b/TopoCorrectIC.m new file mode 100644 index 0000000..2585c01 --- /dev/null +++ b/TopoCorrectIC.m @@ -0,0 +1,86 @@ +function data_nir = TopoCorrectIC(surf_band,index_exclude_cloud,sun_zenith,sun_azimuth, slope,aspect,ndvi) + data_nir = surf_band; + jidim = size(surf_band); + + step = 50; + buffer = 25; +% for i = 5000 +% for j = 5000 + for i = 1:step:jidim(1) + for j = 1:step:jidim(2) + % extract window images + data_tmp = double(data_nir(max(i-buffer,1):min(i+step+buffer-1,jidim(1)),max(j-buffer,1):min(j+step+buffer-1,jidim(2)))); + data_nir_ori = surf_band(max(i-buffer,1):min(i+step+buffer-1,jidim(1)),max(j-buffer,1):min(j+step+buffer-1,jidim(2))); + index_exclude_cloud_water = index_exclude_cloud(max(i-buffer,1):min(i+step+buffer-1,jidim(1)),max(j-buffer,1):min(j+step+buffer-1,jidim(2))); + sun_zenith_deg = sun_zenith(max(i-buffer,1):min(i+step+buffer-1,jidim(1)),max(j-buffer,1):min(j+step+buffer-1,jidim(2))); + sun_azimuth_deg = sun_azimuth(max(i-buffer,1):min(i+step+buffer-1,jidim(1)),max(j-buffer,1):min(j+step+buffer-1,jidim(2))); + slope_data = slope(max(i-buffer,1):min(i+step+buffer-1,jidim(1)),max(j-buffer,1):min(j+step+buffer-1,jidim(2))); + aspect_data = aspect(max(i-buffer,1):min(i+step+buffer-1,jidim(1)),max(j-buffer,1):min(j+step+buffer-1,jidim(2))); + ndvi_ind = ndvi(max(i-buffer,1):min(i+step+buffer-1,jidim(1)),max(j-buffer,1):min(j+step+buffer-1,jidim(2))); + ndvi_gt = (ndvi_ind>=0.5&ndvi_ind<=1); + ndvi_lt = (ndvi_ind<0.5&ndvi_ind>=0); + % IC topographic correction + tmp_ndvi = TopoCorrectICWindow(data_nir_ori,index_exclude_cloud_water,sun_zenith_deg,sun_azimuth_deg, slope_data,aspect_data,ndvi_gt); + data_tmp(ndvi_gt==1) = double(tmp_ndvi); + tmp_ndvi = TopoCorrectICWindow(data_nir_ori,index_exclude_cloud_water,sun_zenith_deg,sun_azimuth_deg, slope_data,aspect_data,ndvi_lt); + data_tmp(ndvi_lt==1) = double(tmp_ndvi); + data_nir(max(i-buffer,1):min(i+step+buffer-1,jidim(1)),max(j-buffer,1):min(j+step+buffer-1,jidim(2))) = data_tmp; + end + end +end + +% CSC+C stratied on cos i with 0.1 increasement. a total of 50,000 pixels. +function [band_data] = TopoCorrectICWindow(band_ori,index_exclude_cloud_water,sun_zenith_deg,sun_azimuth_deg, slope_data,aspect_data,ndvi_ind) +% 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) + +% data_nir = data_nir_ori; +% ori = data_nir_ori; + [ndvi_ind1] = find(ndvi_ind==true); + if isempty(ndvi_ind1) + band_data = []; + return + end + ndvi_ind(1:2:end,1:2:end) = false; + [ndvi_ind2] = find(ndvi_ind==true); + loc_ind = ismember(ndvi_ind1,ndvi_ind2)==1; + + + % extract the corresponding ndvi>0.5 or ndvi<=0.5 pixels + band_ori = band_ori(ndvi_ind1); + % make sure the output image size is the same with input + band_data = band_ori; + index_exclude_cloud_water = index_exclude_cloud_water(ndvi_ind1); + sun_zenith_deg = sun_zenith_deg(ndvi_ind1); + sun_azimuth_deg = sun_azimuth_deg(ndvi_ind1); + slope_data = slope_data(ndvi_ind1); + aspect_data = aspect_data(ndvi_ind1); + + sun_zenith_rad=deg2rad(double(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(single(sun_azimuth_deg)-aspect_data)); + + + pix_ind = ((double(cos_sita)-double(sun_zenith_cos))<=-0.05)|((double(cos_sita)-double(sun_zenith_cos))>=0.05); + + clear aspect_data; + cos_sita_exclude_cloud=cos_sita(index_exclude_cloud_water&cos_sita>=0&loc_ind&pix_ind); + data_nir_ori_tmp=band_ori(index_exclude_cloud_water&cos_sita>=0&loc_ind&pix_ind); + + if length(cos_sita_exclude_cloud)<=10 + band_data = band_ori; + else + + c_fitted_1=polyfit(double(cos_sita_exclude_cloud),double(data_nir_ori_tmp),1); + band_data(pix_ind)=double(band_ori(pix_ind))-c_fitted_1(1,1)*(cos_sita(pix_ind)-sun_zenith_cos(pix_ind)); + + end + +end + diff --git a/TopoCorrectSCS.m b/TopoCorrectSCS.m new file mode 100644 index 0000000..72ce17f --- /dev/null +++ b/TopoCorrectSCS.m @@ -0,0 +1,19 @@ +function band_corrected = TopoCorrectSCS(band_ori,sun_zenith_deg,sun_azimuth_deg, slope_data,aspect_data) +%TOPOCORRECT makes a topographic correction. +% Detailed explanation goes here + band_corrected = band_ori; + sun_zenith_rad=deg2rad(double(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(single(sun_azimuth_deg)-aspect_data)); + clear aspect_data; + + % when -0.05 < cos_sita - sun_zenith_cos < 0.05, DO NOT make + % correction. Ref. Tan et al. RSE (2013) + cor_mask = ((double(cos_sita)-double(sun_zenith_cos))<=-0.05)|... + ((double(cos_sita)-double(sun_zenith_cos))>=0.05); + + band_corrected(cor_mask)=double(band_ori(cor_mask)).*(cos(deg2rad(slope_data(cor_mask))).*sun_zenith_cos(cor_mask))./(cos_sita(cor_mask)); +end + diff --git a/TopoCorrectSCSplusC.m b/TopoCorrectSCSplusC.m new file mode 100644 index 0000000..c8dc716 --- /dev/null +++ b/TopoCorrectSCSplusC.m @@ -0,0 +1,173 @@ +function [sr_cor,c] = TopoCorrectSCSplusC(sr_ori,clr_mask,sun_zenith_deg,sun_azimuth_deg, slope_data,aspect_data) + +% TOPOCORRECTCSCPLUSC SCS + C stratied on cos i with 0.1 increasement. +% a total of 40,000 pixels. +% +% Ref. Qiu, Shi, et al. "Improving Fmask cloud and cloud shadow detection +% in mountainous area for Landsats 48 images." Remote Sensing of +% Environment 199 (2017): 107-119. +% +% +% Inputs: +% sr_ori: orginal surface reflectance +% +% clr_mask: clear sky mask; 1 is clear; 0 is non-clear, such as cloud and +% cloud shadow. +% +% sun_zenith_deg: the zenith angle of sun; unit: decimal degree +% +% sun_azimuth_deg: the azimuth angle of sun; unit: decimal degree +% +% slope: slope derived from DEM; unit: decimal degree +% +% aspect: aspect derived from DEM; unit: decimal degree +% +% +% History: +% 1. Create this function. (1. January, 2017 by Shi Qiu) +% 2. A total of samples are revised as 40,000 from 50,000. (8. March, 2018 by Shi Qiu) +% 3. When c is calculated as Nan, this function will not make the topo +% correction. (8. March, 2018 by Shi Qiu) +% 4. when -0.05 < cos_sita - sun_zenith_cos < 0.05, DO NOT make +% correction. (27. July, 2018 by Yukun Lin) + + sr_cor = sr_ori; + + sun_zenith_rad=deg2rad(double(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(single(sun_azimuth_deg)-aspect_data)); + clear aspect_data; + cos_sita_exclude_cloud=cos_sita(clr_mask); + % random stratied sampling + cos_sita_min=min(cos_sita_exclude_cloud); + cos_sita_max=max(cos_sita_exclude_cloud); + + % when -0.05 < cos_sita - sun_zenith_cos < 0.05, DO NOT make + % correction. Ref. Tan et al. RSE (2013) + cor_mask = ((double(cos_sita)-double(sun_zenith_cos))<=-0.05)|... + ((double(cos_sita)-double(sun_zenith_cos))>=0.05); + +% % % do not use the pixels what we do not need to correct +% % cos_sita_exclude_cloud = cos_sita_exclude_cloud & cor_mask; + +% total_sample=50000; + total_sample=40000; % enough derived from Fmask 4.0 (Qiu at al., 2018) + cos_sita_interval=0.1; + samples_ids= stratiedSampleHanlder(cos_sita_exclude_cloud,cos_sita_min,cos_sita_max,total_sample,cos_sita_interval); + + 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; + +% sr_ori_tmp=sr_ori(clr_mask&cor_mask); + data_samples_sr=sr_ori(samples_ids); + clear sr_ori_tmp; + c_fitted=polyfit(double(cos_sita_samples),double(data_samples_sr),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) + sr_cor=sr_ori; + clear sr_ori; + else + sr_cor(cor_mask)=double(sr_ori(cor_mask)).*(cos(deg2rad(slope_data(cor_mask))).*sun_zenith_cos(cor_mask)+c)./(cos_sita(cor_mask)+c); + clear sr_ori; + end +end + +function samples_ids=stratiedSampleHanlder(data_dem_clear,dem_b,dem_t,total_sample,ele_strata) +%STRATIEDSAMPLEHANLDER Further select clear sky (land) pixels to normalize +%BT. +% +% Syntax +% +% samples_ids= +% stratiedSampleHanlder(data_dem_clear,dem_b,dem_t,dim,total_sample,ele_strata,distance) +% +% Description +% +% Stratied sampling method is used by Qiu et al., (2017). +% History: +% 1. Create this function. (1. January, 2017) +% 2. This stratied sampling method sometimes results in no enough +% samples. (8. March, 2018) +% +% Input arguments +% +% data_dem_clear Clear sky pixels' DEM. +% dem_b Basic elevation. +% dem_t Highest elevation. +% dim Dim for data. +% total_sample All clear sky samples (40,000). +% ele_strata Stratied sampling along elevation (300 meters). +% distance Minmum distance among different samples(450 meters). +% Distance rule will be removed if distance = 0. +% resolution Spatial resolution (Landsat 30 meters; Sentinel-2 20 meters). +% Output arguments +% +% Tempd Nomalized Temperature (BT). +% +% +% Author: Shi Qiu (shi.qiu@ttu.edu) +% Date: 8. March, 2018 + + +% % num_clear_sky_pixels=numel(data_dem_clear); +% % % no enough clear-sky pixels afther removing the pixels out of the upper +% % % and lower levels.more than 1/4 total samples (10,000) can be used for +% % % estimating lapse rate and c in topo correction. +% % if num_clear_sky_pixels=i_dem&data_dem_clear0 +% num_strata=num_strata+1; + strata_avail=[strata_avail;1]; + else + strata_avail=[strata_avail;0]; + end + clear dem_clear_index_tmp; + end + % equal samples in each stratum + num_sam_each_strata=round(total_sample/sum(strata_avail)); + clear total_sample; + samples_ids=[];% to restore samples selected. + % loop each strata and select samples + loop_i=1; + for i_dem=dem_b:ele_strata:dem_t + if strata_avail(loop_i) + % find all clear-sky pixels in subgroup. + samples_ids_tmp=find(data_dem_clear>=i_dem&data_dem_clearnum_max_tmp + num_tmp=num_max_tmp; + end + clear num_max_tmp; + samples_ids_rand_tmp=samples_ids_rand(1:num_tmp); + clear num_tmp; + % store data + samples_ids=[samples_ids;samples_ids_rand_tmp]; + clear samples_ids_rand samples_ids_rand_tmp; + end + loop_i=loop_i+1; + end + clear samp_dis num_sam_each_strata; + +% % % no enough clear-sky pixels afther removing the pixels out of the upper +% % % and lower levels.more than 1/4 total samples (10,000) can be used for +% % % estimating lapse rate and c in topo correction. +% % if numel(samples_ids)