简介：这份代码包围绕MCD19A2遥感数据处理展开，面向具备一定IDL与Python基础、从事地球科学、环境监测及植被分析的科研人员与学习者，重点解决批量预处理、影像融合、镶嵌、重采样、裁剪和栅格点值提取等流程的代码实现问题。压缩包共12个文件，约11KB，以3个py脚本和3个pro程序为主体，分别承担重投影、裁剪、AOD栅格提取点值以及MCTK调用、均值统计、批量镶嵌等任务，另含xml、iml等IDE项目配置，便于在对应环境中查看与复用。代码覆盖MCD19A2从影像组织到专题值提取的关键环节，可帮助读者理解MCTK批量处理的脚本结构、参数设置与流程衔接，并据此改造自己的自动化处理方案，对刚接触该产品、需要快速搭建批量流程的读者尤其友好。已有1248人学习，适合作为遥感数据批量处理的学习参考与代码模板。

1. 为什么一份气溶胶数据的处理代码值得单独拆开看

第一次拿到 MCD19A2 的人，八成会卡在"打开文件"这一步：HDF4 容器、二维正弦投影、带比例因子的整型子数据集、藏在一串二进制位里的质量标记，四件事叠在一起，能直接跑通的完整代码反而不多。这份资源就是围绕 MCD19A2 数据处理过程整理的一组代码，覆盖读取、质量筛选、重投影、裁剪镶嵌到时间聚合的整条链路。它不是那种装完就能点的软件/插件，而是一批可读写、可改参数、可挂到批处理任务上的脚本，适合要批量出 AOD 栅格的研究生、做大气环境分析的工程师，以及手头有一堆 HDF 文件却不知道怎么变成 GeoTIFF 的人。

2. 环境与数据准备：先把 HDF4 容器和正弦网格摸清楚

2.1 依赖装不上的那半天：pyhdf 与 GDAL 的编译依赖

这套代码要跑起来，绕不开两个"重量级"依赖：读 HDF4 需要pyhdf，重投影和写出需要gdal。这两个包最恶心的地方在于，它们都不是纯 Python，底层各自带着 C 库。用pip install pyhdf直接装，在没有预编译轮子的机器上大概率会卡在编译阶段报错，提示找不到hdf.h。血泪经验是先装系统的 HDF4 开发库，再装 Python 包。

# 用 conda 装是最省事的路径，它会自动把底层 C 库一起解决 conda create -n maiac python=3.10 conda activate maiac conda install -c conda-forge pyhdf gdal rasterio pyproj xarray numpy -y # 装完先验一下 GDAL 有没有把 HDF4 驱动编进去 python -c "from osgeo import gdal; print(gdal.GetDriverByName('HDF4') is not None)"

上面这段命令做了三件事：建环境、装依赖、验证 GDAL 是否带 HDF4 驱动。最后一行打印True才算过关，打印False说明这个 GDAL 是精简编译的，读不了 MCD19A2，得换一个 conda-forge 的包。这个检查我一般放在所有脚本之前，因为驱动缺失的表现是"文件明明存在，gdal.Open却返回 None"，排查起来特别费时间。

第二条命令里pyhdf和gdal装在同一个 conda 环境里，是为了避免两边链接到不同版本的 HDF4 库导致段错误。如果团队里有人用 pip 混装，建议让所有人都用同一份环境文件，别各装各的。

2.2 一个文件里到底装了什么：SDS 清单与属性字段

MCD19A2 是网格化的二级产品，每个 HDF4 文件里装着十几个子数据集（SDS），常用的其实就三组：气溶胶光学厚度（Optical_Depth_047、Optical_Depth_055）、质量标记（AOD_QA）、观测几何（Solar_Zenith、View_Zenith、Relative_Azimuth）。这些 SDS 全是整型，真实物理量要靠scale_factor和add_offset反算，无效像元用_FillValue填充。

在动手写正式代码前，先跑一段探测脚本把文件结构打出来，比翻文档快得多：

