MODIS 二级气溶胶产品指北（with Python）

MODIS 二级气溶胶产品指北（with Python）

笔者初次使用 MODIS 二级气溶胶产品时，一下就被密密麻麻一堆变量搞懵了：很多变量名字里带个 Optical_Depth，这我能猜到，就是气溶胶光学厚度，但各种 Corrected、Effective、Best、Average、Small、Large 的前后缀鬼知道是什么。看过的论文基本不说具体用的哪个变量，各种教程也不会告诉你这些亲戚间的差异，似乎这件事一点也不重要。本着 know your data 的心态，我在翻阅了 MODIS 的几个官网后总算从反演的原理中稍微体会到了这些前后缀的意义。现将学习经验总结归纳如下，希望能帮到和我一样疑惑的小伙伴。同时本文还会提供简单的 Python 示例代码。

如果嫌正文太啰嗦，可以直接跳到文末的总结部分，那里直接给出了各个变量的使用建议。

MODIS 全称 MODerate resolution Imaging Spectroradiometers（中分辨率成像光谱仪），基本信息为：

• 搭载平台：Terra（2000 年至今）和 Aqua（2002 年至今）。

• 轨道：太阳同步，上午 10:30 降轨通过赤道（Terra）或下午 1:30 升轨通过赤道（Aqua）。

• 波段：0.41 ~ 14.5 μm的 36 个波段（即可见光到远红外）。

• 刈幅：2330 km。

• 空间分辨率：250 m（bands 1 - 2）、500 m（bands 3 - 7）、1000 m（bands 8 - 36）。

上图为大气层顶（绿线）和海平面（红线）处的太阳入射光谱，以及 MODIS 的各个波段（橙线）。

MODIS 的优势可以概括为两点：

• 多波段的辐射测量能揭示关于陆面、海洋和大气的丰富信息。
• 高空间分辨率、太阳同步轨道和宽广的刈幅能覆盖宽广的时空范围

气溶胶光学厚度（Aerosol Optical Depth，AOD）反映到达地表的太阳辐射受气溶胶的衰减程度，在无云的晴空下，太阳光度计可以通过直接测量太阳方向的辐射来计算 AOD。星载的 MODIS 与之不同，如上图所示，其镜头朝向地表，测量的是经地表反射和气溶胶散射后的太阳辐射。利用地表和气溶胶在不同波段的辐射性质差异，有机会从测量结果中分离出气溶胶的信号，进而推测气溶胶的物理性质，这便是 MODIS 反演气溶胶的思路。当然现实是全球的地表类型和气溶胶组成不尽相同，为了简化问题，需要根据观测资料提前对地表性质和气溶胶组成做出一些假设。目前 MODIS 反演气溶胶的业务算法有两种：

• 暗目标算法（Dark Target，DT）：洋面和植被覆盖的陆面在可见光波段因为反射率较低显得偏暗，被称作“暗目标”，而气溶胶在这一背景上显得偏亮，由此可以忽略或估算地表的信号，进而分离出气溶胶的信号。因为陆面的地表类型更为复杂，所以 DT 算法在实现上又细分为针对洋面的 DT-ocean 和针对陆面的 DT-land。
• 深蓝算法（Deep Blue，DB）：干旱和沙漠类型的地表即便在可见光波段也显得很亮，DT 算法成立的假设失效。而在近紫外的“深蓝”波段这些地表又变得偏暗，所以借助这一波段可以扩展反演的覆盖范围。目前深蓝算法仅应用于陆面。

MODIS 的二级气溶胶产品名为 MOD04_L2 和 MYD04_L2，其中 MOD 和 MYD 分别代表 Terra 和 Aqua 卫星。每个文件对应一段 5 分钟长度的 MODIS 轨道片段（称作 granule），星下点的空间分辨率为 10 x 10 km，共含 203 x 135 个像元。主要提供的变量为 0.55 μm 的 AOD 和细粒子比（Fine Mode Fraction，FMF），这些变量还会按算法分为 DT-ocean、DT-land 和 DB 等多个版本。MODIS 每次改进和更新算法后重新生产的产品集被称作一个 collection，目前的最新版是 2016 年发布的 C6.1，这一版本还引入了 3 x 3 km 高空间分辨率的 MOD04_3K 和 MYD04_3K 产品。

后面将会以 C6.1 版本的 MYD04_L2 产品为例，介绍每种算法提供的主要变量，以及用 Python 进行读取的方法。选择 MYD 是因为 Terra 发射时间更早仪器衰减更大，Aqua 的产品可能更为可靠（参考 Dark Target FAQ）。同时因为 DT 算法是气溶胶产品的主要算法，并且其文档更为详尽，所以 DB 算法的部分会比较简略。

MYD04_L2 文件名的意义如下图所示

