MCD19A2气溶胶AOD批处理：HDF4读取、QA筛选与GDAL重投影

简介：面向遥感影像处理与地球科学数据分析的学习者，这套代码围绕MCD19A2产品展开，覆盖重投影、镶嵌、重采样、裁剪以及AOD栅格提取点值等常见环节，适合已具备一定IDL与Python基础、希望把零散处理流程整合为批量作业的初中级用户参考。压缩包共12个文件、约11KB，以3个py脚本、3个pro程序为主，另有xml与iml等工程配置文件和gitignore，其中脚本与程序负责数据处理逻辑，配置文件描述IDL工程环境，便于直接在已有工具链中复用。已有1248人学习下载，说明其在MODIS数据处理场景中具备一定参考价值。读者可从中获得MCTK批量处理思路、均值统计与栅格镶嵌的具体实现，以及Python与IDL混合协作的组织方式，用于对照排查自己流程中的坐标转换、重叠区处理等细节。

1. 从一堆积灰的 HDF4 说起：MCD19A2 数据处理代码到底解决什么问题

第一次拿到 MCD19A2 数据的人，十有八九会在 GIS 软件里点开一个.hdf文件，然后看到一片黑——不是数据坏了，是这套 1km 的 MAIAC 气溶胶产品躺在正弦投影里，还套着一层int16缩放和一堆看不见的 QA 位域。这份 MCD19A2 数据处理代码要干的事很朴素：把“一天一个 tile、一个 tile 十几个科学数据集”的原始文件，批量变成能直接参与分析的反演 AOD 栅格，以及按站点、按区域汇总的时序表。

它适合三类人：做气溶胶时序分析、需要长序列 AOD 的同学，要把 AOD 接进大气校正或辐射传输流程的算法工程师，以及手上已经攒了几百景文件、却还在手动点 GIS 菜单的同行。整条链路用 Python 串起来，pyhdf负责把 HDF4 掰开，GDAL 负责投影，numpy 和 pandas 负责筛选与统计，做成可复用的小工具之后，一天几千景也只是排队跑的问题。

2. 工程结构与依赖：从文件命名规则到可复现的目录骨架

2.1 先认清原料：文件名、图层与网格约定

在写任何一行代码之前，先把原料摸清楚，后面一大半的坑都能提前绕开。MCD19A2 的文件名是高度结构化的，典型形态长这样：

MCD19A2.A2020047.h26v04.006.2020052134512.hdf

拆开看：A2020047是年积日（2020 年第 47 天），h26v04是正弦投影的瓦片行列号，006是集合版本号。瓦片号直接决定了这个文件在地球上的位置，所以解析文件名不是一个“顺手写写”的小函数，它是整条链路的地理定位入口。我在项目里会把它单独抽出来，因为后面算仿射变换、做瓦片拼接、按 tile 分组批处理，全靠这两个数字。

文件内部是一组科学数据集（SDS），常用的几个：

网格约定是另一个必须记住的常量：全球正弦网格被切成 36×18 个瓦片，每个瓦片 1200×1200 像素，单像元边长 926.625433 m，一个瓦片的跨度约 1111950.52 m。整套数字后面会反复用到，建议直接定义成模块级常量，别散落在各处。

2.2 技术栈选型：为什么是 pyhdf + GDAL 而不是一把梭

处理 HDF4 的库有好几个，选择上踩过的坑比想象中多。rasterio底层也是 GDAL，但如果本机的 GDAL 在编译时没开 HDF4 驱动，rasterio.open()打开.hdf会直接抛not recognized as a supported file format，报错信息还完全看不出根因。h5py更不行，MCD19A2 是 HDF4 而不是 HDF5，两者除了名字像没有任何关系。

最后的组合是这样的：pyhdf负责把 SDS 读成 numpy 数组，这个库读 HDF4 最稳，属性和填充值都能拿到；GDAL 只负责重投影和写出，不参与读取，这样就绕开了 GDAL 版本对 HDF4 支持不一致的问题。中间产物是一个在内存里手工设好地理参考的数据集，投影和仿射变换由我们自己算，不依赖 GDAL 从 HDF-EOS 结构里推断——推断出来的结果在不同 GDAL 版本上偶有半像素级的差异，做长时序叠加时这就是条带。

