简介：面向遥感、生态监测与气候变化方向的学习者和科研人员，这套代码围绕 MCD19A2 地表反射率与 AOD 产品组织批处理流程，重点解决多幅影像在重投影、镶嵌、裁剪、重采样与栅格取值等环节的重复劳动。压缩包共 12 个文件、约 11KB，以 3 个 Python 脚本、3 个 IDL 工程文件（.pro）和 4 个 XML 配置为主，另有工程索引与忽略规则文件；Python 部分对应重投影、TIF 裁剪和 AOD 栅格提取点值，IDL 脚本则借助 MCTK 完成批量镶嵌与均值统计，模块按功能拆分，便于按需调用或串联成完整流水线。已有 1248 人学习下载。对于掌握基础 Python 与 IDL、希望把零散操作整理成可复用模板的读者，可从中参考任务拆分方式、参数组织与批量处理写法，并迁移到其他 MODIS 产品的日常处理中。

1. 从一堆 HDF 文件说起：MCD19A2 代码包到底该先跑哪一步

MCD19A2 数据处理过程的代码，真正难的不是把 HDF 转成 TIFF，而是转完以后敢不敢拿去算均值。前阵子 A 同学拿了一批 MCD19A2 文件，直接批量转栅格，结果 AOD 图层一半是 0，一半是 -28672，均值算出来像玄学。它不是一个双击就开的桌面软件，更像一套能塞进批处理流程的脚本/插件式代码：列子数据集、读 AOD 与 QA、做投影裁剪、再按时间合成。适合手里已经有一堆 HDF、需要把气溶胶光学厚度落到统一网格的人；如果你只想看单景快视图，先别上这套，容易把简单事做重。

2. 把 HDF4 拆开看：MCD19A2 的 AOD、QA 与角度波段怎么读

MCD19A2 的 HDF4 文件不是“一个文件一个波段”的结构，里面挂着一组 EOS_GRID 子数据集，AOD 还分 047 和 055 两个通道，QA 是位掩码，角度波段的量纲又和 AOD 不一样。读错一个环节，后面投影和合成都会跟着翻车。我一般先把“读什么、怎么读、读完是什么量纲”这三件事固定下来，再谈批量。

2.1 HDF4 子数据集不是普通波段，先列清单再读

很多人上来就手拼HDF4_EOS:EOS_GRID:...路径，拼错一个 grid 名就返回 None，或者读出全 0。稳妥做法是用 GDAL 枚举子数据集，按关键字筛出 AOD、QA、角度层，再把完整路径交给gdal.Open。这样换一批文件、换一个 grid 名，脚本不用跟着改。

from osgeo import gdal import re def list_subdatasets(hdf_path): ds = gdal.Open(hdf_path, gdal.GA_ReadOnly) if ds is None: raise FileNotFoundError(hdf_path) subs = ds.GetSubDatasets() out = [] for name, desc in subs: # name 是完整子数据集路径，desc 里带层名 out.append((name, desc)) ds = None return out if __name__ == "__main__": path = "MCD19A2.A2021001.h00v00.000.2021003.hdf" for name, desc in list_subdatasets(path): if re.search(r"Optical_Depth|AOD_QA|Solar_Zenith|View_Zenith", desc, re.I): print(desc) print(name)

逻辑上，GetSubDatasets()返回的是(完整子数据集路径, 描述)，描述里通常带层名，直接拿完整路径去gdal.Open最稳。参数上，hdf_path是单个 HDF 文件路径；筛选关键字不要只写AOD，因为AOD_QA也会命中，最好把Optical_Depth和AOD_QA分开匹配。读子数据集时不要重复打开父文件，父文件用完显式置None，批量跑的时候能少一堆文件句柄问题。

注意：HDF4 驱动没启用时，gdal.Open会直接返回 None。先跑gdalinfo --formats看有没有 HDF4，再排查路径。

2.2 QA 位运算：别把云掩膜当成普通阈值

AOD_QA 不是“大于某个值就保留”的连续量，它是按位打包的分类信息。常见错误是把 QA 转成 float 再做阈值，位运算一参与就全乱。正确做法是先按uint16读入，再用移位和掩码提取目标位段。下面这段代码把“提取位段”和“保留晴空像元”拆开，位偏移和宽度留成参数，方便按手头产品的 QA 表改。

