陆地卫星辐射校准 python 陆地卫星 5 如何进行辐射校准
前言
一、辐射定标
辐射定标：将记录的原始DN值转换为辐射值，也就是将图像的亮度灰度值转换为绝对辐射亮度。
如下便是通过imtool函数查看Landsat_5 TM band5的原始亮度灰度值，其整体区间在0~255左右，数据格式为uint8类型。
那该如何将像元亮度灰度值转换为绝对辐射亮度呢？
辐射标定的公式：
从上面的公式不难看出，DN是已知的，但是却缺少k和c，在遥感领域通常将其分别称为Gain 和 bais，这两个值其实在头文件中可以找到RANDIANCE_MULT~ 以及 RADIANCE_ADD~ 分别对于的增益和偏移参数。
知道上面的值，我们在MATLAB中实现以下，代码如下
%辐射标定
%读取头文件
imtool(image{5});
head_data=readtable("LT05_L1TP_127039_19860802_20170221_01_T1_MTL.txt");
gain=table2array(head_data(145:151,'L1_METADATA_FILE'));
bias=table2array(head_data(152:158,'L1_METADATA_FILE'));
fs_image={};
for ii=1:length(image)-1
for jj=1:length(gain)
fs_image{ii,1}=double(image{ii}).*gain(jj)+bias(jj);
end
end二、计算大气层表观反射率
至于大气层表观反射率的概念，在此就不叙述了，大家可参考别的博客，下面直接上公式：
式中：
L表示辐射亮度值；
d表示日地距离
θ表示太阳高度角
ESUM表示太阳均辐射度，可以在官网查询
中间两个都是已知值，可以在头文件中找到：
而ESUM需要去权威网站查询，如下所示：
将以上值应用于公式可以得到：
%计算大气层表观反射率
%读取日地天文单位距离以及太阳高度角
sun_elevation=table2array(head_data(61,"L1_METADATA_FILE"));
earth_sun_distance=table2array(head_data(62,"L1_METADATA_FILE"));
%太阳辐照度均值
ESUN=[1997,1812, 1533, 1039 ,230.8 , 0 ,84.90];
%计算大气层顶表观反射率
Reflectance_value={};
for oc=1:length(image)-1
Reflectance_value{oc}=(double(fs_image{oc}).*(earth_sun_distance ...
.^2).*pi)./(ESUN(oc).*sind(sun_elevation));
end三、查看定标结果
我们在图像中随机选取两个点，注意该图像只作为底图，是为了获取点的坐标而设立的参考图，实际点的属性值与其无关。可以看见，其两个点的辐射率以及大气表观反射率
程序清单>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
msgbox函数：弹出提示框，自定义提示内容
ginput函数：从当前图框中选择点，并返回点相对于图框的坐标值
diag函数：提取矩阵的对角线元素
cell2mat函数：将元胞数组转换为矩阵
end>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
msgbox('温馨提示：请查看弹出框图的左上角标题，按标题操作！！！！');
imshow(hend.*3);
[xx,yy]=ginput;
hold on
scatter(xx,yy,1000,'r','+');
len_data=length(xx);
if len_data>2
msgbox('输入的点大于2，请重新运行！！！')
end
band_info={};
for jr =1:length(fs_image)
map_data1=double(fs_image{jr});
cc1=map_data1(int16(yy),int16(xx));
band_info1{jr}=diag(cc1);
map_data3=double(Reflectance_value{jr});
cc3=map_data3(int16(yy),int16(xx));
band_info3{jr}=diag(cc3);
end
feshelv_matrix1=cell2mat(band_info1);
feshelv_matrix3=cell2mat(band_info3);
feshelv_matrix3(feshelv_matrix3==Inf)=0;
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
for ta=1:2
x=[1,2,3,4,5,6,7];
y=feshelv_matrix1(ta,:);
figure;
plot(x,y,'Color','r');
xlabel('bands');
ylabel('Data Value');
name=strcat('the Radiance--value with gain and bias of--','point--',num2str(ta));
title(name);
axis([0 7 0 10]);
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
y3=feshelv_matrix3(ta,:);
figure;
plot(x,y3,'Color','r');
xlabel('bands');
ylabel('Data Value');
name=strcat('the Reflectance--value of --','point--',num2str(ta));
title(name);
axis([0 7 0 0.2]);
end