aodkit 0.3.0

imutum's packages for aerosol optical depth retrieval

aodkit

AOD (Aerosol Optical Depth) 卫星产品验证工具库。覆盖验证全流程：光谱插值、时间匹配、空间匹配、精度评估与科学制图。

安装

pip install aodkit

AERONET 地面真值数据库

验证所需的 AERONET 站点 AOD 数据需预先下载并处理为结构化格式。数据处理脚本位于 scripts/aeronet/。

• 数据库搭建指南 — 下载 → 拆分 → 550nm 插值 → 提取，含 Level 1.5 / 2.0 不可合并的实验依据
• 数据库使用指南 — 时间匹配与时空匹配示例

验证流程概览

卫星反演 AOD 与地基 AERONET 观测在光谱、时间、空间三个维度上存在不一致。标准验证流程 (Ichoku et al., 2002) 依次消除这三个差异，生成可对比的验证对 (matchup)：

$$\text{光谱插值};\xrightarrow{\text{§1}};\text{时间匹配};\xrightarrow{\text{§2}};\text{空间匹配};\xrightarrow{\text{§3}};\text{精度评估}$$

本库各模块对应关系：

步骤	模块	核心函数
光谱插值	aodkit.interp	interp_aod_at_wavelength()
时间匹配	aodkit.interp	time_match()
空间匹配	aodkit.matchup	WindowRegion / RadiusRegion + aggregate_*()
精度评估	aodkit.metrics / aodkit.comparison	calc_comparison_stats() 等
制图	aodkit.plot	plot_default_density_kernel_chart() 等

---

§1 光谱插值

