Fmask4.0代码学习记录

Fmask简介

Fmask是一种提取遥感图像中云和云阴影区域的方法，最开始由Zhu et al(2012)提出，只针对Landsat4~5、7的TM/ETM+图像，后来经过Qiu et al()的不断改进，可同时用于Landsat4~5、7、8的TM/ETM+/OLI以及Sentinel2 的MSI图像。

Fmask最新版本为4.0(Qiu et al, 2019), 代码可以从https://github.com/GERSL/Fmask下载，为MATLAB代码，同时Fmask也有Python的代码，可以安装python-fmask包。两个版本的代码应用于Sentinel-2 MSI图像时的一些问题，可以在https://forum.step.esa.int/t/sentinel-2-cloud-mask-with-fmask/4152查看。这里主要简单介绍一些MATLAB代码

Fmask代码介绍

“autoFmask.m”为程序的主入口，从代码介绍中可以看到，作者编写的Readme中可以看出，Fmask程序是通过命令行运行的，从第57行可以看到。

%% get parameters from inputs
p = inputParser;

暂且不管这些，直接跳到第99行，这里为程序计算的开始，即读入数据。

[data_meta,data_toabt,angles_view,trgt] = LoadData(path_data,sensor,InputFile,main_meta);

其中 LoadData的第一个参数 path_data，在第54行通过 path_data=pwd;赋值，pwd 在整个代码中是找不到的。

LoadData.m文件

转到 LoadData函数的定义文件，从代码注释中可以看到，这个函数即可以读取 Landsat数据，又可以读取Sentinel-2数据，考虑到Landsat和Sentinel-2波段之间的差异，LoadData需要判断不同的情况，为了抓住重点，这里主要看LoadData中的几个返回值：

1. data_meta
2. data_toabt
3. trgt
   其中 data_meta和data_toabt最为重要，它们在第49行代码定义，可以看出两个变量分别是类 ObjMeta和ObjTOABT的实例。

	%% Input data class defination
    data_meta=ObjMeta;
    data_toabt=ObjTOABT;

从变量名来看，data_meta表示图像的元数据信息，从ObjMeta类的定义可以看出，data_meta包含的是 名字、遥感器类型、分辨率、图像行列数、太阳天顶角、方位角等信息；data_toabt表示TOA反射率图像和BT亮温图像。
赋值给 data_meta和data_toabt的变量都来自与函数nd2toarbt_msi和nd2toarbt的返回结果，两个函数分别执行了对Sentinel2 MSI和Landsat图像以及元数据的读取。

Landsat原始数据包含的是各波段的 *.tif 文件和 _MTL.txt 元数据文件，Sentinel-2的原始数据包含在 .SAFE 文件夹下，元数据在 MTD_MSIL1C.xml 中，读取起来比Landsat要复杂，因此先看nd2toarbt文件中的内容。

nd2toarbt.m文件

转到nd2toarbt.m文件，从注释中可以看到，这个函数实现了Landsat图像从DN值到TOA反射率和BT亮温的转换，并且读取元数据。
函数中首先通过 lndhdrread解析了 _MTL.txt 文件，对于Landsat4~5、7和Landsat8分情况进行处理，获取由DN值转为TOA反射率和BT亮温的系数，以及其他系数：

% Outputs:
% 1) Lmax = Max radiances;  		对应 RADIANCE_MAXIMUM_BAND_*
% 2) Lmin = Min radiances; 			对应 RADIANCE_MINIMUM_BAND_*
% 3) Qcalmax = Max calibrated DNs; 	对应 QUANTIZE_CAL_MAX_BAND
% 4) Qcalmin = Min calibrated DNs; 	对应 QUANTIZE_CAL_MIN_BAND
% 5) ijdim_ref = [nrows,ncols]; % dimension of optical bands
% 6) ijdim_thm = [nrows,ncols]; % dimension of thermal band
% 7) reo_ref = 28/30; % resolution of optical bands
% 8) reo_thm = 60/120; % resolution of thermal band
% 9) ul = [upperleft_mapx upperleft_mapy];
% 10) zen = solar zenith angle (degrees); 太阳天顶角
% 11) azi = solar azimuth angle (degrees); 太阳方位角
% 12) zc = Zone Number
% 13) Lnum = 4,5,or 7; Landsat sensor number
% 14) doy = day of year (1,2,3,...,356);
% 15) pro_type=PRODUCT_TYPE (L1G: DEM bias; others: no DEM bias)