环境搭建我一般用 conda，pyhdf用 pip 装经常卡在编译 HDF4 上，conda-forge 有现成二进制包：

conda create -n mcd19 python=3.10 -y conda activate mcd19 conda install -c conda-forge pyhdf gdal rasterio geopandas pandas numpy tqdm pyyaml -y

装完先验一下，两行命令能省掉后面两小时的困惑：

python -c "from pyhdf.SD import SD, SDC; print('pyhdf ok')" python -c "from osgeo import gdal; print(gdal.__version__); print('HDF4' in gdal.GetDriverByName('HDF4').GetMetadata()['DMD_LONGNAME'])"

2.3 目录骨架与配置：把路径、tile、时间范围赶出代码

代码写到第二个脚本就会开始出现硬编码路径，第三个脚本就会出现“我记得改过啊”的惨案。我的习惯是第一天就把目录和配置定死：

mcd19a2-pipeline/ ├── config.yaml # 路径、tile 列表、时间范围、QA 策略 ├── src/ │ ├── naming.py # 文件名解析、tile 号与仿射变换 │ ├── reader.py # pyhdf 读 SDS、缩放、填充值处理 │ ├── qa.py # QA 位域解析与筛选规则 │ ├── warp.py # 内存数据集构造 + 重投影 │ ├── aggregate.py # 日/月/年合成、站点取值 │ └── cli.py # 命令行入口 ├── data/ │ ├── raw/ # 原始 hdf，按年/日存放 │ ├── interim/ # 单景中间 tif，可随时删 │ └── out/ # 最终产品 ├── tests/ │ └── fixtures/ # 一景已知晴天的样本，用于回归 └── scripts/ └── run_batch.sh

配置文件用 YAML 而不是 Python 常量，理由是二次开发的人往往不是写这套代码的人，改一个 YAML 比读settings.py心理负担小得多：

paths: raw: data/raw interim: data/interim out: data/out study_area: tiles: ["h26v04", "h26v05"] start: "2020-01-01" end: "2020-12-31" bbox: [110.0, 30.0, 120.0, 40.0] # 目标区域，WGS84 warp: dst_crs: "EPSG:4326" resolution: 0.01 resample: bilinear qa: cloud_bits: [0, 2] # 保留 clear(0) 和 probably_clear(2) reject_adjacency: true

这套结构的好处是每一层都能单独跑、单独测。reader.py只负责把 hdf 变成带掩膜的数组，warp.py只负责把数组加上地理参考扔到目标坐标系，中间不互相依赖，出了问题一眼能定位在哪一层。

3. 核心链路拆解：一次 MCD19A2 从 HDF4 到可用 AOD 栅格的完整流转

3.1 读取科学数据集：缩放因子、填充值、数据类型三件套

先看读取这一段。很多人第一次读出来的 AOD 全是 -28 或者 0，就是在这里翻的车。

# src/reader.py import numpy as np from pyhdf.SD import SD, SDC FILL_RAW = -28672 # MCD19A2 AOD 类 SDS 的原始填充值 SCALE = 0.001 # 缩放因子，int16 -> 物理量 def read_sds(hdf_path: str, sds_name: str, scale: float = SCALE): """读取一个 SDS，返回 (float32 数组, 有效掩膜)。 注意顺序：先用原始整型比较填充值，再乘缩放因子。 反过来做会因为浮点精度把 -28672 变成 -28.6720001 之类的值而漏判。 """ hdf = SD(hdf_path, SDC.READ) try: sds = hdf.select(sds_name) raw = sds.get() # 原始整型，别急着 astype attrs = sds.attributes() fill = attrs.get('_FillValue', FILL_RAW) sds.endaccess() finally: hdf.end() raw = raw.astype('int32') # 防止 uint16 相减溢出 valid = raw != fill # 部分 SDS 还有独立的 valid_range，超范围一样当无效 vmin, vmax = attrs.get('valid_range', (None, None)) if vmin is not None: valid &= (raw >= vmin) & (raw <= vmax) data = np.where(valid, raw.astype('float32') * scale, np.nan) return data, valid