import numpy as np from osgeo import gdal def read_subdataset(hdf_path, key): ds = gdal.Open(hdf_path) subs = ds.GetSubDatasets() hit = [(n, d) for n, d in subs if key.lower() in d.lower()] if not hit: raise KeyError(key) sub = gdal.Open(hit[0][0]) arr = sub.ReadAsArray().astype(np.uint16) sub = None ds = None return arr def extract_bits(qa, offset, width): mask = (1 << width) - 1 return (qa >> offset) & mask qa = read_subdataset("MCD19A2.A2021001.h00v00.000.2021003.hdf", "AOD_QA") cloud = extract_bits(qa, 0, 3) # 低 3 位云掩膜，偏移和宽度按 QA 表改 keep = np.isin(cloud, [0, 1]) # 0 晴空、1 可能晴空，编码以 QA 表为准 print(qa.shape, keep.mean())

逻辑上，extract_bits先构造width位全 1 的掩码，再右移取目标段，最后按位与。参数上，offset是位段起始位，width是位段宽度，keep里的保留值必须对照产品 QA 表，不要照抄示例里的 0 和 1。QA 数组必须用无符号整型读，转成 float 再移位会得到不可预期的结果。批量处理时，我一般先把keep存成布尔掩膜，后面 AOD 合成直接复用，避免每个文件重复算位运算。

2.3 角度波段与尺度因子：AOD 为什么要乘 0.001

MCD19A2 的 AOD 和角度波段大多以整数存储，读出来不乘尺度因子，数值范围会大得离谱。AOD 常见尺度因子是 0.001，角度常见是 0.01 度，填充值常见 -28672，但不同层、不同版本要以元数据为准。下面这段读取函数把“填值转 NaN”和“乘尺度因子”放在同一步，后面做统计就不会把填充值当真实值。

import numpy as np from osgeo import gdal def read_scaled(hdf_path, key, scale, fill=-28672): ds = gdal.Open(hdf_path) sub_path = [n for n, d in ds.GetSubDatasets() if key.lower() in d.lower()][0] sub = gdal.Open(sub_path) arr = sub.ReadAsArray().astype(np.float32) arr[arr == fill] = np.nan arr = arr * scale sub = None ds = None return arr aod = read_scaled("MCD19A2.A2021001.h00v00.000.2021003.hdf", "Optical_Depth_047", 0.001) sza = read_scaled("MCD19A2.A2021001.h00v00.000.2021003.hdf", "Solar_Zenith", 0.01) print(np.nanmin(aod), np.nanmax(aod), np.nanpercentile(aod, [5, 50, 95]))

逻辑上，先读成float32，再把填充值替换成NaN，最后乘尺度因子。参数上，scale按层给，AOD 用 0.001，角度用 0.01；fill不要硬编码成全局常量，最好从子数据集元数据里读，读不到再用常见值兜底。检查时先看 AOD 的 5%、50%、95% 分位数，合理范围通常在 -0.1 到 5 之间，整层超过 5 或大面积 -0.1 以下，先回头查填充值和尺度因子，不要急着做合成。

3. 投影与网格对齐：正弦投影转目标 CRS 的三种落地方式

MCD19A2 原生网格是正弦投影，全球等积，做全球统计很舒服；但一旦要叠矢量、裁区域、和别的栅格做像元级运算，直接拿正弦投影往 GIS 里一扔，错位是常事。投影这一步不是“转一下就行”，目标 CRS、分辨率、范围、重采样方式四个参数要一起定，否则后面合成时网格对不上，前面读得再对也白搭。

3.1 正弦投影不是“错”，但直接给 GIS 用会挨骂：先定目标网格

正弦投影本身没问题，问题在于它和常见 GIS 图层的坐标网格不一致。区域分析我一般先定目标网格四件事：目标 CRS、分辨率、输出范围、重采样方式。做区域统计优先选等面积投影，做叠加展示选经纬度或本地投影；分辨率不要一上来就放大到 30 米，源数据是 1 公里，放大只会让像元看起来更细，信息量不会变多。下面这张表是我批处理前会先填的配置，填完再写代码，比边跑边改省事。