然后定义了一个 Tab_ES_Dist 表示了一年中某一天DOY对应的日地距离比例 $d$

然后对于Landsat4~5、7和Landsat8分情况进行处理，开始读取各波段图像、定标的处理过程

1. Landsat4~5、7
   图像读取利用了 imread，对于Band2，还利用了GRIDobj，如下：

    % Band2
    % Also used to be target or referenced image when resampling and reprojecting to the image extent and resolution.
    n_B2=dir(fullfile(path_top,'L*B2*'));
    if isempty(n_B2)
        n_B2=dir(fullfile(path_top,'L*b2*'));
    end
    im_B2=single(imread(fullfile(path_top,n_B2.name)));
    trgt = GRIDobj(fullfile(path_top,n_B2.name));
%     im_B2=single(trgt.Z);
    trgt.Z=[];% remove non-used values to save memory.

根据代码注释部分，trgt目的是为了在读取tiff的时候，同时获取投影等信息，进而在保存云掩模图像的时候使用，但是注释还提到GRIDobj读取方法会存在问题，实验的时候也发现报错。

% Revisions:
% Cannot use the GRIDobj to read our band, because it may result in Nan data.
% (Shi 2/26/2018)
% Added a target image based on blue band, which will be helpfull to
% resample and reproject the auxiliary data.  (Shi 2/16/2018)

定标需要将反射波段和热红外波段分别定标为TOA反射率和亮温
将反射率波段从DN值定标为TOA反射率的公式如下：
$TO{A}_{Rad}=\frac{L_{max}-L_{min}}{Qca{l}_{max}-Qca{l}_{min}}\times (DN-Qca{l}_{min})+Lmin$

$TO{A}_{Ref}=\frac{\pi \times TO{A}_{Rad}\times d^2}{E_{SUN}\times \cos ({\theta }_s)}$
其中${\theta }_s$为太阳天顶角，以弧度为单位计算，MTL文件中的SUN_ELEVATION为太阳高度角，与天顶角互余（相加等于90度）。$E_{SUN}$在代码中定义为：

    % Landsat 4, 5, and 7 solar spectral irradiance
    % see G. Chander et al. RSE 113 (2009) 893-903
    esun_L7=[1997.000, 1812.000, 1533.000, 1039.000, 230.800, -1.0, 84.90];
    esun_L5=[1983.0, 1796.0, 1536.0, 1031.0, 220.0, -1.0, 83.44];
    esun_L4=[1983.0, 1795.0, 1539.0, 1028.0, 219.8, -1.0, 83.49];

代码中还对反射率做了10000倍的放大

热红外波段定标时，先类似反射波段那样转换为辐亮度，再由辐亮度转为亮温，公式如下，温度单位为℃，代码中做了100倍放大
$BT=\frac{K_2}{\ln ((K_1/DN)+1)}-273.15$

2. Landsat8 OLI

Landsat8的图像读取过程与Landsat4~5、7相同，不同之处在于定标的公式存在差异。
T O A R e f = { R e f m a x − R e f m i n Q c a l m a x − Q c a l m i n ∗ ( D N − Q c a l m i n ) + R e f m i n } / cos ⁡ ( θ s ) TOA_{Ref} = \{ \frac{Ref_{max}-Ref_{min}}{Qcal_{max}-Qcal_{min}}*(DN-Qcal_{min})+Ref_{min} \} /\cos(\theta_s) TOARef​={Qcalmax​−Qcalmin​Refmax​−Refmin​​∗(DN−Qcalmin​)+Refmin​}/cos(θs​)

其中$Re{f}_{max}$和$Re{f}_{min}$也由 lndhdrread返回，分别对应MTL文件中的“REFLECTANCE_MAXIMUM_BAND_”和“REFLECTANCE_MINIMUM_BAND_”；
亮温计算与Landsat4、5、7方法一样。