逻辑上就三步：先拿到原始整型数组，用_FillValue和valid_range双重判无效，最后才做缩放。为什么非要这个顺序？因为-28672 * 0.001 = -28.672，浮点比较本身没问题，但一旦中间做过别的运算（比如先减add_offset），填充值就可能变成-28.671999，==判断直接失效。血泪经验是：所有掩膜判断都在整型域里完成，进了浮点域就只信np.nan。

参数说明：sds_name传 SDS 全名，注意AngstromExp_470-780这种带连字符的名字在拼 GDAL 子数据集路径时要加引号；scale默认 0.001，但Column_WV之类的 SDS 缩放因子不同，从attrs['scale_factor']读更保险，我这里留了参数让你能覆盖。

3.2 QA 位域筛选：把云、邻域效应和气候态 AOD 挡在门外

AOD_QA是一个 uint16 的位域，不是一个“质量分数”。把它当数值用——比如qa > 0或者qa < 100——是这条链路上最常见的翻车点，因为数值小不代表质量好，位的位置才决定含义。

# src/qa.py import numpy as np def extract_bits(qa: np.ndarray, start: int, length: int) -> np.ndarray: """从 uint16 位域里取出连续 length 位，返回 0~(2^length-1) 的整数。""" mask = (1 << length) - 1 return ((qa.astype('uint16') >> start) & mask).astype('uint8') def aod_qa_mask(qa: np.ndarray, cloud_keep=(0, 3), reject_adjacency: bool = True) -> np.ndarray: """按云状态和邻域标志生成布尔掩膜，True 表示可用。 cloud_keep: 允许的云状态编码。 0 = clear, 1 = possibly cloudy, 2 = cloudy, 3 = probably clear 默认保留 0 和 3，做定量反演时很多人只留 0。 reject_adjacency: 是否剔除受邻域效应污染的像元。 这些像元本身可能标 clear，但周边有云，AOD 会被抬高。 """ cloud = extract_bits(qa, start=0, length=3) keep = np.isin(cloud, cloud_keep) if reject_adjacency: adj = extract_bits(qa, start=5, length=3) # 邻域标志位 keep &= (adj == 0) return keep

位域起始位置和长度做成了参数，不是因为我喜欢抽象，而是因为这个产品的不同版本在位上确实调整过，硬编码进去以后换一版数据就得回来改代码。云状态占低 3 位这点比较稳定，中高位分别是邻域效应、水陆掩膜、气候态 AOD 替代值、耀斑和气溶胶模型，具体顺序以手头版本的产品说明为准，但extract_bits这个工具函数是不变的。

reject_adjacency这个开关值得单独说一句。默认打开会掉 20%~40% 的有效像元，但留下来的值更干净。做城市尺度的 AOD 年均值，我一般打开；做全球尺度的稀疏采样，我会关掉，因为再剔就没数据了。这个取舍没有标准答案，但你必须知道自己按下了哪个开关。

3.3 正弦投影对齐与重投影：warp 之前必须先做的事

这一步是整条链路最容易出错的地方。有两种走法：一是让 GDAL 直接打开HDF4_EOS子数据集然后gdal.Warp，二是自己读数组、手动设地理参考、再 warp。我选第二种，原因是第一种在不同 GDAL 版本上读出来的 geotransform 会有细微差别，做多景叠加时会看到一两个像元的偏移，肉眼看不出来，统计上却是系统性偏差。

手动设参考的关键是瓦片号到仿射变换的换算：

# src/naming.py import re import numpy as np PIXEL = 926.625433055833 # 单像元边长（米） TILE_PX = 1200 # 瓦片像素数 TILE_SPAN = PIXEL * TILE_PX # 1111950.5196669... X_MIN = -20015109.354 # 全球网格最左 Y_MAX = 10007554.677 # 全球网格最上 _HV = re.compile(r'\.h(\d{2})v(\d{2})\.') def parse_tile(hdf_path: str): m = _HV.search(hdf_path) if not m: raise ValueError(f'文件名里找不到瓦片号: {hdf_path}') return int(m.group(1)), int(m.group(2)) def tile_geotransform(h: int, v: int): """返回 GDAL 顺序的 (x0, dx, 0, y0, 0, -dy)。""" x0 = X_MIN + h * TILE_SPAN y0 = Y_MAX - v * TILE_SPAN return (x0, PIXEL, 0.0, y0, 0.0, -PIXEL) SINU_WKT = ( 'PROJCS["MODIS_Sinusoidal",' 'GEOGCS["Unknown",DATUM["Unknown",SPHEROID["Unknown",6371007.181,0]],' 'PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],' 'PROJECTION["Sinusoidal"],PARAMETER["longitude_of_center",0],' 'PARAMETER["false_easting",0],PARAMETER["false_northing",0],' 'UNIT["metre",1]]' )