表格里最容易被忽略的是“输出范围外扩半像元”。裁剪时如果范围卡得太死，边缘像元会被切掉一半，统计时有效像元比例偏低，看起来像数据缺失，其实是网格没对齐。另一个是 QA 重采样，很多人图省事把 AOD 和 QA 放同一个gdal.Warp里用 bilinear，结果 QA 出现 1.5、2.7 这种值，后面位运算直接失效。QA 必须单独用最近邻重采样，或者先重投影 AOD，再把 QA 用最近邻对齐到同一网格。

3.2 用 gdal.Warp 做重投影与重采样：参数怎么配

gdal.Warp是批处理里最稳的重投影入口，命令行和 Python 都能用。关键参数不是dstSRS，而是srcNodata和dstNodata：源填充值不声明，双线性会把 -28672 和真实 AOD 一起平均，边缘出现一圈异常低值；目标填充值用 0，后面统计又会把 0 当真实观测。下面这段把常用参数一次性配齐，输出 Float32 GeoTIFF，带压缩和分块。

from osgeo import gdal def warp_aod(src_path, dst_path, dst_srs="EPSG:4326", xres=0.01, yres=0.01): opts = gdal.WarpOptions( format="GTiff", dstSRS=dst_srs, xRes=xres, yRes=yres, resampleAlg="bilinear", srcNodata=-28672, dstNodata=-9999, outputType=gdal.GDT_Float32, multithread=True, warpMemoryLimit=512, creationOptions=["COMPRESS=LZW", "TILED=YES", "BIGTIFF=IF_SAFER"] ) ds = gdal.Warp(dst_path, src_path, options=opts) ds = None warp_aod("aod_047_sub.tif", "aod_047_wgs84.tif")

逻辑上，gdal.Warp先按srcNodata把源填充值排除，再按resampleAlg重采样，最后用dstNodata填输出空洞。参数上，dstSRS填目标 CRS，xRes/yRes填目标分辨率，srcNodata必须和源数据一致，dstNodata用 -9999 这类明显偏离真实范围的值，别用 0。warpMemoryLimit=512是 MB，批量跑大文件时别设太小，否则会频繁落临时文件；multithread=True对多核机器有用，但内存吃紧时先关掉。需要固定输出范围时再加outputBounds，否则默认按源数据范围走。

3.3 按矢量裁剪：掩膜、all_touched 与边界像元

重投影完再裁剪，比先裁剪再重投影更稳，因为裁剪后的边缘像元在重投影时容易被二次重采样。裁剪用rasterio.mask时，all_touched是最容易改变结果的参数：False只保留中心点落在 AOI 内的像元，True会把边界擦到的像元也带进来，后者适合做覆盖统计，不适合做均值统计。下面这段先对齐 CRS，再裁剪并写回 GeoTIFF。

import rasterio from rasterio.mask import mask import geopandas as gpd def clip_by_aoi(src_tif, aoi_shp, dst_tif): with rasterio.open(src_tif) as src: gdf = gpd.read_file(aoi_shp) if gdf.crs != src.crs: gdf = gdf.to_crs(src.crs) geoms = [g.__geo_interface__ for g in gdf.geometry] out_image, out_transform = mask( src, geoms, crop=True, all_touched=False, filled=True, nodata=-9999 ) out_meta = src.meta.copy() out_meta.update({ "driver": "GTiff", "height": out_image.shape[1], "width": out_image.shape[2], "transform": out_transform, "nodata": -9999, "compress": "lzw" }) with rasterio.open(dst_tif, "w", **out_meta) as dst: dst.write(out_image) clip_by_aoi("aod_047_wgs84.tif", "aoi.shp", "aod_047_aoi.tif")

逻辑上，mask按矢量几何裁栅格，crop=True会缩小输出范围，filled=True用nodata填洞。参数上，all_touched=False适合均值统计，True适合覆盖范围检查；nodata=-9999要和前面重投影的输出保持一致。裁剪前必须确认矢量 CRS 和栅格 CRS 一致，不一致先to_crs，不要指望mask自动帮你转。边界像元差异在小区域上很明显，同一块 AOI 用两种all_touched跑出来的均值可能差几个百分点，报告里要写清楚用的是哪一种。

4. 时间序列合成：日尺度 AOD 去噪、填补与月均值落盘

单景 AOD 只能看一天，真正要回答季节变化、污染过程、长期趋势，必须把日文件叠成时间序列。这一步的难点不是np.nanmean，而是日期对齐、QA 权重、最小有效样本数三件事。日期解析错一天，合成结果会整体偏移；QA 权重给得太粗，云边缘会混进来；最小样本数不设，月均值会被一两天的异常值带跑。