总结数据读取过程

整个数据读入和定标处理的过程如下图所示

Fmask 代码的修改与运行

Fmask方法具体的计算过程暂时不做分析，现在把数据载入的过程理清楚了，就对程序做一下测试。
复制一份autoFmask.m的代码做修改，先把path_data赋一个文件路径，这里采用了Landsat8 的一景图像测试。

% 原来第54行代码如下
% path_data=pwd;
% 修改为：
path_data='D:\LC81250372016203LGN00';

Fmask4.0增加了对DEM、全球水体掩模的使用，作者也给出了这些数据的下载方式，可以下载之后把文件路径添加到代码中，但是作者提供的百度网盘下载地址打不开了，为了简单期间，这部分代码暂时不修改，也不使用DEM这些辅助数据。

直接运行发现在GRIDobj读取数据时报错，因此把nd2toarbt.m文件这部分代码注释掉：

    % Band2
    % Also used to be target or referenced image when resampling and reprojecting to the image extent and resolution.
    n_B2=dir(fullfile(path_top,'L*B2*'));
    if isempty(n_B2)
        n_B2=dir(fullfile(path_top,'L*b2*'));
    end
    im_B2=single(imread(fullfile(path_top,n_B2.name)));
  
% 注释掉下面三行代码，这三行代码分别在Landsat4、5、7和Landsat8图像读取的位置都出现
% 所有两个位置的代码都要注释掉
%    trgt = GRIDobj(fullfile(path_top,n_B2.name));
%    im_B2=single(trgt.Z);
%    trgt.Z=[];% remove non-used values to save memory.

但是为了保持nd2toarbt.m返回值的正常，在其代码开头加一句 trgt =[];

再返回autoFmask.m文件，最后结果保存的位置使用了 trgt，因此也需要进行修改。把 trgt的代码注释掉，利用 write_envi.m将结果保存为envi格式，保存的时候需要注意对最终得到的图像 cs_final 做一下转置transpose。

    fmask_name=[data_meta.Name,'_Fmask4'];
    fmask_outdir=fullfile(path_data,data_meta.Output)
    if ~exist(fmask_outdir)
        mkdir(fmask_outdir)
    end
    fmask_output=fullfile(fmask_outdir,[fmask_name]);
%    % output as geotiff.
%    trgt.Z=cs_final;     
%    trgt.name=fmask_name;
%    GRIDobj2geotiff(trgt,fmask_output);
    write_envi(transpose(cs_final),fmask_output);

实验结果如下图所示：
蓝色为水体，红色为云，黄色为云阴影

应用到Sentinel-2 MSIL1C数据

Sentinel-2 MSIL1C数据是经过定标的TOA反射率图像，它包含4个10m波段，6个20m波段，3个60m波段，Fmask中为了便于处理，将分辨率统一重采样到20m，最后输出的云掩模也是20m分辨率，重采样过程见nd2toarbt_msi.m中第89行for循环，利用了MATLAB的imresize函数，从10m重采样到20m时，利用box方法，60m到20m时，利用nearest方法。

        if Ratio(iB)<1 % box-agg pixel
            TOAref(:,:,iB) = imresize(dum,Ratio(iB),'box');
        elseif Ratio(iB)>1 % split pixel
            TOAref(:,:,iB) = imresize(dum,Ratio(iB),'nearest');
        elseif Ratio(iB)==1
            TOAref(:,:,iB) = dum;
        end

下载的Sentinel-2图像，解压后是一个 “.SAFE” 文件夹，下面包含多个文件夹和文件：

将Fmask用于Sentinel-2图像时，文件路径要写到 “GRANULE” 文件夹下，即

path_data='E:\temp\S2A_MSIL1C_20200222T030731_N0209_R075_T50SKE_20200222T060244.SAFE\GRANULE\L1C_T50SKE_A024382_20200222T031414'

最后得到的云掩模图像，也在 “GRANULE” 文件夹下：

.\GRANULE\L1C_T50SKE_A024382_20200222T031414\FMASK_DATA\L1C_T50SKE_A024382_20200222T031414_Fmask4.img

实验结果如下：