aodkit 0.2.0

imutum's packages for aerosol optical depth retrieval

aodkit

AOD (Aerosol Optical Depth) 卫星产品验证工具库。提供卫星-地面站时空匹配、精度指标计算、科学制图等功能。

安装

pip install aodkit
# 可选：CJK 字体自动提取
pip install aodkit[plot]

卫星-地面站匹配 (matchup)

基于 Du et al. 2025 (Remote Sens. 17:1235) 提出的三种空间匹配方法，用于评估 AOD 反演算法精度。

原理

卫星产品与地面站观测在光谱、时间、空间上存在差异，需要配准后才能对比。完整流程：

光谱配准 → 时间匹配 (±30min) → 空间匹配 → 精度评估

空间匹配是核心环节。对同一组卫星像素和地面站观测，不同的空间匹配策略会产生不同的验证结论。 本库实现了论文中的三种方法：

直接法 (Direct, Eq.1) — 选取距站点最近的单个像素：

$$x_{\text{direct}} = x_j, \quad j = \arg\min_{i \in W} d(\text{pixel}_i, \text{site})$$

不做 QA 筛选，直接记录该像素的 AOD 和 QA 值，后续再进行质量过滤。 最简单、受单像素噪声影响最大。

均值法 (Average, Eq.2) — 空间窗口内 QA 筛选后取均值：

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

$W$ 为空间窗口与 QA 筛选的交集。NASA 标准验证方法（DT: QA=3, DB: QA≥2）。 通过空间平均降低随机误差，但亚像素云污染会引入系统偏差。

最优法 (Optimal, Eq.3) — 空间窗口内 QA 筛选后，选与地面真值绝对误差最小的像素：

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

后验评估方法，需要已知站点真值 $x_{\text{site}}$（来自时间匹配后的 AERONET 观测）。 结果代表反演算法的精度上限——如果能完美选像素，算法能有多准。 不能替代均值法，但可以揭示 QA 筛选机制的改进空间。

三步流水线

原始 2D 数组 (aod, qa, lat, lon)
  │
  ├── Step 1: pos_mask = region.make_mask(shape)     → 空间范围 mask
  │
  ├── Step 2: qa_mask = 用户自定义 QA 条件            → 质量筛选 mask
  │            (直接法跳过此步)
  │
  └── Step 3: aggregate_*(aod, qa, pos_mask, ...)    → (aod值, 代表性qa, 像素数)

每一步只做一件事：缩小有效像素集合。聚合函数在最终的有效像素上操作。

空间区域

两种数据类型统一为 make_mask() 接口：

from aodkit.matchup import WindowRegion, RadiusRegion

# ── 网格化产品 (MCD19A2)：已知站点像素索引，无需 lat/lon ──
region = WindowRegion(600, 800, "5x5")

# ── 网格化产品：通过经纬度查找最近像素 ──
region = WindowRegion.from_latlon(lat, lon, site_lat, site_lon, "5x5")

# ── L2 刈幅产品 (DT/DB/VIIRS)：半径匹配 ──
region = RadiusRegion(lat, lon, site_lat, site_lon, radius_km=25)

# 统一接口
pos_mask = region.make_mask(aod.shape)  # bool 数组，True = 在空间范围内

QA 筛选

用户在外部生成布尔掩码，聚合函数不关心 QA 编码方式：

# DT 产品 (标量 QA, 0-3, 越大越好)
qa_mask = qa >= 3               # 陆地推荐
qa_mask = qa >= 2               # 海洋推荐

# MCD19A2 (16-bit 位掩码)
cloud = (qa_bits >> 0) & 0b111   # bits 0-2: CloudMask
adj   = (qa_bits >> 5) & 0b111   # bits 5-7: AdjacencyMask
qa_mask = (cloud == 1) & (adj == 0)  # Clear + Normal

聚合

from aodkit.matchup import (
    aggregate_direct, aggregate_direct_grid,
    aggregate_average, aggregate_optimal,
)

# ── 直接法 ──
# L2 刈幅：haversine 距离找最近邻
aod_val, qa_val, n = aggregate_direct(aod, qa, pos_mask, lat, lon, site_lat, site_lon)

# 网格化：窗口中心像素
aod_val, qa_val, n = aggregate_direct_grid(aod, qa, pos_mask, center_row=600, center_col=800)

# ── 均值法 ──
mean_aod, rep_qa, n = aggregate_average(aod, qa, pos_mask, qa_mask)
# rep_qa = 最接近均值的像素的 QA (argmin|x_i - mean|)

# ── 最优法 ──
best_aod, best_qa, n = aggregate_optimal(aod, qa, pos_mask, qa_mask, site_aod=0.35)
# site_aod 来自时间匹配后的 AERONET 观测

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

完整示例：DT 产品验证

from aodkit.matchup import RadiusRegion, aggregate_average

# 读取卫星数据
aod = ...   # 2D AOD 数组 (e.g., Corrected_Optical_Depth_Land)
qa = ...    # 2D QA 数组 (e.g., Land_Ocean_Quality_Flag)
lat, lon = ...  # 2D 经纬度

# Step 1: 空间提取 (25km 半径)
region = RadiusRegion(lat, lon, site_lat=40.0, site_lon=116.0, radius_km=25)
pos_mask = region.make_mask(aod.shape)

# Step 2: QA 筛选 (DT 陆地推荐 QA=3)
qa_mask = qa >= 3

# Step 3: 均值法聚合
mean_aod, rep_qa, n_pixels = aggregate_average(aod, qa, pos_mask, qa_mask)

完整示例：MCD19A2 验证

from aodkit.matchup import WindowRegion, aggregate_optimal

# 读取 MCD19A2 数据 (正弦投影网格，无需 lat/lon)
aod = ...       # Optical_Depth_055
qa_bits = ...   # AOD_QA (16-bit)

# Step 1: 空间提取 (已知站点像素位置)
region = WindowRegion(center_row=600, center_col=800, window="5x5")
pos_mask = region.make_mask(aod.shape)

# Step 2: QA 筛选 (Best quality = CloudMask=Clear AND AdjacencyMask=Normal)
cloud = (qa_bits >> 0) & 0b111
adj   = (qa_bits >> 5) & 0b111
qa_mask = (cloud == 1) & (adj == 0)

# Step 3: 最优法聚合 (需要站点真值)
site_aod = 0.35  # 来自 AERONET 时间匹配
best_aod, best_qa, n = aggregate_optimal(aod, qa_bits, pos_mask, qa_mask, site_aod)
# best_qa 为该像素的原始 16-bit QA 值，可按需解码各位域

参考文献

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.