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常用的几个SDS 名称含义数据类型备注Optical_Depth_047470nm 气溶胶光学厚度int16缩放因子 0.001Optical_Depth_055550nm 气溶胶光学厚度int16缩放因子 0.001AOD_QA质量控制位域uint16按位解析不是数值大小AOD_UncertaintyAOD 不确定度int16同样带缩放AngstromExp_470-780angstrom 指数int16名字里带连字符拼接路径要小心Column_WV大气可降水量int16常被顺手拿来做水汽分析网格约定是另一个必须记住的常量全球正弦网格被切成 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 版本上偶有半像素级的差异做长时序叠加时这就是条带。环境搭建我一般用 condapyhdf用 pip 装经常卡在编译 HDF4 上conda-forge 有现成二进制包conda create -n mcd19 python3.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, start0, length3) keep np.isin(cloud, cloud_keep) if reject_adjacency: adj extract_bits(qa, start5, length3) # 邻域标志位 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, bboxNone, 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( formatGTiff, dstSRSdst_crs, xResres, yResres, resampleAlgresample, outputBoundsbbox, targetAlignedPixelsTrue, # 关键保证多景输出网格严格一致 multithreadTrue, creationOptions[COMPRESSLZW, TILEDYES, BIGTIFFIF_SAFER], ) gdal.Warp(dst_path, mem, optionsopts) mem None # 及时释放批处理时很吃内存 return dst_pathtargetAlignedPixelsTrue这一行是整个函数里最重要的一行。没有它GDAL 会按源数据范围自动算输出网格原点不同瓦片、不同日期算出来的原点可能差几十米后面做时序平均的时候就会看到明显的网格错位和边缘锯齿。加上它之后所有输出都对齐到res的整数倍叠加才安全。参数说明xRes/yRes用 0.01 度大约对应 1km跟产品原生分辨率相当再细就是插值放大了bbox在dst_crs下给跨投影的 bbox 换算不要自己手算交给 GDAL 更稳resample别用nearAOD 虽然是 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, scale1.0) # QA 不缩放 qa qa.astype(uint16) mask aod_qa_mask(qa, cloud_keeptuple(cfg[qa][cloud_bits]), reject_adjacencycfg[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_crscfg[warp][dst_crs], rescfg[warp][resolution], bboxcfg[study_area][bbox], resamplecfg[warp][resample])这里有个容易被忽略的优化if not np.isfinite(aod).any(): return 。一整景全被云覆盖的情况在雨季非常常见提前返回能省掉一次完整的 warp批处理跑下来能快 15% 左右。同时它也让下游知道“这景没有有效数据”而不是傻乎乎地写一个全空的 tif 出来占空间。导出之后记得把处理参数写进 TIFF 的元数据哪怕只是AOD_SDSOptical_Depth_047;QAcloud0,3;adj1这么一行。三个月后你回头看一个 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_workersworkers) as ex: futures [ex.submit(_worker, t) for t in tasks] for fut in tqdm(as_completed(futures), totallen(tasks), descwarp): 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 failedworkers开多少合适这不是越多越好。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(dtypefloat32, count1, nodatanp.nan, compresslzw, tiledTrue) with rasterio.open(out_path, w, **profile) as dst: dst.write(mean, 1) dst.update_tags(n_fileslen(files)) return out_path用float64累积而不是float32是因为几百景累加之后float32的有效位会被吃光尤其是 AOD 值域本来就在 0~2 这种小范围相对误差会被放大。多占一倍内存换个安心值。第二种是“先落盘再算”适合跨年、跨 tile 的大规模场景把每景的中间 tif 存下来再用xarraydask做延迟计算的resample(time1M).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, dtypefloat32) 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, indexFalse)双线性这里有个细节要注意四个角里只要有一个是 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显式指定范围targetAlignedPixelsTrue强制对齐聚合前再校验一遍所有输入的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的语义变了下游所有统计都得重新校准。扩展点要改的地方风险增加 SDS如Column_WV配置文件 process_one多读一次低注意缩放因子不同更换重投影目标config.yaml的dst_crs/bbox/resolution中res单位随坐标系变化自定义 QA 位域规则qa.py的aod_qa_mask高会改变有效像元总体换聚合时间粒度aggregate.py的文件筛选逻辑低前提是网格已对齐接入新数据源naming.py的瓦片解析与仿射变换中投影参数要重算验证方法上我强烈建议保留一景“黄金样本”挑一个天气干净、AOD 空间分布均匀的日子把当前版本的输出统计量有效像元数、均值、P10/P90写进tests/任何一次改动之后跑一遍对比偏差超过 1% 就得停下来查。这比任何单元测试都管用因为它验的是整条链路而不只是某个函数。从那以后我每加一个新 SDS、每改一次 QA 规则都强制先跑一遍黄金样本对比再放进批处理队列。这个习惯让我躲过了至少三次“看起来没问题、跑到一半发现全错”的返工。希望帮到你。本文还有配套的精品资源点击获取