from pyhdf.SD import SD, SDC path = "MCD19A2.A2020150.h26v05.006.2020158031655.hdf" hdf = SD(path, SDC.READ) # datasets() 返回 {SDS 名称: (维度, 类型, 维度索引)} for name, (dims, dtype, idx) in hdf.datasets().items(): sds = hdf.select(name) attrs = sds.attributes() print(f"{name:24s} dims={dims} dtype={dtype}") # 只打印和数值换算相关的关键属性，属性全打会刷屏 for key in ("scale_factor", "add_offset", "_FillValue", "valid_range"): if key in attrs: print(f" {key} = {attrs[key]}") sds.endaccess() hdf.end()

第一个datasets()给出的是全部 SDS 的元信息，维度一般是(1200, 1200)，也就是单块瓦片 1200 行 1200 列。注意dtype字段，气溶胶厚度通常是int16，QA 可能是 8 位也可能是 16 位，不同批次的数据集版本会不一样，所以代码里不要写死位宽，先打印出来再决定掩码怎么写。

第二段attributes()是重点。scale_factor决定了原始整数乘多少才是物理量，气溶胶一般都是 0.001，也就是原始值 1000 对应 AOD 1.0；_FillValue常见是 -28672，所有等于这个值的像元都必须先剔除再做运算，否则参与平均会把结果拉得一塌糊涂。valid_range则是产品自己声明的有效区间，AOD 允许轻微负值（代表低于探测限的洁净大气），所以别一上来就把负数全砍掉。

还有一点：每次select之后必须endaccess()，文件句柄会累积，批量处理几千个文件时不释放很快就会撞上系统句柄上限。这个细节在单文件测试里完全看不出来，一上批量就翻车。

2.3 目录约定与配置入口：让批处理脚本认识你的数据盘

原始数据一般按"年/日序"或者"瓦片号"分层存放，脚本如果每换一批数据就改一次路径，维护成本会越来越高。比较常见的做法是抽一个 JSON 配置文件出来，把所有可变的东西集中到一个地方：

{ "input_dir": "/data/maiac/2020", "output_dir": "/data/out/aod_550", "boundary_shp": "/data/shp/region.shp", "dataset_name": "Optical_Depth_055", "qa_datasets": ["AOD_QA"], "qa_accept_bits": [0, 1], "scale": 0.001, "fill_value": -28672, "valid_min": -100, "valid_max": 5000, "out_nodata": -9999.0, "resampling": "bilinear", "target_epsg": 4326, "max_workers": 8, "resume_log": "/data/out/processed.jsonl" }

字段分成三类看：input_dir/output_dir/boundary_shp是路径；dataset_name、scale、fill_value、valid_*是数据解码参数，换了产品版本或者换了波段就要动这里；resampling、target_epsg、max_workers是处理策略。把这些全部外置的好处是，同一份代码既能出 550nm 也能出 470nm，改一个字段就行，不用去翻源码。

配置里特意留了resume_log，它是断点续跑的依据，后面第 4 章会用到。目录上建议输出按"年/月"再分一层，单目录塞几万个 GeoTIFF，Windows 上资源管理器会直接卡死，Linux 上ls也要等半天。

3. 核心读取链路：从 HDF4 子数据集到带地理信息的数组

3.1 读一个 SDS 的完整动作：填充值、比例因子与无效值

读取这一步看起来简单，实际上顺序错了结果就全废。正确的顺序是：先取原始整型数组 → 剔除填充值 → 转浮点 → 乘比例因子 → 用valid_range之外的再筛一遍。反过来先乘比例因子再比较填充值，会因为浮点精度问题漏掉一部分无效像元。