拿到仿射变换和 WKT 之后，就能把 numpy 数组塞进一个内存数据集，交给 GDAL 做重投影。这样做还有个附带好处：AOD 和 QA 掩膜可以在同一个数组里先合并，重投影只做一次，天然对齐，不会出现“AOD 和 QA 各 warp 一遍结果差半个像元”的问题。

# src/warp.py import numpy as np from osgeo import gdal, osr gdal.UseExceptions() def warp_array(arr: np.ndarray, geotransform, src_wkt: str, dst_path: str, dst_crs: str = 'EPSG:4326', res: float = 0.01, bbox=None, resample: str = 'bilinear') -> str: """把带地理参考的 numpy 数组重投影并写成 GeoTIFF。 arr: float32，无效值用 np.nan bbox: (minx, miny, maxx, maxy)，目标坐标系下的裁剪范围，None 表示不裁 resample: AOD 是连续量，用 bilinear；分类产品才用 near """ ny, nx = arr.shape mem = gdal.GetDriverByName('MEM').Create('', nx, ny, 1, gdal.GDT_Float32) mem.SetGeoTransform(geotransform) srs = osr.SpatialReference() srs.ImportFromWkt(src_wkt) mem.SetProjection(srs.ExportToWkt()) band = mem.GetRasterBand(1) band.SetNoDataValue(float('nan')) band.WriteArray(arr) band.FlushCache() opts = gdal.WarpOptions( format='GTiff', dstSRS=dst_crs, xRes=res, yRes=res, resampleAlg=resample, outputBounds=bbox, targetAlignedPixels=True, # 关键：保证多景输出网格严格一致 multithread=True, creationOptions=['COMPRESS=LZW', 'TILED=YES', 'BIGTIFF=IF_SAFER'], ) gdal.Warp(dst_path, mem, options=opts) mem = None # 及时释放，批处理时很吃内存 return dst_path

targetAlignedPixels=True这一行是整个函数里最重要的一行。没有它，GDAL 会按源数据范围自动算输出网格原点，不同瓦片、不同日期算出来的原点可能差几十米，后面做时序平均的时候就会看到明显的网格错位和边缘锯齿。加上它之后，所有输出都对齐到res的整数倍，叠加才安全。

参数说明：xRes/yRes用 0.01 度大约对应 1km，跟产品原生分辨率相当，再细就是插值放大了；bbox在dst_crs下给，跨投影的 bbox 换算不要自己手算，交给 GDAL 更稳；resample别用near，AOD 虽然是 1km 格网，但重投影后每个输出像元落在不同源像元上，最近邻会引入明显的块状伪影。

3.4 裁剪、导出与元数据落盘

把上面三步串成一个单景处理函数，逻辑就清楚了：

# src/pipeline.py import os import numpy as np from .reader import read_sds from .qa import aod_qa_mask from .naming import parse_tile, tile_geotransform, SINU_WKT from .warp import warp_array def process_one(hdf_path: str, cfg: dict) -> str: aod, _ = read_sds(hdf_path, 'Optical_Depth_047') qa, _ = read_sds(hdf_path, 'AOD_QA', scale=1.0) # QA 不缩放 qa = qa.astype('uint16') mask = aod_qa_mask(qa, cloud_keep=tuple(cfg['qa']['cloud_bits']), reject_adjacency=cfg['qa']['reject_adjacency']) aod = np.where(mask & np.isfinite(aod), aod, np.nan).astype('float32') if not np.isfinite(aod).any(): return '' # 整景无效，直接跳过，省一次 IO h, v = parse_tile(hdf_path) gt = tile_geotransform(h, v) name = os.path.basename(hdf_path).replace('.hdf', '_aod047.tif') dst = os.path.join(cfg['paths']['interim'], name) return warp_array(aod, gt, SINU_WKT, dst, dst_crs=cfg['warp']['dst_crs'], res=cfg['warp']['resolution'], bbox=cfg['study_area']['bbox'], resample=cfg['warp']['resample'])