文件格式是基于 HDF4 的 HDF-EOS2，因为自带经纬度坐标，所以不需要像陆面产品那样需要用 MRT 或 HEG 进行重投影，直接读了用就行。Python 中可以用 GDAL 或 pyhdf 包进行读取，这里使用后者。反演得到的参数以科学数据集（Scientific Data Set，SDS）的形式存储在 HDF 文件中，所谓 SDS 不过是带有元数据的数组罢了，而 pyhdf 的 SD 模块能提供关于 SDS 的 API。其使用方法详见 官方文档 和参考资料里的链接，不过写代码之前，可以先用 NASA 开发的 Panoply 软件直接打开 MYD04_L2 文件，查看 SDS 的元数据，并进行简单的可视化，如下面两张图所示

以 Optical_Depth_Land_And_Ocean 变量为例，可以看到其维度是沿卫星轨道的 203 个像元和横向的 135 个像元，这样一个数组画在等经纬度地图上呈边缘带弧度的矩形区域，同时有很多像元因为缺测而没有标出颜色。值得注意的是，AOD 本应是 double 类型（对应于 np.float64）的数据，但在 Panoply 中显示为 short 类型（对应于 np.int16），这是因为卫星产品常常会利用偏移和缩放操作将浮点类型的数据转换为整数类型，牺牲一点精度以压缩存储空间。Panoply 在画图时会自动换算回去，但使用 pyhdf 时我们需要根据 add_offset 和 scale_factor 属性手动进行换算。HDF 文件的元数据中记录有换算公式

float value = scale_factor * (stored integer - add_offset)

不过换算前需要小心 _FillValue 属性标记的缺测部分。用 Panoply 探索一番后，就可以动手实现一个读取 MYD04_L2文件的类

from pyhdf.SD import SD, SDC
import numpy as np
import pandas as pd

class ReaderMYD:
    '''读取MYD04_L2文件的类.'''
    def __init__(self, filepath):
        '''读取文件.'''
        self.sd = SD(filepath, SDC.READ)

    def __enter__(self):
        return self

    def search_sds(self, keyword):
        '''搜索SDS名,无视大小写.'''
        sdsnames = sorted(self.sd.datasets().keys())
        for sdsname in sdsnames:
            if keyword.lower() in sdsname.lower():
                print(sdsname)

    def read_lonlat(self):
        '''读取经纬度.'''
        lon = self.sd.select('Longitude')[:]
        lat = self.sd.select('Latitude')[:]

        return lon, lat

    def read_sds(self, sdsname, scale=True):
        '''读取SDS.'''
        sds = self.sd.select(sdsname)
        data = sds[:]

        # 按需进行缩放.
        if scale:
            attrs = sds.attributes()
            fill_value = attrs['_FillValue']
            scale_factor = attrs['scale_factor']
            add_offset = attrs['add_offset']
            data = np.where(
                data == fill_value,
                np.nan,
                (data - add_offset) * scale_factor
            )

        return data

    def read_time(self):
        '''读取沿swath的TAI时间,返回DatetimeIndex.'''
        second = self.sd.select('Scan_Start_Time')[:, 0]
        time = pd.to_datetime(second, unit='s', origin='1993-01-01')

        return time

    def close(self):
        '''关闭文件'''
        self.sd.end()

    def __exit__(self, *args):
        self.close()

这里使用 NaN 标记缺测值（参考 NumPy 系列：缺测值处理），用 __enter__ 和 __exit__ 实现上下文管理器（即可以使用 with 语句），读取的 TAI 时间与 UTC 时间仅相差数十秒，偷懒直接当作 UTC 时间用。

DT-ocean 算法及其变量

DT-ocean 的基本思路是，忽略可见光波段洋面贡献的反射率，那么大气顶的反射率可以表示为细模态（fine mode）和粗模态（coarse mode）气溶胶产生的反射率之和。所谓细模态指有效半径在 0.1 ~ 0.25 μm 范围内的细粒子，粗模态指有效半径在 1.0 ~ 2.5 μm 范围内的粗粒子。利用辐射传输模式模拟各种细模态气溶胶和粗模态气溶胶组合产生的反射率，用模拟值去逼近 MODIS 的实测值，最后求解出总 AOD、气溶胶的组合，以及细模态气溶胶的贡献比例。

具体反演过程基于查找表（Look-Up Table，LUT）法。首先根据以往的观测资料预定义四种典型的细模态气溶胶和五种典型的粗模态气溶胶，如下表所示