import numpy as np from pyhdf.SD import SD, SDC def read_sds(path, name, scale=0.001, fill=-28672, vmin=-100, vmax=5000): """读取单个 SDS，返回 float32 数组，无效像元为 NaN。""" hdf = SD(path, SDC.READ) try: sds = hdf.select(name) raw = sds.get() # 原始整型，形状 (1200, 1200) attrs = sds.attributes() # 属性里有的以属性为准，没有就用传入的默认值 scale = attrs.get("scale_factor", scale) fill = attrs.get("_FillValue", fill) sds.endaccess() finally: hdf.end() arr = raw.astype(np.float32) # 第一步：填充值置 NaN，必须在乘比例之前做 arr[raw == fill] = np.nan # 第二步：整型原始值按 valid_range 过滤，超出即无效 arr[(raw < vmin) | (raw > vmax)] = np.nan # 第三步：换算成物理量 arr = arr * float(scale) return arr

read_sds里有几个点值得说。sds.get()拿到的是完整的二维数组，1200×1200 的int16大约 2.8 MB，一次读进内存没问题；但如果是镶嵌后的更大范围，就得考虑分块读，pyhdf支持按行列切片读取，可以后面再优化。

attrs.get(...)这种写法是给不同批次数据留后路。同一个产品不同版本，scale_factor偶尔会变，写死之后换数据就静默出错——数值看着正常，量级差 10 倍，这种错误最难发现，因为图还是能画出来的。用属性优先、参数兜底的方式，两边的坑都能避开。

第三步的量纲换算放在最后，是为了让前两步的比较都在整型域完成，语义最清晰。返回float32而不是float64，是因为后面要做多年合成，内存占用差一倍，float32对 AOD 这种精度完全够用。

还有个小细节：try/finally保证异常时也关闭文件。批处理跑几千个文件，中间只要有一个损坏文件抛异常没关句柄，跑到后面就会莫名其妙地失败。

3.2 AOD_QA 位运算：把一团数字拆成三档可信度

AOD 的 QA 不是"好/坏"两档，而是一串二进制位，不同位段分别表示云掩膜邻域、地表类型、反演质量等级、耀斑影响等。直接把 QA 当数值比大小（比如qa > 2）是典型误用，因为数值大小和位段含义完全不对应。

稳妥的做法是用位掩码逐段提取：

import numpy as np def decode_aod_qa(qa, quality_shift=5, quality_mask=0b111, accept_levels=(2, 3)): """从 AOD_QA 中提取反演质量等级，返回布尔掩膜。""" qa = qa.astype(np.uint16) # 右移到目标位段起始位置，再与掩码取出该段的原始值 level = (qa >> quality_shift) & quality_mask # 只保留质量等级在可接受集合内的像元 good = np.isin(level, accept_levels) # 顺带把未反演（level==0）的像元单独标出来，便于统计 return good, level

quality_shift和quality_mask是可调参数，不同批次产品对 QA 位段的定义可能不一样，改之前一定要对着产品文档的位定义表核一遍，别照抄别人博客里的数字。这也是我把 shift 和 mask 做成函数参数而不是写死常量的原因。

np.isin那一步等价于"只接受质量等级为 2 和 3 的像元"，如果你想更宽松，把accept_levels改成(1, 2, 3)即可。返回的level数组别扔掉，后面做质量统计、画 QA 分布图都用得上，比重新解一遍省事。

实际用的时候，通常会再叠一层几何筛选：太阳天顶角过大（低太阳高度角）时反演误差明显上升，加一条solar_zenith < 70的条件能砍掉不少边缘数据。几何 SDS 的比例因子一般是 0.01，读法跟 AOD 一样。

最终掩膜是"质量合格 & 几何合格 & 非 NaN"三者取与。这一步做完，一整块瓦片大概只剩 30%～60% 的有效像元，看着心疼，但这是让结果能用的前提。

3.3 正弦投影到 WGS84：自己算仿射参数还是交给 GDAL

MCD19A2 用的是正弦投影（Sinusoidal），坐标单位是米，原点在赤道与本初子午线交点。每个瓦片覆盖固定的地理范围，只要知道瓦片号和网格尺寸，就能手算出左上角坐标和像元大小：