这里有个容易被忽略的优化：if not np.isfinite(aod).any(): return ''。一整景全被云覆盖的情况在雨季非常常见，提前返回能省掉一次完整的 warp，批处理跑下来能快 15% 左右。同时它也让下游知道“这景没有有效数据”，而不是傻乎乎地写一个全空的 tif 出来占空间。

导出之后记得把处理参数写进 TIFF 的元数据，哪怕只是AOD_SDS=Optical_Depth_047;QA=cloud0,3;adj=1这么一行。三个月后你回头看一个 tif，是不会记得当时用的是哪套 QA 策略的，这个后悔药成本极低。

4. 批处理与时空聚合：把单景脚本改成能跑满一整个 tile 的流水线

4.1 多进程调度与失败重试：几百景文件的排队策略

单景跑通只是开始，真正的工作量是几千景。这里第一个要处理的问题是进程模型。用multiprocessing.Pool直接传 GDAL 数据集对象会炸，因为 GDAL 的 C 句柄不能跨进程序列化；正确的做法是每个 worker 自己打开自己关闭，池子里只传文件路径。

# src/batch.py import os import traceback from concurrent.futures import ProcessPoolExecutor, as_completed from tqdm import tqdm from .pipeline import process_one def _worker(args): """worker 内部独立完成打开-处理-关闭，不共享任何 GDAL 对象。""" hdf_path, cfg = args try: out = process_one(hdf_path, cfg) return hdf_path, out, None except Exception as e: return hdf_path, None, f'{type(e).__name__}: {e}' def run_batch(hdf_files, cfg, workers: int = 8): tasks = [(f, cfg) for f in hdf_files] ok, empty, failed = 0, 0, [] with ProcessPoolExecutor(max_workers=workers) as ex: futures = [ex.submit(_worker, t) for t in tasks] for fut in tqdm(as_completed(futures), total=len(tasks), desc='warp'): path, out, err = fut.result() if err: failed.append((path, err)) elif not out: empty += 1 else: ok += 1 if failed: with open('failed.log', 'w') as f: for p, e in failed: f.write(f'{p}\t{e}\n') print(f'成功 {ok} / 全无效 {empty} / 失败 {len(failed)}') return failed

workers开多少合适？这不是越多越好。GDAL warp 本身是多线程的，单个进程就可能吃掉两三个核，进程数开到CPU 核数 / 2通常比开满更快。我一般先用 8 跑一小批测吞吐，再往上调，直到 IO 成为瓶颈为止——大量小文件场景下，机械硬盘的随机读往往比 CPU 先撑不住。

失败重试的策略要分开看：Exception里如果是文件损坏或者下载不完整，重试一百次也没用，记进failed.log人工处理；如果是临时性的磁盘满或者句柄耗尽，退避重试一次就够。无差别重试只会把真正的错误淹掉，这是我在一个 3000 景的批次上跑了六个小时才发现日志里全是同一个报错的教训。

4.2 日合成、月均值与年均：时间维聚合的两种写法

MCD19A2 本身是日尺度的，但大多数分析要的是月均或者年均。聚合有两种走法，选哪种取决于数据量。

第一种是“边处理边累积”，适合几百景以内的规模：维护一个sum数组和count数组，每处理完一景就累加进去，最后除一次。内存占用是固定的两个float32栅格，非常省。

# src/aggregate.py import glob import numpy as np import rasterio def monthly_mean(tif_glob: str, out_path: str): """同网格的日产品叠加成月均，NaN 不参与计数。""" files = sorted(glob.glob(tif_glob)) if not files: raise FileNotFoundError(tif_glob) acc = None cnt = None profile = None for fp in files: with rasterio.open(fp) as ds: arr = ds.read(1).astype('float64') if acc is None: acc = np.zeros_like(arr) cnt = np.zeros_like(arr) profile = ds.profile.copy() valid = np.isfinite(arr) acc[valid] += arr[valid] cnt[valid] += 1 mean = np.where(cnt > 0, acc / np.maximum(cnt, 1), np.nan).astype('float32') profile.update(dtype='float32', count=1, nodata=np.nan, compress='lzw', tiled=True) with rasterio.open(out_path, 'w', **profile) as dst: dst.write(mean, 1) dst.update_tags(n_files=len(files)) return out_path