4.1 日文件清单与日期解析：别让文件名排序骗了你

文件名里带年积日，字符串排序在单一年份内看起来没问题，跨年、多 tile、缺文件时就会乱。我一般用正则提取年份和年积日，转成datetime再排序，tile 单独分组。下面这段只做日期解析和排序，先把清单跑通，再进合成。

import re from datetime import datetime, timedelta from pathlib import Path def parse_date(path): m = re.search(r"MCD19A2\.A(\d{4})(\d{3})", Path(path).name) if not m: return None year, doy = int(m.group(1)), int(m.group(2)) return datetime(year, 1, 1) + timedelta(days=doy - 1) files = sorted(Path("data").glob("MCD19A2.A*.hdf"), key=parse_date) for f in files[:3]: print(parse_date(f), f.name)

逻辑上，正则抓A后面的 4 位年份和 3 位年积日，再用datetime(year, 1, 1) + timedelta(days=doy-1)转成日期。参数上，doy-1不能少，年积日从 1 开始；key=parse_date让排序按真实日期走，不按文件名字典序。多 tile 场景下，glob出来会混着不同 tile，合成前按文件名里的 tile 字段再分组，否则同一网格会叠进不同区域的数据。

4.2 QA 加权合成：晴天多给权重，云边缘别硬留

日合成不是把所有有效像元一视同仁。晴空像元权重高，可能晴空权重低，云边缘要么给低权重，要么直接剔除。下面这段把 QA 位段转成权重，再做加权平均；权重为 0 的像元不参与分母，避免除零。

import numpy as np def qa_weight(qa, cloud_offset=0, cloud_width=3): cloud = (qa >> cloud_offset) & ((1 << cloud_width) - 1) w = np.zeros_like(qa, dtype=np.float32) w[cloud == 0] = 1.0 w[cloud == 1] = 0.6 w[cloud == 2] = 0.2 return w def weighted_composite(aod_stack, qa_stack): w = np.zeros_like(aod_stack, dtype=np.float32) for i in range(aod_stack.shape[0]): w[i] = qa_weight(qa_stack[i]) w[~np.isfinite(aod_stack)] = 0 num = np.nansum(aod_stack * w, axis=0) den = np.nansum(w, axis=0) out = np.where(den > 0, num / den, np.nan) return out.astype(np.float32)

逻辑上，先按 QA 位段给每个时次生成权重，再把无效 AOD 的权重置 0，最后加权求和除以权重和。参数上，cloud_offset和cloud_width按 QA 表改，权重 1.0/0.6/0.2 是经验值，不是标准答案；云边缘权重给 0.2 是为了保留信息，如果做污染过程分析，可以直接把cloud == 2置 0。np.nansum会忽略 NaN，但aod_stack * w里 NaN 乘 0 仍是 NaN，所以先nansum再判断分母，不要用普通sum。

4.3 缺失值填补与月均值：滑动窗口和最小样本数

月均值不是把一个月所有日文件直接平均，先要设最小有效天数。某个像元一个月只有 2 天有效，平均出来的值代表性很差，不如留空。下面这段按年月分组，先算有效天数，再算均值，低于min_valid的像元置 NaN。

import pandas as pd import numpy as np def monthly_mean(cube, dates, min_valid=5): # cube: (time, y, x) s = pd.Series(range(len(dates)), index=pd.to_datetime(dates)) out = [] for period, idx in s.groupby([s.index.year, s.index.month]): sub = cube[idx.values] valid = np.isfinite(sub).sum(axis=0) mean = np.nanmean(sub, axis=0) mean[valid < min_valid] = np.nan out.append((period, mean)) return out

逻辑上，pd.Series把时间索引和数组位置绑在一起，groupby按年月拿到位置索引，再切cube。参数上，min_valid=5是最低有效天数，区域小、云多的时候可以降到 3，但要在报告里写明；np.nanmean遇到全 NaN 会报警告，可以先算valid再决定要不要算均值。缺失值填补我一般只在展示时用滚动中值，做统计时不填补，填出来的值容易把趋势做平。

5. 避坑与排查：MCD19A2 批处理里最容易翻车的五个环节