表中给出了每种气溶胶的折射系数和谱分布参数，利用 Mie 散射模型可以计算气溶胶在多个波段的 AOD 性质，利用辐射传输模式又可以计算出气溶胶存在时的大气顶反射率、大气透过率、大气后向散射比等参数，这些预先计算好的信息便构成了 LUT。接下来从四种细模态气溶胶和五种粗模态气溶胶中各选一种构成二十对组合，任意一对组合的总反射率表示为 ρλLUT(τ0.55)=η0.55ρλf(τ0.55)+(1−η0.55)ρλc(τ0.55) \rho_{\lambda}^{LUT}(\tau_{0.55}) = \eta_{0.55}\rho_{\lambda}^{f}(\tau_{0.55}) + (1 - \eta_{0.55})\rho_{\lambda}^{c}(\tau_{0.55}) ρλLUT​(τ0.55​)=η0.55​ρλf​(τ0.55​)+(1−η0.55​)ρλc​(τ0.55​) 其中 ρλf\rho_{\lambda}^{f}ρλf​ 和 ρλc\rho_{\lambda}^{c}ρλc​ 分别表示细模态和粗模态气溶胶的反射率，二者的值可以用 0.55 μm AOD τ0.55\tau_{0.55}τ0.55​ 在 LUT 中查询到。细粒子比 η0.55\eta_{0.55}η0.55​ 用于衡量细模态和粗模态气溶胶的贡献，值为 0 时表示全为细粒子，值为 1 时表示全为粗粒子。ρλLUT\rho_{\lambda}^{LUT}ρλLUT​ 即模拟出的气溶胶的总反射率，已知 MODIS 在 0.55、0.65、0.86、1.24、1.63 和 2.11 μm 六个波段测得大气层顶反射率 ρλm\rho_{\lambda}^{m}ρλm​，通过最小化 ρλLUT\rho_{\lambda}^{LUT}ρλLUT​ 和 ρλm\rho_{\lambda}^{m}ρλm​ 之间的拟合误差即可得到一对 τ0.55\tau_{0.55}τ0.55​ 和 η0.55\eta_{0.55}η0.55​ 的值——这便是反演的主要产品：0.55 μm 的气溶胶光学厚度和细粒子比。之后在 LUT 中很容易查询到其它波段 τ\tauτ 和 η\etaη 的值，这些值又可以用来得出其它变量，例如用 0.55 和 0.86 μm 的 τ\tauτ 计算 Ångström 指数（AE）。

二十对气溶胶组合能反演出二十组结果，我们肯定要进行挑选，其中拟合误差最小的结果称作 best solution，而误差较小的所有组合平均之后的结果称作 average solution。于是最后的产品里含有 best 和 average 两种类型的反演结果，官方推荐使用后者。

为了表示反演结果的可信度，每个像元被赋予了 QA（Quality Assuracne）标记，其中包含像元状态、云的影响、反演中遇到的问题等多种信息，不过我们一般只需要用到最基本的 QA Confidence（QAC）值：

• QAC = 3 表示结果非常好（very good）。
• QAC = 2 表示结果很好（good）。
• QAC = 1 表示结果有些勉强（marginal）。
• QAC = 0 表示结果很差（bad）。

画图时为了数据的完整性可以不考虑 QA，但进行定量统计时官方建议仅采用 QAC > 0 的像元。

简叙完必要的概念，下面介绍 DT-Ocean 的变量，并用颜色标出了比较重要的几个

Effective_Optical_Depth_Average_Ocean：Average solution 得到的 0.47、0.55、0.65、0.86、1.24、1.63 和 2.11 μm 七个波段的 AOD。被云遮挡的部分、陆面部分等都默认缺测。虽然表中没提，但文件中存在已经提取好的 0.55 μm AOD Effective_Optical_Depth_0p55um_Ocean。你可能会问为什么有个 Effective 的前缀？我也不清楚，可能是自老版本产品传承下来的约定。

Optical_Depth_Ratio_Small_Ocean_0.55micron：Average solution 得到的 0.55 μm 的细粒子比 FMF，表示细模态 AOD 与总 AOD 的比值，可以用来计算细模态 AOD 和粗模态 AOD τ0.55f=η0.55τ0.55τ0.55c=(1−η0.55)τ0.55 \tau_{0.55}^{f} = \eta_{0.55}\tau_{0.55} \newline \tau_{0.55}^{c} = (1-\eta_{0.55})\tau_{0.55} τ0.55f​=η0.55​τ0.55​τ0.55c​=(1−η0.55​)τ0.55​ Optical_Depth_Small_Average_Ocean：Average Solution 得到的七个波段的细模态 AOD，相当于提前帮你准备好了细模态 AOD，不需要用 FMF 手动进行计算。

Optical_Depth_Large_Average_Ocean：同上，提前准备好的粗模态 AOD。

Quality_Assurance_Ocean：QA 标记，数据类型为 short，QA 相关的信息需要通过位运算从 bit 中提取（详见 MODIS Atmosphere QA Plan），例如 QAC 可以这样提取