import numpy as np R = 6371007.181 # 正弦投影球体半径（米） TILE = 1111950.519667 # 单块瓦片边长（米） GRID = 1200 # 每块瓦片行列数 PIXEL = TILE / GRID # 像元边长，约 926.625 米 X_ORIGIN = -20015109.354 # 全球网格最左边界 Y_ORIGIN = 10007554.677 # 全球网格最上边界 def tile_transform(h, v): """由瓦片行列号推出 GDAL 风格的仿射六参数。""" x_min = X_ORIGIN + h * TILE y_max = Y_ORIGIN - v * TILE return (PIXEL, 0.0, x_min, 0.0, -PIXEL, y_max)

tile_transform返回的六个数字就是 GDAL 的 GeoTransform：(x方向像元宽, 旋转, 左上角x, 旋转, y方向像元高(负), 左上角y)。h和v从文件名里解析出来，第 4 章会给正则。注意 y 方向的像元高是负数，因为栅格行号向下增长而地理坐标向上增长，符号写错的结果就是图像上下颠倒，但看上去"还挺像那么回事"，属于隐蔽性很强的错误。

手算仿射参数的好处是不依赖外部库，速度快，适合先落一个中间 GeoTIFF 再统一重投影的两段式流程。缺点是只对标准网格有效，如果数据被裁过或者重采样过，公式就对不上了。

另一种做法是直接用 GDAL 一步到位：

from osgeo import gdal, osr import numpy as np def warp_to_wgs84(src_path, dst_path, arr, transform, epsg=4326): """把内存数组写成带正弦投影的临时栅格，再重投影到 WGS84。""" driver = gdal.GetDriverByName("GTiff") tmp = "/vsimem/tmp_sinu.tif" # 内存文件，不落磁盘 ds = driver.Create(tmp, arr.shape[1], arr.shape[0], 1, gdal.GDT_Float32) ds.SetGeoTransform(transform) srs = osr.SpatialReference() srs.ImportFromProj4( f"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a={R} +b={R} +units=m +no_defs") ds.SetProjection(srs.ExportToWkt()) band = ds.GetRasterBand(1) band.SetNoDataValue(np.nan) band.WriteArray(arr) band.FlushCache() ds = None gdal.Warp(dst_path, tmp, dstSRS=f"EPSG:{epsg}", resampleAlg="bilinear", srcNodata=np.nan, dstNodata=np.nan, multithread=True)

/vsimem/是 GDAL 的内存文件系统，临时栅格不落磁盘，批量处理时能省掉大量 IO。SetNoDataValue(np.nan)配合srcNodata/dstNodata保证无效像元在重采样时不会污染邻域——这一点很关键，双线性插值会把 NaN 周围的有效值一起拉低。

multithread=True在大范围镶嵌时提升明显，但单块 1200×1200 的瓦片上收益有限，可以按需关掉。

两段式（先手算仿射写中间文件，再gdal.Warp）和一步式（内存数组直接 warp）各有场景：前者适合需要保留中间结果的调试，后者适合流水线里只想拿到最终产物。我在批量任务里用的是后者，因为磁盘上多出来的中间文件清理起来很烦。

4. 批处理流水线：把单文件脚本变成能跑一整年的任务

4.1 从文件名解析瓦片号与日期

MCD19A2 的文件名是高度结构化的，年份、日序、瓦片行列号、版本、处理时间全在里面。与其一个个手动填参数，不如用正则一次性抽出来：

import re from datetime import datetime, timedelta PATTERN = re.compile( r"^MCD19A2\.A(?P<year>\d{4})(?P<doy>\d{3})" r"\.h(?P<h>\d{2})v(?P<v>\d{2})", re.IGNORECASE) def parse_name(path): """从文件名解析日期与瓦片号，解析失败返回 None。""" name = path.rsplit("/", 1)[-1] m = PATTERN.match(name) if not m: return None year, doy = int(m.group("year")), int(m.group("doy")) date = datetime(year, 1, 1) + timedelta(days=doy - 1) return { "date": date, "year": year, "month": date.month, "h": int(m.group("h")), "v": int(m.group("v")), "tile": f"h{m.group('h')}v{m.group('v')}", }

