遥感入门 | 使用matlab读取Landsat8

matlab中文翻译过来叫做“矩阵实验室”，众所周知，遥感图像是图像， 图像是矩阵，所以matlab处理遥感图像很在行。

matlab一直是科研圈子里验证算法或者研究算法的利器，以前听师兄说，matlab除了不会生孩子之外，什么都能做。

下面是一段matlablandsat 8 原始影像的代码.

代码大概的流程是：

1.分别读取第一波段到第七波段，得到七个矩阵，就是七个图像。

2.读取元数据xml，各自波段进行辐射定标，得到辐亮度数据。

3.进行系统级大气校正，得到TOA大气顶部反射率。（不是SR地表反射率）

4.将原始DN值量化为1，因为原始的Landsat 8 是16bit。进行原始影像真彩色合成， 并可视化。

5.利用TOA进行水体指数NDWI计算， 并可视化。

clc;clear all;close all;

for i=1:7
    %设置待处理读取文件的路径
    loadpath = 'D:\DDD\a2.28\input\landsat 8\LC81220442017296LGN00';
    
    txt='*_MTL.txt';    
    file_txt = dir([loadpath,'\',txt]).name;
    targetfile = dir([loadpath,'\','*B',num2str(i),'.TIF']).name;
    
    %读取头文件的辐射定标的信息
    MDATA = parseLandSat8MetaData(file_txt).L1_METADATA_FILE(1,1).RADIOMETRIC_RESCALING(1,1);
    % 增益
    RADIANCE_MULT = ['RADIANCE_MULT_BAND_' num2str(i) '=MDATA.RADIANCE_MULT_BAND_' num2str(i) ';' ];
    eval(RADIANCE_MULT);
    % 偏移
    RADIANCE_ADD = ['RADIANCE_ADD_BAND_' num2str(i) '=MDATA.RADIANCE_ADD_BAND_' num2str(i) ';' ];
    eval(RADIANCE_ADD);
    
    %原始影像
    Image=  ['B' num2str(i) '=double(imread(targetfile));'];
    eval(Image);
    
    %过滤无效值
    n=['B' num2str(i) '(B' num2str(i) '<0)=nan;'];
    eval(n);
    
    %辐射定标
    Rad= ['B' num2str(i) '_rad' '=B' num2str(i) '*RADIANCE_MULT_BAND_' num2str(i) '+RADIANCE_ADD_BAND_' num2str(i) ';'];
    eval(Rad);
    %过滤无效值
    n=['B' num2str(i) '_rad' '(B' num2str(i) '_rad' '<0)=nan;'];
    eval(n);
    
    %% 大气顶部反射率
    %漏=pi*辐亮度*日地天文距离^2/(E*cos(theta))
    %漏是大气顶部反射率；E是太阳辐照度；theta是太阳天顶角；
    
    %如果存在REFLECTANCE_MULT_BAND 反射增益、偏移参数，
    %则可以使用公式   漏=（M*Q+A）/sin(b)
    %漏是大气顶部反射率；M是反射增益参数；Q是DN值；A是反射偏移参数；b是太阳高度角

    MDATA1 = parseLandSat8MetaData(file_txt);
    
    %% 计算儒略日时间
    date=datetime(MDATA1.L1_METADATA_FILE(1,1).PRODUCT_METADATA(1:1).DATE_ACQUIRED);
    ds=Julian_Day(date.Year,date.Month,date.Day);
    
    %%
    %太阳高度角
    b=90-MDATA1.L1_METADATA_FILE(1,1).IMAGE_ATTRIBUTES.SUN_ELEVATION;
        % 增益
    REFLECTANCE_MULT = ['REFLECTANCE_MULT_BAND_' num2str(i) '=MDATA.REFLECTANCE_MULT_BAND_' num2str(i) ';' ];
    eval( REFLECTANCE_MULT);
    % 偏移
     REFLECTANCE_ADD = ['REFLECTANCE_ADD_BAND_' num2str(i) '=MDATA.REFLECTANCE_ADD_BAND_' num2str(i) ';' ];
    eval(REFLECTANCE_ADD);

    %
        %大气校正（大气顶层反射率）
       
    Ref= ['B' num2str(i) '_ref' '=(B' num2str(i) '*REFLECTANCE_MULT_BAND_' num2str(i) '+REFLECTANCE_ADD_BAND_' num2str(i) ')/sin(b);'];
    eval(Ref);
%     %过滤无效值
    n2=['B' num2str(i) '_ref' '(B' num2str(i) '_ref<0)=nan;'];
    eval(n2);
    
end
%将原始DN值量化为1，因为Landsat 8 是16bit








B2=B2/(2^16);
B3=B3/(2^16);
B4=B4/(2^16);
B5=B5/2^16;
% b1=b/255;

%原始影像真彩色合成   波段4 3 2
G(:,:,1)=B4_ref(:,:,1);
G(:,:,2)=B3_ref(:,:,1);
G(:,:,3)=B2_ref(:,:,1);

%原始影像真彩色合成   波段5 4 3
L(:,:,1)=B5;
L(:,:,2)=B4;
L(:,:,3)=B3;
figure,
subplot(121),imshow(G);
subplot(122),imshow(L);

%水体指数计算
ndwi=(B2_ref-B5_ref)./(B2_ref+B5_ref);

%画图
figure;
hold on;
imshow(ndwi,[-1 1]);
%添加颜色棒
colorbar;
colormap(jet(256));

以上代码是早些年写的

小结

1.matlab对传统算法复现很快，而且便于人类去理解计算机的算法运行逻辑，让人加深对某个算法的理解，从而进一步去创新。

2.但是matlab不开源，且不免费，也不利于后续的工程化。