QAC = (Quality_Assurance_Ocean[:, :, 0] & 14) >> 1

这种操作还是挺麻烦的，所以文件中的 Land_Ocean_Quality_Flag 变量直接提供 0 ~ 3 的 QAC 值，其中还含有后面会介绍的 DT-land 的 QAC。例如筛去低质量的 AOD

Effective_Optical_Depth_0p55um_Ocean[Land_Ocean_Quality_Flag == 0] = np.nan

下面以 2015 年 2 月 11 日扫过中国东南沿海地区的 MYD04_L2.A2015042.0535.061.2018047210001.hdf 文件为例，画出变量 0.55 μm 的 AOD 和 FMF。首先在 EOSDIS Worldview 网站在线查看当天 Aqua MODIS 的真彩色图

可见沿海地区被“灰霾”覆盖，说明有气溶胶存在。接着在 Python 中用 Matplotlib 和 cartopy 包画图

import copy

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator
import cartopy.crs as ccrs

from reader import ReaderMYD
from map_funcs import set_map_extent_and_ticks

# 读取经纬度,时间和AOD变量.
filepath = './MYD04_L2.A2015042.0535.061.2018047210001.hdf'
with ReaderMYD(filepath) as f:
    lon, lat = f.read_lonlat()
    time = f.read_time()
    aod = f.read_sds('Effective_Optical_Depth_0p55um_Ocean')
    fmf = f.read_sds('Optical_Depth_Ratio_Small_Ocean_0.55micron')

# 设置地图.
extents = [100, 130, 10, 40]
proj = ccrs.PlateCarree()
fig, axes = plt.subplots(
    1, 2, figsize=(10, 6), subplot_kw={'projection': proj}
)
for ax in axes:
    ax.coastlines(resolution='10m', lw=0.3)
    set_map_extent_and_ticks(
        ax, extents,
        xticks=np.arange(-180, 190, 10),
        yticks=np.arange(-90, 100, 10),
        nx=1, ny=1
    )
    ax.tick_params(labelsize='small')

# 缺测设为灰色.
cmap = copy.copy(plt.cm.jet)
cmap.set_bad('lightgray')

# 绘制AOD.
im = axes[0].pcolormesh(
    lon, lat, aod[:-1, :-1], cmap=cmap, vmin=0, vmax=2,
    shading='flat', transform=proj
)
cbar = fig.colorbar(
    im, ax=axes[0], pad=0.1,
    ticks=MultipleLocator(0.5), orientation='horizontal'
)
cbar.set_label('AOD', fontsize='medium')
cbar.ax.tick_params(labelsize='small')
axes[0].set_title('Effective_Optical_Depth_0p55um_Ocean', fontsize='medium')

# 绘制FMF
im = axes[1].pcolormesh(
    lon, lat, fmf[:-1, :-1], cmap=cmap, vmin=0, vmax=1,
    shading='flat', transform=proj
)
cbar = fig.colorbar(
    im, ax=axes[1], pad=0.1,
    ticks=MultipleLocator(0.2), orientation='horizontal'
)
cbar.set_label('FMF', fontsize='medium')
cbar.ax.tick_params(labelsize='small')
axes[1].set_title(
    'Optical_Depth_Ratio_Small_Ocean_0.55micron', fontsize='medium'
)

# 用平均过境时间作为标题.
time_str = time.mean().strftime('%Y-%m-%d %H:%M')
fig.suptitle(time_str, fontsize='large')

plt.show()

代码中 set_map_extent_and_ticks 函数的定义请见 Cartopy 系列：从入门到放弃 一文。图中 Aqua 卫星于 UTC 时间 05:37，即北京时间下午 13:37 过境，MODIS 在黄海反演出 1.0 左右的 AOD 和 0.8 以上的 FMF，暗示主要为粒径偏小的人为排放气溶胶。同时台湾北侧因为有云遮挡而全部缺测。

DT-land 算法及其变量

DT-land 算法的基本思路与 DT-ocean 一致，但因为陆面在可见光和短波红外波段的反射率明显不为零，所以为了排除地表信号的影响，需要额外估算地表反射率。同时因为用到的观测波段变少了，采用的气溶胶模型也不同于 DT-ocean。

首先，DT-land 算法主要使用 0.47、0.65 和 2.12 μm 三个波段的观测数据，波段太少难以像 DT-ocean 那种动态选取细模态和粗模态气溶胶进行组合。取而代之的是一种更为宽泛的表示方法：fine model 和 coarse model。一个 model 由多个模态组成，反映了多种气溶胶的综合效应，fine model 即细模态气溶胶占主导地位的 model，coarse model 同理。DT-land 算法中将气溶胶概括以下五种 model