这套流程跑通一次不难，难的是批量跑一百个文件时结果还能对得上。下面五条是血泪经验里出现频率最高的，每条按“现象 → 原因 → 解决”写，遇到问题先对照，不要一上来就怀疑算法。

5.1 读取与投影阶段的三个坑

现象：gdal.Open返回 None，或者子数据集读出来全是 0。原因通常是 HDF4 驱动没启用，或者手拼的子数据集路径里 grid 名写错。解决分三步：先跑gdalinfo --formats确认 HDF4 在列表里；再用GetSubDatasets()把完整路径打印出来，不要凭记忆拼；最后拿打印出来的路径单独gdal.Open，能读出数组再进批处理。批量脚本里加一个“父文件打开失败就记录文件名并跳过”的分支，比跑完几百个文件再回头找漏文件省时间。

现象：重投影后 AOD 边缘出现一圈 0 或异常低值。原因基本是srcNodata没设，-28672 被双线性重采样平均进真实像元，边缘就被拉低。解决是把srcNodata和dstNodata都显式写进gdal.WarpOptions，源填充值按元数据填，目标填充值用 -9999，不要用 0。已经转过的旧结果，可以按边缘缓冲区做一次掩膜，但更好的做法是回退到重投影那一步重跑，边缘污染一旦混进时间序列，后面月均值会跟着偏。

现象：QA 图层重采样后出现 1.5、2.7 这种小数。原因是对分类量用了 bilinear 或 cubic。解决是 QA 单独走resampleAlg="nearest"，重投影到和目标 AOD 一样的网格，不要和 AOD 共用一次双线性。如果已经混在一起，把 QA 按最近邻重新对齐，再用位运算提取云掩膜；不要试图对小数 QA 做四舍五入，位段边界会被舍入弄错。

5.2 合成与输出阶段的两个坑

现象：合成结果大面积 0，统计均值被拉低。原因是 NoData 被当成 0 参与计算，读入时没把填充值转 NaN，或者转 NaN 之后又用np.nanmean之前做了普通加法。解决是读入即转 NaN，权重数组里无效像元权重置 0，输出统一用 -9999 表示空值。检查时先打印有效像元比例和 5% 分位数，有效比例低得离谱先查填充值，不要先调合成权重。

现象：批处理跑几十个文件后内存暴涨，或者文件句柄耗尽。原因是循环里gdal.Open后没释放，数组越堆越多。解决是每个文件处理完显式把ds、sub置None，大数组用完del掉；能用rasterio.open的上下文管理器就用上下文管理器。分块处理比一次性读整幅稳，尤其是多时次叠 cube 的时候，先落盘再合并，不要在内存里硬撑。

提示：批处理脚本先拿 3 个文件跑通，再扩到 30 个，最后全量。直接全量跑，出错定位成本会高很多。

6. 验证与进阶：用 QA 分层和交叉检查把结果钉死

结果写完不算完，我一般先跑三张检查表：有效像元比例、AOD 分位数、QA 分层均值。有效像元比例低于 10%，先回头查云掩膜和裁剪范围；AOD 的 5% 到 95% 分位数如果整体贴着 0 或者超过 5，先查尺度因子和填充值；QA 分层均值如果晴空和云边缘差不多，说明权重没起作用。下面这段快速检查代码不依赖额外库，读完栅格就能打印。

import numpy as np import rasterio def quick_check(tif_path): with rasterio.open(tif_path) as src: aod = src.read(1).astype(np.float32) nodata = src.nodata if nodata is not None: aod[aod == nodata] = np.nan valid = np.isfinite(aod) print("valid_ratio", valid.mean()) print("p5_p50_p95", np.nanpercentile(aod, [5, 50, 95])) print("mean", np.nanmean(aod)) quick_check("aod_047_aoi.tif")

逻辑上，先把nodata转 NaN，再算有效比例和分位数；参数上，p5_p50_p95比单一均值更能看出异常，valid_ratio低于阈值时不要继续合成。进阶一点的做法是把 QA 分层和 AOD 均值交叉统计，晴空像元的均值应该明显低于云边缘，如果两者几乎一样，说明 QA 位段取错了或者权重没写进合成。

我现在的习惯是每批结果先跑quick_check，再看 QA 分层表，不通过就不回写月均值。以前图省事跳过这一步，后面返工把投影、裁剪、合成全重跑了一遍，那个后悔药真不好吃。希望帮到你。

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