日序转日期那行要注意：datetime(year, 1, 1) + timedelta(days=doy - 1)，减 1 是因为 1 月 1 日对应日序 001。忘了减 1 会导致全年日期整体偏移一天，做时间序列分析时和地面站点对不上，查起来得从头翻。

正则只匹配到瓦片号就停，不去匹配版本和处理时间，因为这两个字段在后续流程里用不上，匹配得越长越容易因为格式差异失配。返回None而不是抛异常，是为了让批处理能跳过不符合命名的杂文件（比如产品附带的元数据文件）。

解析出来的tile直接用于 3.3 的tile_transform，date用于输出路径分年和时间聚合的分组键。

4.2 多进程 + 断点续跑：中断之后不用从头再来

单文件处理大概要几百毫秒到一两秒，一年的全国范围数据动辄几千个文件，串行跑要几个小时。用进程池并发是常规操作，但真正让批处理"能长期用"的是断点续跑——跑到一半机器重启、磁盘满、被踢下线，重新跑的时候不用从头开始。

import json, os from concurrent.futures import ProcessPoolExecutor, as_completed def load_done(log_path): """读取已完成记录，返回已处理文件路径集合。""" done = set() if os.path.exists(log_path): with open(log_path, "r", encoding="utf-8") as f: for line in f: line = line.strip() if not line: continue try: done.add(json.loads(line)["src"]) except (json.JSONDecodeError, KeyError): continue # 中断时最后一行可能是半截，直接跳过 return done def run_batch(files, worker, log_path, max_workers=8): done = load_done(log_path) todo = [f for f in files if f not in done] print(f"total={len(files)} done={len(done)} todo={len(todo)}") with ProcessPoolExecutor(max_workers=max_workers) as pool, \ open(log_path, "a", encoding="utf-8") as log: futures = {pool.submit(worker, f): f for f in todo} for fut in as_completed(futures): src = futures[fut] try: out = fut.result() log.write(json.dumps({"src": src, "out": out}) + "\n") log.flush() # 每条立刻落盘，进程被杀也不丢 except Exception as e: log.write(json.dumps({"src": src, "error": str(e)}) + "\n") log.flush()

load_done里对 JSON 解析失败的行直接跳过，是因为进程被强杀时，最后一行大概率只写了半截。如果这里抛异常，整批任务就再也启动不了了，非常要命。

log.flush()那两行是整个断点续跑的关键。Python 的文件缓冲默认几 KB，不主动 flush 的话，程序被kill -9时缓冲区里的记录全丢，重跑时又得把已经做完的文件再做一遍。加上 flush 之后，最多损失正在处理的那几个文件。

max_workers不要设得太大。每个进程要加载 GDAL 和 numpy，内存占用不低，8 个进程大概吃掉 3～4 GB，普通工作站开 16 个很容易触发内存交换，速度反而下降。一般取 CPU 物理核数的一半左右比较稳。

日志用的是 JSON Lines 格式，一行一条记录，追加写，不用读全文件就能续跑，比 CSV 更适合这种场景。

4.3 镶嵌、按边界裁剪与时间维度聚合

单景处理完之后，常见需求有三个：把相邻瓦片拼成一整块、按研究区边界裁出来、按时间维度做合成。这三件事各自有坑。

镶嵌用 GDAL 的Warp一次性传多个输入最省事：