前四种是 fine model，第五种即沙尘是单独的 coarse model。其中 continental model 仅用于后面将会提到的地表偏亮情况时的反演。利用聚类分析等统计方法可以从 AERONET 站点的观测数据中得出全球不同区域各季节占主导地位的 fine model，之后的反演会参考这一结果预先选取合适的 fine model，再与沙尘 model 相组合。通过辐射传输模式计算的每种 model 产生的大气顶反射率等信息会被存入 LUT 中。

其次，地表反射率的估计有赖于 VISvs2.12 关系。Kaufmann 等通过观测和模拟研究发现，在植被覆盖或深色土壤的地表，可见光与短波红外波段的地表反射率的比值近乎定值，例如 0.47vs2.12 约为 0.25，0.65vs2.12 约为 0.5，这被称作 VISvs2.12 关系。DT-land 使用的 VISvs2.12 关系用散射角和短波短波红外 NDVI（衡量地表绿度）进行了修正，若已知 2.12 μm 的地表反射率，便可以根据该关系计算出 0.47 和 0.65 μm 波段的值。

对于暗地表，fine model、coarse model 和地表三者产生的大气顶反射率表示为 ρλ∗=ηρλ∗f+(1−η)ρλ∗c \rho_{\lambda}^{*} = \eta\rho_{\lambda}^{*f} + (1-\eta)\rho_{\lambda}^{*c} ρλ∗​=ηρλ∗f​+(1−η)ρλ∗c​ 其中 ρλ∗f\rho_{\lambda}^{*f}ρλ∗f​ 是 fine model 的反射率 ρλf\rho_{\lambda}^{f}ρλf​ 与地表反射率 ρλs\rho_{\lambda}^{s}ρλs​ 的综合效果，其数值可以用 τ0.55\tau_{0.55}τ0.55​ 在 LUT 中查询得到，ρλ∗c\rho_{\lambda}^{*c}ρλ∗c​ 同理。细粒子比 η\etaη（或 η0.55\eta_{0.55}η0.55​） 用于衡量 fine model 和 coarse model 的贡献。需要注意，除了 τ0.55\tau_{0.55}τ0.55​ 和 η0.55\eta_{0.55}η0.55​，ρλs\rho_{\lambda}^{s}ρλs​ 也是要求解的量，不同波段的值由 VISvs2.12 关系相关联。求解方程为 ρ0.47∗−ρ0.47m=0ρ0.66∗−ρ0.66m=ϵρ2.12∗−ρ2.12m=0 \rho_{0.47}^{*} - \rho_{0.47}^{m} = 0 \newline \rho_{0.66}^{*} - \rho_{0.66}^{m} = \epsilon \newline \rho_{2.12}^{*} - \rho_{2.12}^{m} = 0 ρ0.47∗​−ρ0.47m​=0ρ0.66∗​−ρ0.66m​=ϵρ2.12∗​−ρ2.12m​=0 即让 0.47 和 2.12 μm 波段的模拟值与实测值恰好相等，然后找出最小化 0.66 μm 波段误差 ϵ\epsilonϵ 的 τ0.55\tau_{0.55}τ0.55​、η0.55\eta_{0.55}η0.55​ 和 ρ2.12s\rho_{2.12}^{s}ρ2.12s​——此即反演的主要产品：0.55 μm 的 AOD 和 FMF，以及 2.12 μm 的地表反射率。之后可以利用 LUT 查询其它波段的值，并计算一系列的衍生量。

为了确保解的唯一性，η0.55\eta_{0.55}η0.55​ 仅取 -0.1 ~ 1.1 之间间隔为 0.1 的离散值，非物理的 -0.1 和 1.1 是为了在反演时容纳误差，反演结束后会被设为 0 和 1。另外当 τ0.55<2\tau_{0.55} < 2τ0.55​<2 时，η0.55\eta_{0.55}η0.55​ 准确度非常差，所以直接设为缺测。

对于不那么暗、稍微有点亮的地表也能进行反演，此时会设定 η0.55=1\eta_{0.55} = 1η0.55​=1，气溶胶模型仅由 continental model 构成，然后直接求解 0.47 和 2.12 μm 波段模拟值与实测值的方程。因为 VISvs2.12 关系在这种地表的误差会增大，导致反演结果不准，所以走这一流程的像元会被打上 QAC = 0 的标记。

DT-land 一个很特殊的地方在于，允许反演的 AOD 存在负值。因为反演在 0 附近存在误差，如果把值截断在 0 反而会使数据有偏，所以统计意义上有负值会更好。不过值最低也只允许到 -0.05 罢了。

因为陆面远比洋面复杂多变，所以 DT-land 的误差要比 DT-ocean 大。C5 版本的验证结果表示全球尺度上 DT-ocean AOD 的误差为 ±(0.03 + 5%)，DT-land AOD 的误差为 ±(0.05 + 15%)。验证还发现，陆面上的 FMF 和 AE 不太靠谱，仅有定性的价值，用其进行论证时需要留心。