用float64累积而不是float32，是因为几百景累加之后float32的有效位会被吃光，尤其是 AOD 值域本来就在 0~2 这种小范围，相对误差会被放大。多占一倍内存换个安心，值。

第二种是“先落盘再算”，适合跨年、跨 tile 的大规模场景：把每景的中间 tif 存下来，再用xarray+dask做延迟计算的resample(time='1M').mean()。这条路的好处是能利用 chunk 并行，坏处是对文件命名和坐标一致性要求极高，文件名里必须能解析出时间，所有 tif 的网格必须完全相同——这也是上一章targetAlignedPixels那么重要的原因。

4.3 站点/区域取值的接口设计与 CSV 输出

做地面验证或者跟站点观测对比的时候，需要把栅格值取到点位上。这一步的接口设计我想多说两句，因为它决定了这个工具能不能被别人用起来。

# src/sample.py import numpy as np import pandas as pd import rasterio def sample_points(tif_path: str, lons, lats) -> np.ndarray: """按经纬度双线性采样，返回一维数组，落在无效区返回 NaN。 用 rasterio.transform.rowcol 拿到浮点行列号，自己做双线性， 比 nearest 稳，也不会像 GDAL 的 gdal.LocationInfo 那样受无效值影响。 """ with rasterio.open(tif_path) as ds: arr = ds.read(1).astype('float64') inv = ~ds.transform # 逆变换：地理坐标 -> 像素坐标 out = np.full(len(lons), np.nan, dtype='float32') for i, (lon, lat) in enumerate(zip(lons, lats)): col_f, row_f = inv * (lon, lat) r, c = int(np.floor(row_f)), int(np.floor(col_f)) if r < 0 or c < 0 or r + 1 >= arr.shape[0] or c + 1 >= arr.shape[1]: continue dr, dc = row_f - r, col_f - c block = arr[r:r + 2, c:c + 2] if not np.isfinite(block).all(): # 四邻域里只要有 NaN，就退化成最近邻，避免污染整片 out[i] = arr[r, c] continue out[i] = (block[0, 0] * (1 - dr) * (1 - dc) + block[0, 1] * (1 - dr) * dc + block[1, 0] * dr * (1 - dc) + block[1, 1] * dr * dc) return out def batch_sample(tif_list, stations: pd.DataFrame, out_csv: str): rows = [] for fp in tif_list: vals = sample_points(fp, stations['lon'].values, stations['lat'].values) date = fp.split('_')[-2] # 视你的命名规则而定 for sid, v in zip(stations['id'], vals): rows.append({'station': sid, 'date': date, 'aod': v}) pd.DataFrame(rows).to_csv(out_csv, index=False)

双线性这里有个细节要注意：四个角里只要有一个是 NaN 就退回最近邻。如果不这么处理，直接用np.nanmean或者忽略 NaN 做加权，边缘像元的值会被高估——AOD 的边缘往往是云边界，本来值就偏高。这个处理方式算不上完美，但比“整片区域全丢”和“边缘值失真”都强。

输出 CSV 的字段设计建议留station、date、aod，再加一列valid_pixels记录该点位附近有多少有效像元，方便后面做质量过滤。少了这一列，别人用你的数据时没法判断某个点是真实观测还是边缘插值出来的。

5. 避坑与排查：QA 位域、投影对齐和填充值上的五处常见翻车

5.1 现象：读出来的 AOD 全是 -28 或者 0

原因：跳过了缩放和填充值处理，把原始int16直接当物理量用了。填充值 -28672 乘以 0.001 就是 -28.672，看起来像“异常低的 AOD”。

解决：所有 SDS 读取统一走read_sds，在整型域里做valid判断，最后才乘scale。写完第一版之后，找一景已知晴天的数据，打印一下有效像元比例——陆地上通常在 20%~60% 之间，如果出来 99% 或者 0.5%，基本可以确定掩膜写错了。

5.2 现象：重投影后 AOD 和 QA 差半个像元，叠加出条纹