from osgeo import gdal def mosaic_and_clip(inputs, out_path, shp_path, epsg=4326): """多景镶嵌后用边界裁切，一步出结果。""" options = gdal.WarpOptions( format="GTiff", dstSRS=f"EPSG:{epsg}", cutlineDSName=shp_path, # 直接按矢量边界裁 cropToCutline=True, # 输出范围紧贴边界，不留空白 resampleAlg="bilinear", srcNodata=float("nan"), dstNodata=-9999.0, multithread=True, creationOptions=["COMPRESS=LZW", "TILED=YES", "BIGTIFF=IF_SAFER"], ) gdal.Warp(out_path, inputs, options=options)

cutlineDSName直接吃矢量文件，省掉"先镶嵌再裁剪"两步。cropToCutline=True让输出范围紧贴边界外接框，不保留整景范围，文件能小很多。

BIGTIFF=IF_SAFER是给大范围数据留的后路，超过 4 GB 的 GeoTIFF 不加这个参数会直接写失败，报错信息还很不直观。COMPRESS=LZW对浮点栅格压缩率一般，但聊胜于无。

时间聚合建议用xarray，它对带时间维的数据立方体处理起来比手写循环清晰得多：

import xarray as xr import glob def monthly_composite(pattern, out_nc): """按天文件聚合成月均值，输出 NetCDF。""" files = sorted(glob.glob(pattern)) ds = xr.open_mfdataset( files, combine="by_coords", preprocess=lambda d: d.where(d > -9000), # 剔除 nodata parallel=True, ) monthly = ds.resample(time="1MS").mean(skipna=True) monthly.to_netcdf(out_nc, encoding={"__xarray_dataarray_variable__": {"zlib": True, "complevel": 4}})

preprocess里那行where(d > -9000)是专门对付 nodata 的。如果输入文件的 nodata 设成了 -9999 而xarray没识别出来，月均值会被拉低一大截，图上看就是大片异常的暗区。这一步加上去之后就干净了。

resample(time="1MS")表示按月起算，skipna=True保证有效像元少的日子不会把整月拖成 NaN。编码里的complevel=4是在压缩率和速度之间取的折中，设到 9 会慢好几倍。

5. 避坑记录：五个让我重跑整夜的排查现场

5.1 读出来的数组全是 -28672

现象：脚本没报错，输出的 GeoTIFF 打开后一片均匀值，统计一下全是 -28672 或者 0。

原因：sds.get()返回的原始数组里，无效像元就是-28672，代码里只乘了比例因子没有先剔除填充值，或者剔除用了np.isnan(raw)——原始数组是整型，isnan判断永远返回 False。

解决：先确认原始数组的 dtype，整型的话用raw == fill_value做布尔索引，不要用isnan。填充值从_FillValue属性里取，取不到再用默认值。这个坑我在第一次跑批的时候踩过，一整晚的机时全浪费在生成垃圾文件上，从那以后每个新数据集我都先跑一遍小样本统计。

5.2 重投影后图像整体偏移了半个像元

现象：把结果和研究区矢量叠在一起，图像整体往左上角偏了差不多一个像元，边缘尤其明显。

原因：犯的是 GeoTransform 的经典错误——把像元中心坐标当成了像元左上角坐标，或者把正弦投影的球体半径写成了 WGS84 的椭球参数。MCD19A2 用的是球体，a和b都应该是 6371007.181，混用椭球参数会让整幅图产生系统性位移。

解决：检查ImportFromProj4里的+a和+b是否都写成了球半径，检查 GeoTransform 的第三个和第六个参数是不是左上角外边界。验证方法很简单：取一个已知经纬度的地面站点，看它在栅格上的位置和图上的位置是否一致，偏半个像元在图上肉眼就能看出来。

5.3 多进程跑到一半内存爆掉

现象：批处理跑到几百个文件后系统开始疯狂交换，速度骤降，最后被 OOM Killer 干掉。

原因：一是max_workers开太大，每个进程都加载了完整 GDAL 和 numpy，基础内存就有几百 MB；二是sds.endaccess()或者ds = None漏了，GDAL 数据集对象没释放，内存只增不减；三是用open_mfdataset一次性打开了上千个文件，每个都保持句柄。

解决：把max_workers降到物理核数的一半，用psutil监控每个进程的 RSS；所有 HDF 和 GDAL 对象都放在try/finally里或者显式置 None；open_mfdataset加chunks={"time": 30}做分块读取，让xarray按需加载而不是全读进内存。

5.4 镶嵌结果出现明显的接缝

现象：相邻两景拼起来之后，交界处有一条颜色突变的线，AOD 值在接缝两侧差 0.1 以上。

原因：接缝处的像元在两景里都是边缘像元，重采样时双线性插值把 nodata 混了进来，导致边缘值被拉低。另一种可能是两景的观测时间不同，气溶胶本身就在变化，属于真实差异而不是处理错误。

解决：先确认是不是 nodata 污染，把srcNodata正确设上，或者在镶嵌前对每景做一次边缘腐蚀（把最外圈几个像元置为无效）。如果确认是时间差异，就改成分时间段镶嵌，或者接受这个差异并在论文里说明。我一般会做一个接缝剖面图，沿垂直于接缝的方向采样，一眼就能看出是插值问题还是真实变化。

5.5 时间聚合后数值整体偏低

现象：月均值和单日值对比，均值明显偏低，而且有效像元比例越高，偏低越明显。

原因：xarray在做mean时如果没有正确识别 nodata，会把-9999当成真实值参与平均；另一种情况是老版本open_mfdataset的combine="by_coords"在时间坐标有微小差异时产生了重复时间点，导致权重错乱。

解决：在preprocess里显式把 nodata 转成 NaN，并且用ds.where(ds != nodata)而不是ds.where(ds > -9000)这种模糊判断；聚合前打印一下ds.time的长度和文件数是否一致，不一致说明合并出了问题。还有一个习惯：聚合完成后随机抽一个像元，把它所有日的原始值和均值手动算一遍对一下，这招比看统计量管用得多。

6. 二次开发：把脚本改成一个能长期维护的处理库

脚本跑通只是起点，真正要长期用，得把它拆成"配置 + 处理类 + 命令行入口"三层。我一般的做法是定义一个AodProcessor类，把解码参数和处理策略都在__init__里接住，process_one只负责单个文件，run负责调度。这样换产品、换波段、换研究区都只动配置，不动逻辑。

class AodProcessor: def __init__(self, cfg: dict): self.cfg = cfg self.cache = {} # 缓存边界矢量，避免每个进程重复读 def process_one(self, src): meta = parse_name(src) if meta is None: return None aod = read_sds(src, self.cfg["dataset_name"], scale=self.cfg["scale"], fill=self.cfg["fill_value"]) qa_ok, _ = decode_aod_qa(self._read_qa(src)) valid = np.isfinite(aod) & qa_ok arr = np.where(valid, aod, np.nan) transform = tile_transform(meta["h"], meta["v"]) dst = self._out_path(meta) warp_to_wgs84(src, dst, arr, transform, self.cfg["target_epsg"]) return dst

改造完之后，验证正确性我固定走三步。第一步拿一个月的数据跑新旧两版代码，逐像元做差，最大绝对差应该为 0；第二步抽三景结果和手算的单像元值对照，确保量纲没变；第三步把输出拖进 QGIS 和研究区矢量叠一下，看边界有没有偏移。这三步做完基本能排除 90% 的改造事故。

最容易翻车的地方是"顺手优化"。看到read_sds里有三个循环就想合并成一个，看到float32想改成float64提精度，结果要么改变了无效值的判断逻辑，要么内存翻倍。我的血泪经验是：重构只做结构搬运，不改数值路径，任何涉及计算顺序的改动都单独提一次，用上面三步验一遍。

还有一个实用技巧：给处理类加一个dry_run开关，只解析文件名、算输出路径、检查目标文件是否已存在，不真正读数据。批量任务前面先跑一遍 dry run，能在几分钟内发现路径拼错、目录没权限、瓦片号越界这类问题，比跑到一半再报错省事太多。

从那以后我每次改这套代码，都强制先跑一遍 dry run 加一个月的小样本对照，确认没问题再挂到全量任务上。希望帮到你。

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