DT-land 的 QA 同 DT-ocean，官方推荐使用 QAC = 3 的像元进行定量研究，这要求就比 DT-ocean 高很多，也可以理解为是陆面反演误差大，需要更强的限制。

因为 DT-land 的坑很多所以上面写的零零散散的。接着具体介绍变量

Corrected_Optical_Depth_Land：0.47、0.55 和 0.66 μm 三个波段的 AOD。跟 DT-ocean 类似，有个意义不明的 Corrected 的前缀。

Optical_Depth_Ratio_Small_Land：0.55 μm 的 FMF，数值是间隔为 0.1 的离散值。同时如前所述，因为陆面的 FMF 和 AE 不靠谱，所以最新的 C6.1 版本里陆面的 AE 和 fine model AOD 惨遭移除，如果有需求必须手动计算。

Quality_Assurance_Land：QA 标记，推荐直接使用 Land_Ocean_Quality_Flag，例如

Corrected_Optical_Depth_Land[:, Land_Ocean_Quality_Flag < 3] = np.nan

还是 DT-ocean 中那个例子，画出 0.55 μm 的 AOD 和 FMF

可以看到江苏省区域 AOD 较大FMF 偏高。同时注意到西边 AOD 较低的区域 FMF 直接缺测。

DT-ocean 和 DT-land 联合的变量

为了便于使用，MYD04_L2 提供两个把 DT-ocean 和 DT-land 结果合并的变量：

• Image_Optical_Depth_Land_And_Ocean：把 0.55 μm 的陆面 AOD 和 average solution 的洋面 AOD 直接拼在一起。因为没做任何质量控制，所以有值的像元很多，适合做图，这也呼应了其名字中的 image 一词。
• Optical_Depth_Land_And_Ocean：类似于前者，但洋面要求 QAC > 0，陆面要求 QAC = 3。虽然会使很多像元缺测，但更适合用于定量研究。

以前面的例子做图

这下终于能在地图上连续展示 AOD 的水平分布了。右图中陆面有值的像元数量有所缩水，但还可以接受。

DB 算法及其变量

沙漠和干旱地区的地表反射率在可见光波段很高，但在 412 nm 的深蓝波段又变得很低，DB 算法便是利用这一信息实现了亮地表的反演，同时也适用于一般的植被覆盖区域。具体算法详见 Hsu 等的一系列论文，类似于 DT-land，同样需要选取气溶胶模型、构建 LUT、估算地表反射率、模拟大气层顶的反射率等。DB 算法主要基于 412、470 和 650 nm 波段的观测，反演的主要变量是 AOD（550 nm） 和 AE（亮地表对应 412/470，暗地表对应 470/670）。后者定义为 α=ln⁡(τλ1/τλ2)ln⁡(λ1/λ2) \alpha = \frac{\ln(\tau_{\lambda_1}/\tau_{\lambda_2})}{\ln(\lambda_1/\lambda_2)} α=ln(λ1​/λ2​)ln(τλ1​​/τλ2​​)​ 一般来说 AE 小于 1 暗示粗粒子，大于 1 暗示细粒子。下面这张图展示了常见气溶胶的 AOD 和 AE

DB 同样有 QA 标记，与 DT 的区别在于 0 直接表示缺测。

肯定会有人问，既然 DT 和 DB 都能反演陆面，那么选哪个好？据 Sayer 等的验证研究，全球陆面上二者反演的 AOD 及其与 AERONET 观测的一致性都很相近，虽然 DT 与 AERONET 的相关性更好，但 DB 可以反演的空间范围更广，并且在低 AOD 情况下的误差更小。也许在一些特定的场合某个算法表现会更好，但一般任选一个使用就行，详细的区别还请自行查阅。

DB 相关的主要变量如下表所示

Deep_Blue_Aerosol_Optical_Depth_550_Land：550 nm 的 AOD。

Deep_Blue_Aerosol_Optical_Depth_550_Land_Best_Estimate：同上，但是仅选取 QA > 1 的值，官方将其描述为“大部分用户应该使用的量”。

Deep_Blue_Aerosol_Optical_Depth_550_Land_QA_Flag：QA 标记，不需要解码直接使用即可。

Deep_Blue_Angstrom_Exponent_Land：AE，亮地表对应于 412/470，暗地表对应于 470/670。

接下来画出 AOD 和 AE 看看

AOD 的水平分布与 DT-land 一致，但有值的像元增加了很多。山东和江苏区域的 AOD 值很高，甚至能超过 2.0，而 1.7 左右的 AE 暗示气溶胶粒径偏小。

DT 和 DB 联合的变量

