Metadata-Version: 2.4
Name: aodkit
Version: 0.2.3
Summary: imutum's packages for aerosol optical depth retrieval
License: MIT License
        
        Copyright (c) 2023 imutum
        
        Permission is hereby granted, free of charge, to any person obtaining a copy
        of this software and associated documentation files (the "Software"), to deal
        in the Software without restriction, including without limitation the rights
        to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
        copies of the Software, and to permit persons to whom the Software is
        furnished to do so, subject to the following conditions:
        
        The above copyright notice and this permission notice shall be included in all
        copies or substantial portions of the Software.
        
        THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
        IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
        FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
        AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
        LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
        OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
        SOFTWARE.
        
Classifier: Development Status :: 3 - Alpha
Classifier: Programming Language :: Python :: 3.10
Classifier: Intended Audience :: Developers
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: MacOS
Classifier: Operating System :: POSIX :: Linux
Classifier: Operating System :: Microsoft :: Windows
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: matplotlib
Requires-Dist: cycler
Requires-Dist: numpy
Requires-Dist: pandas
Requires-Dist: pyarrow
Requires-Dist: hdfkit
Requires-Dist: fonttools
Requires-Dist: scipy>=1.10
Requires-Dist: shapely>=2.0.0
Requires-Dist: geopandas>=0.13.2
Dynamic: license-file

# aodkit

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

## 安装

```bash
pip install aodkit
```

## 验证流程概览

卫星反演 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)。

```python
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)。

```python
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` |

```python
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 筛选

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

```python
# 整数 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)

```python
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)

```python
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.
