Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@qsly09
qsly09 committed May 16, 2019
1 parent 7032ec4 commit e94d56e
Show file tree
Hide file tree
Showing 4 changed files with 419 additions and 0 deletions.
141 changes: 141 additions & 0 deletions BRDFAdjust.m
View file
@@ -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

86 changes: 86 additions & 0 deletions TopoCorrectIC.m
View file
@@ -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

19 changes: 19 additions & 0 deletions TopoCorrectSCS.m
View file
@@ -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

0 comments on commit e94d56e

Please sign in to comment.