介绍完前面五花八门的变量，终于到了集大成的时刻——将 DT-ocean、DT-land 和 DB 算法反演出来的最好的 0.55 μm AOD 合并成一个变量。这一联合变量于 C6 版本引入，目前为 Aqua 卫星专供，合并的策略并不复杂：

• 洋面直接使用 QA_DT > 0 的 DT-ocean AOD。
• 对于陆面，当 NDVI < 0.2，即植被稀疏时，采用 QA_DB > 1 的 DB AOD；当 NDVI > 0.3，即植被茂密时，采用 QA_DT = 3 的 DT-land AOD；当 NDVI 介于中间时：
  • QA_DB > 1 且 QA_DT = 3，则取两种 AOD 的平均值。
  • 若有一方的 QA 不满足上一条，则选取满足的那一方。
  • 若双方的 QA 都不达标，则设为缺测。

合并得到的 AOD 也有 QA，数值基本继承自 QA_DT 和 QA_DB，和 QA_DB 一样用 0 表示缺测。不过既然合并策略里已经进行了质量控制，这个 QA 略显多余。另外需要注意这个量的可靠性目前还没有得到系统性的验证。

AOD_550_Dark_Target_Deep_Blue_Combined：顾名思义，DT 和 DB 合并得到的 0.55 μm AOD。

AOD_550_Dark_Target_Deep_Blue_Combined_Algorithm_Flag：像元具体采用了哪个算法，0 表示 DT，1 表示 DB，2 表示混合（即平均）。

下面画出这两个量

可以看出中国东南区域的合并 AOD 主要基于 DT，到了北边才开始使用 DB。

格点化处理

本节以 2021 年 3 月 15 日中国发生的沙尘暴天气为例，演示如何把一系列 MYD04_L2 文件里的 AOD 数据格点化并画出来。AOD 就选取上一节介绍的 AOD_550_Dark_Target_Deep_Blue_Combined，在南方植被茂密的区域会选用成熟的 DT 算法，而在北方的干旱区域又能发挥 DB 算法的优势，并且只含高 QA 的结果。下图为当天的 Aqua MODIS 真彩色图

中国北方一带明显可见浓厚的黄色沙尘，甚至还卷入了东边的气旋中，不难预测这些区域会有特别高的 AOD。下载 MYD04_L2 产品可以到 NASA 的 EARTHDATA 网站，选中这一天，经纬度范围设为 70°E ~ 140°E，10°N ~ 60°N，最后会筛选出十多个文件。代码和结果如下

from pathlib import Path

import numpy as np
from scipy.stats import binned_statistic_2d
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from matplotlib.ticker import MultipleLocator
import cartopy.crs as ccrs

from reader import ReaderMYD
from map_funcs import add_Chinese_provinces, set_map_extent_and_ticks

def grid_data(x, y, data, xbins, ybins):
    '''
    利用平均的方式格点化一维散点数据.

    没有数据的格点记作NaN.若格点内的数据全部缺测会产生warning,可以无视.
    为了便于画图,结果的维度为(ny, nx).
    '''
    avg, ybins, xbins, _ = binned_statistic_2d(
        y, x, data, bins=[ybins, xbins], statistic=np.nanmean
    )
    xc = (xbins[1:] + xbins[:-1]) / 2
    yc = (ybins[1:] + ybins[:-1]) / 2

    return xc, yc, avg

def aod_cmap():
    '''制作适用于AOD的cmap.'''
    rgb = np.loadtxt('./NEO_modis_aer_od.csv', delimiter=',')
    cmap = ListedColormap(rgb / 256)

    return cmap