原因：AOD 和 QA 分别 warp 了两次，两次的源数据集虽然一样，但 GDAL 在做重采样时的边界处理有细微差异；更常见的情况是两次 warp 之间网格原点算得不一样。

解决：所有需要空间对齐的图层，在源坐标系下先合并成一个多波段数组，只做一次 warp。如果业务上确实要分开出图，就强制给两次 warp 传同一个outputBounds和同一个xRes/yRes，并打开targetAlignedPixels。

5.3 现象：天气很好的区域也是大片空值

原因：QA 位域取错。把AOD_QA当数值比较，比如qa < 8想筛出晴天，结果把大量probably clear（编码 3）之外的像元全删了；另一种情况是位移写错，start从 1 开始数而不是 0。

解决：先用extract_bits把云状态单独导出来做个直方图，看看 0/1/2/3 四类的占比。正常的白天场景，0 和 3 加起来应该占一半以上。如果某一类占 99%，那就是位移错了。

5.4 现象：月度合成结果边缘有锯齿、和原始单景对不上

原因：不同日期的 tif 网格不一致。GDAL 默认会自动计算输出范围和原点，瓦片、日期一变，原点就变，叠加时半像素错位累积成锯齿。

解决：全流程固定一个目标网格，用outputBounds显式指定范围，targetAlignedPixels=True强制对齐，聚合前再校验一遍所有输入的transform和shape是否完全一致，不一致直接报错而不是静默平均。

5.5 现象：多进程跑到一半卡死，或者内存一路涨到被系统杀掉

原因：GDAL 数据集对象泄漏，mem = gdal.GetDriverByName('MEM').Create()之后没释放；或者pyhdf的SD对象没end()。这些句柄不会自动回收，跑几百景之后内存就爆了。

解决：所有SD用try/finally保证end()，所有MEM数据集用完置None。另外把workers从核数降到一半试试，很多时候“卡死”只是 IO 排队排到了天亮。

6. 二次开发：把这条链路封装成可配置插件与新波段的扩展点

跑到这里，单景、批处理、聚合、取点都能用了，剩下的工作往往是“给别人用”。这时候要不要做成插件，取决于使用者是谁。如果对方是写 Python 的同事，click或者typer注册几个子命令就够了；如果对方只认桌面 GIS，那就得挂进 QGIS 的 Processing 框架里，做成一个能被工具箱搜索到的算法。

CLI 的封装成本最低，pyproject.toml里注册一个入口点：

[project.scripts] mcd19 = "mcd19a2_pipeline.cli:main"

装完之后用户直接mcd19 warp --config config.yaml就能跑。子命令按阶段拆成read、warp、monthly、sample四个，每个都能单独执行，这比一个大而全的run-all实用得多——实际使用中，八成的人只需要其中一步。

QGIS 插件那条路要继承QgsProcessingAlgorithm，核心是把process_one包进processAlgorithm里，参数用addParameter声明，进度用feedback.setProgress汇报。这里有个细节：QGIS 的插件运行在它自己的 Python 环境里，pyhdf得装到那个环境里，而不是你原来的 conda 环境。这一点不提前说清楚，交付当天一定会被问“为什么点了没反应”。

扩展点上，我整理了一个表：新加一个 SDS 只需要在配置里加名字、在process_one里多读一次，风险最低；换重投影目标要改dst_crs和bbox，注意不同坐标系下res的含义不一样，投影坐标系里是米，地理坐标系里是度；自定义 QA 策略风险最高，因为aod_qa_mask的语义变了，下游所有统计都得重新校准。

验证方法上，我强烈建议保留一景“黄金样本”：挑一个天气干净、AOD 空间分布均匀的日子，把当前版本的输出统计量（有效像元数、均值、P10/P90）写进tests/，任何一次改动之后跑一遍对比，偏差超过 1% 就得停下来查。这比任何单元测试都管用，因为它验的是整条链路而不只是某个函数。

从那以后我每加一个新 SDS、每改一次 QA 规则，都强制先跑一遍黄金样本对比，再放进批处理队列。这个习惯让我躲过了至少三次“看起来没问题、跑到一半发现全错”的返工。希望帮到你。

本文还有配套的精品资源，点击获取