卫星产品通常报告 550 nm AOD，而 AERONET 不直接观测该波长。需将 AERONET 多波段 AOD（440, 500, 675, 870 nm 等）插值至 550 nm (O'Neill et al., 2001)。

from aodkit.interp import interp_aod_at_wavelength

# dataframe: AERONET 数据，列名含 AOD_440nm, AOD_500nm, AOD_675nm...
aod_550 = interp_aod_at_wavelength(dataframe, method="numpydeg2_polyfit", wavelength=550)

支持 12 种插值方法：

类别	方法	说明
scipy.interpolate	linear, nearest, nearest-up, zero_spline, slinear_spline, quadratic_spline, cubic_spline	interp1d 逐行插值
numpy polyfit (log-log)	numpydeg2_polyfit, numpydeg3_polyfit	ln(AOD) vs ln(λ) 多项式拟合
scipy curve_fit	scipy_curvefit	$\tau(\lambda) = e^{a_0 + a_1\ln\lambda + a_2(\ln\lambda)^2}$
Ångström 指数	Angstrom	$\tau(\lambda) = \tau(440)\cdot(\lambda/440)^{-\alpha}$
复合	cubic+scipy	cubic_spline 与 scipy_curvefit 取均值

推荐使用 numpydeg2_polyfit：在 log-log 空间做二次多项式拟合，精度与计算速度平衡最优。

---

§2 时间匹配

卫星过境为瞬时观测，AERONET 每 5–15 min 记录一次。需在卫星过境时刻附近的时间窗口内聚合 AERONET 测量值 (Ichoku et al., 2002)。

from aodkit.interp import time_match

# srcdata: AERONET AOD 时序 (DatetimeIndex, UTC)
# objtime: 卫星过境时刻表 (含 timestamp 列, UTC 秒)
matched = time_match(srcdata, objtime, method="average", time_range=(-1800, 1800))
# 返回 DataFrame：每个目标时刻对应一行，附 counts 列 (命中的源样本数)

参数	说明
method="average"	窗口内取算术平均（NASA 标准：要求 ≥ 2–3 个有效观测）
method="linear"	在窗口内最近两侧样本间做线性插值
time_range=(-1800, 1800)	以目标时刻为 0 的秒数范围，默认 ±30 min

---

§3 空间匹配

3.1 三种方法

基于 Du et al. (2025) 的系统比较框架，本库实现三种空间匹配方法：

Direct — 直接法

选取距站点空间距离最近的单个卫星像素：

$$x_{\text{direct}} = x_j, \quad j = \arg\min_{i \in W}; d!\left((\varphi_i, \lambda_i),; (\varphi_{\text{site}}, \lambda_{\text{site}})\right) \tag{1}$$

其中 $d(\cdot)$ 为 Haversine 球面距离，$W$ 为空间窗口内的有效像素集合。不做 QA 筛选，直接记录目标像素的 AOD 与 QA 值，质量过滤在后续阶段完成。最简单，但完全依赖单像素，易受亚像素云污染和检索噪声影响。

Average — 均值法

对空间窗口内通过 QA 筛选的像素取算术平均：

$$x_{\text{average}} = \frac{1}{|W'|}\sum_{i \in W'} x_i \tag{2}$$

其中 $W' = W \cap {i \mid \text{QA}(i) \geq \text{threshold}}$。NASA Dark Target 标准验证方法 (Levy et al., 2013)，推荐阈值：DT 陆地 QA = 3，DB 陆地 QA $\geq$ 2。可降低随机误差，但亚像素云污染会引入系统性正偏差。

Optimal — 最优法

在 QA 筛选后的像素集合中，选取与地面真值绝对误差最小的像素：

$$x_{\text{optimal}} = x_j, \quad j = \arg\min_{i \in W'}; |x_i - x_{\text{site}}| \tag{3}$$

其中 $x_{\text{site}}$ 为时间匹配后的 AERONET AOD。后验评估方法：需要已知地面真值，不可用于业务检索，但回答关键问题——若 QA 与空间匹配均理想化，反演算法本身能达到什么精度？Du et al. (2025) 结果表明 MAIAC 1 km 产品在 Optimal 下 =EE 达 94%，揭示高分辨率像素纯度优势。

3.2 实现架构：三步流水线

Step 1  空间提取    pos_mask = region.make_mask(shape)       bool, True = 在空间范围内
Step 2  QA 筛选     qa_mask = 用户条件                        bool, True = 通过质量筛选 (Direct 跳过)
Step 3  聚合计算    aggregate_*(aod, qa, pos_mask, ...)       → (aod_value, qa_value, n_pixels)

空间区域

产品类型	区域类型	构造方式
L2 Swath (DT/DB/VIIRS)	RadiusRegion	经纬度 + 半径 (km)
网格化 (MCD19A2 等)	WindowRegion	像素索引 或 from_latlon

from aodkit.matchup import WindowRegion, RadiusRegion

region = RadiusRegion(lat, lon, site_lat, site_lon, radius_km=25)         # L2 Swath
region = WindowRegion(center_row=600, center_col=800, window="5x5")       # 网格化 (索引)
region = WindowRegion.from_latlon(lat, lon, site_lat, site_lon, "5x5")    # 网格化 (经纬度)

pos_mask = region.make_mask(aod.shape)

QA 筛选

用户根据产品规范生成布尔掩码，聚合函数不感知编码方式：

# 整数 QA — 直接阈值比较
qa_mask = qa >= 3

# 位掩码 QA — 用 hdfkit readbit() 提取 (dataset 名和位范围视产品规范而定)
field_a = reader.readbit("<qa_dataset>", bit_start, bit_end)
field_b = reader.readbit("<qa_dataset>", bit_start, bit_end)
qa_mask = (field_a == 1) & (field_b == 0)

聚合函数

三个函数对应三种方法，通过 region.nearest_index() 多态实现 Swath/Grid 统一：

函数	方法	AOD 输出	QA 输出
aggregate_direct(aod, qa, pos_mask, region)	Eq.1	最近像素值	该像素 QA
aggregate_average(aod, qa, pos_mask, qa_mask)	Eq.2	有效像素算术平均值	argmin|x_i - mean| 像素的 QA
aggregate_optimal(aod, qa, pos_mask, qa_mask, site_aod)	Eq.3	argmin|x_i - x_site| 像素值	该像素 QA

所有函数返回 (aod_value, qa_value, n_pixels)。无有效像素时返回 (nan, nan, 0)。

---

§4 完整示例

4.1 L2 DT 产品验证 (Average)

from aodkit.interp import interp_aod_at_wavelength, time_match
from aodkit.matchup import RadiusRegion, aggregate_average

# ── 光谱插值：AERONET 多波段 → 550nm ──
aeronet_aod_550 = interp_aod_at_wavelength(aeronet_df, method="numpydeg2_polyfit", wavelength=550)
aeronet_df["AOD_550"] = aeronet_aod_550

# ── 时间匹配：±30min 窗口取平均 ──
objtime = pd.DataFrame({"timestamp": [sat_overpass_utc]}, index=[pd.Timestamp(sat_overpass_utc)])
site_matched = time_match(aeronet_df[["AOD_550"]], objtime, method="average", time_range=(-1800, 1800))
site_aod = site_matched.iloc[0]["AOD_550"]

# ── 空间匹配：25km 半径 + QA 筛选 + Average 聚合 ──
region = RadiusRegion(lat, lon, site_lat=40.0, site_lon=116.0, radius_km=25)
pos_mask = region.make_mask(sat_aod.shape)
qa_mask = qa >= 3

mean_aod, rep_qa, n = aggregate_average(sat_aod, qa, pos_mask, qa_mask)

4.2 网格化产品验证 (Optimal)

from aodkit.matchup import WindowRegion, aggregate_optimal

# 空间匹配 (网格化产品，已知站点像素位置)
region = WindowRegion(center_row=600, center_col=800, window="5x5")
pos_mask = region.make_mask(aod.shape)

# QA 筛选 (位掩码 — 具体 dataset 名和位范围视产品规范而定)
field_a = reader.readbit("<qa_dataset>", bit_start, bit_end)
field_b = reader.readbit("<qa_dataset>", bit_start, bit_end)
qa_mask = (field_a == 1) & (field_b == 0)

# Optimal 聚合 (site_aod 来自 §1-§2 的光谱插值 + 时间匹配)
best_aod, best_qa, n = aggregate_optimal(aod, qa, pos_mask, qa_mask, site_aod)

---

参考文献

• Du, B.; Zhong, B.; Cai, H.; Wu, S.; et al. Improving AOD Algorithm Evaluation: A Spatial Matching Method for Minimizing Quality Control Bias. Remote Sens. 2025, 17, 1235.
• Holben, B.N.; Eck, T.F.; Slutsker, I.; et al. AERONET — A Federated Instrument Network and Data Archive for Aerosol Characterization. Remote Sens. Environ. 1998, 66, 1–16.
• Ichoku, C.; Chu, D.A.; Mattoo, S.; et al. A Spatio-Temporal Approach for Global Validation and Analysis of MODIS Aerosol Products. Geophys. Res. Lett. 2002, 29, MOD1-1–MOD1-4.
• Levy, R.; Mattoo, S.; Munchak, L.; et al. The Collection 6 MODIS Aerosol Products over Land and Ocean. Atmos. Meas. Tech. 2013, 6, 2989–3034.
• Lyapustin, A.; Wang, Y.; Korkin, S.; Huang, D. MODIS Collection 6 MAIAC Algorithm. Atmos. Meas. Tech. 2018, 11, 5741–5765.
• O'Neill, N.T.; Eck, T.F.; Holben, B.N.; et al. Bimodal Size Distribution Influences on the Variation of Ångström Derivatives in Spectral and Optical Depth Space. J. Geophys. Res. 2001, 106, 9787–9806.