if __name__ == '__main__':
    # 读取MYD04_L2文件,收集所有granule的经纬度和AOD数据.
    dirpath = Path('./data')
    lon_all, lat_all, aod_all = [], [], []
    for filepath in dirpath.iterdir():
        with ReaderMYD(str(filepath)) as f:
            lon, lat = f.read_lonlat()
            aod = f.read_sds('AOD_550_Dark_Target_Deep_Blue_Combined')
        lon_all.append(lon.ravel())
        lat_all.append(lat.ravel())
        aod_all.append(aod.ravel())
    # 连接为一维数组.
    lon_all = np.concatenate(lon_all)
    lat_all = np.concatenate(lat_all)
    aod_all = np.concatenate(aod_all)

    # 设定网格.
    extents = [70, 140, 10, 60]
    lonmin, lonmax, latmin, latmax = extents
    dlon, dlat = 1, 1
    lon_bins = np.arange(lonmin, lonmax + 0.5 * dlon, dlon)
    lat_bins = np.arange(latmin, latmax + 0.5 * dlat, dlat)
    # 格点化.
    lon_grid, lat_grid, aod_grid = grid_data(
        lon_all, lat_all, aod_all, lon_bins, lat_bins
    )

    # 设置地图.
    proj = ccrs.PlateCarree()
    fig = plt.figure()
    ax = fig.add_subplot(111, projection=proj)
    add_Chinese_provinces(ax, lw=0.3, ec='k', fc='none')
    ax.coastlines(resolution='10m', lw=0.3)
    set_map_extent_and_ticks(
        ax, extents,
        xticks=np.arange(-180, 190, 10),
        yticks=np.arange(-90, 100, 10),
        nx=1, ny=1
    )
    ax.tick_params(labelsize='x-small')

    # 画出格点化的AOD.
    cmap = aod_cmap()
    cmap.set_bad('lightgray')   # 缺测设为灰色.
    im = ax.pcolormesh(
        lon_grid, lat_grid, aod_grid, cmap=cmap, vmin=0, vmax=3,
        shading='nearest', transform=proj
    )
    cbar = fig.colorbar(
        im, ax=ax, ticks=MultipleLocator(0.5),
        extend='both', shrink=0.8, pad=0.1, aspect=30,
        orientation='horizontal'
    )
    cbar.set_label('AOD', fontsize='small')
    cbar.ax.tick_params(labelsize='x-small')

    # 日期作为标题.
    ax.set_title('2021-03-15', fontsize='medium')

    # 保存图片.
    fig.savefig('gridded.png', dpi=200, bbox_inches='tight')
    plt.close(fig)

格点化采用计算落入格子中的数据点的平均值的方式，可以通过写循环数格子来实现，不过这里选择直接调用 scipy.stats.binned_statistic_2d 函数。该函数能统计一维散点数据在二维数据框（bin）中的统计值，在计算平均值时会正确地将空格点设为 NaN。此外为了美观，colormap 选用 Panoply Additional Color Tables 中的 NEO_modis_aer_od，颜色从米黄向深红过渡，个人觉得很适合展示 AOD。不过下载到的 colormap 是 Adobe 的 ACT 格式，需要用这篇 Stack Exchange 回答 中的代码转换为 CSV 格式以方便读取。

画出来的结果不出所料，一条长长的“沙龙”从新疆延伸到东北，AOD 超过了 3.0，这算是非常高的值了，暗示会给沿途的地面带来强烈的污染。中国以北因为存在大片云系所以基本都缺测了。调查资料可知这次沙尘天气的两大成因：一是我国西北和蒙古国地区前期气温偏高、降水稀少，因而地表容易起沙；二是蒙古气旋带来的大风为沙尘的扬起和传输提供了动力条件。

MODIS 气溶胶产品提供的变量众多，抛开前文冗长的讲解，这里给出极简使用指南。

想画 AOD 图：

• Image_Optical_Depth_Land_And_Ocean

想定量研究 AOD：

• Optical_Depth_Land_And_Ocean

专注于洋面上的气溶胶：

• Effective_Optical_Depth_0p55um_Ocean
• Optical_Depth_Ratio_Small_Ocean_0.55micron
• Land_Ocean_Quality_Flag

专注陆面上的气溶胶：

• Corrected_Optical_Depth_Land
• Optical_Depth_Ratio_Small_Land
• Land_Ocean_Quality_Flag

需要研究沙漠和干旱地区：

• Deep_Blue_Aerosol_Optical_Depth_550_Land
• Deep_Blue_Angstrom_Exponent_Land
• Deep_Blue_Aerosol_Optical_Depth_550_Land_QA_Flag

想结合暗目标算法和深蓝算法的结果：

• AOD_550_Dark_Target_Deep_Blue_Combined

如果你需要用到 C6.1 之前的产品，例如 C5 产品，很可能会发现许多变量名字不一样，甚至是找不到。不过，要是你理解了各个算法的基本概念，再多去官网搜搜，一定能推理并找出所需的变量。就算以后 MODIS 产品有了新版本，或者说要面对其它采用了同样算法的卫星产品（例如 VIIRS、NPP）也应该不会有问题。

MODIS 及其气溶胶产品的基本介绍：

Wikipedia: Moderate Resolution Imaging Spectroradiometer

Atmosphere Discipline Team Imager Products: Aerosol (04_L2)

The Collection 6 MODIS aerosol products over land and ocean

Understanding the Details of the MODIS Aerosol Product（某个 Training 的 PPT，我搜不到链接了）

DARK TARGET AEROSOL PRODUCTS USER’S GUIDE

DT 和 DB 算法小组的官网，详细介绍了算法和数据使用的建议：

Dark Target Website

Deep Blue Website

利用 Python 读取 MYD04_L2 文件的教程：

How to read a MODIS HDF4 file using python and pyhdf ?

ARSET - Data Analysis Tools for High Resolution Air Quality Satellite Datasets

MORE PYTHON COMPREHENSIVE EXAMPLES FOR LAADS MOD PRODUCTS