MCD19A2 AOD数据从HDF4到GeoTIFF批量转换与QA掩码处理

简介：面向遥感数据处理初学者与需要批量处理MCD19A2产品的科研人员，这份代码包围绕MODIS地表反射率与植被指数产品的重投影、融合、镶嵌、重采样、裁剪及栅格点值提取等环节，提供可复用的IDL与Python实现思路，适合已有基础地理信息与IDL编程概念、希望减少重复手工操作的用户。压缩包共12个文件，约11KB，主要由3个pro脚本、3个py脚本以及xml、iml等项目配置和忽略文件组成，分别承担IDL批处理逻辑、Python辅助处理与工程环境描述。内容预览显示其涉及AOD栅格提取点值、均值统计、批量镶嵌和MCTK调用等模块，可帮助读者理解多步骤遥感流程如何拆分与串联，并据此改造自己的处理脚本。目前已有1248人学习下载，适合作为小体量参考代码查阅。

1. MCD19A2 数据处理代码：先把一堆 hdf 变成能进 GIS 的栅格

某个下午，A 同学把一块移动硬盘塞过来，里面整整齐齐几千个MCD19A2.A2019***.h27v05.061.*.hdf，说“数据下好了，帮我打开”。双击打不开，拖进 QGIS 只认出第一层，写脚本读出来全是负数。MCD19A2 是 MAIAC 算法反演的 1km 气溶胶光学厚度产品，日尺度、正弦投影、HDF4 容器，一个文件里叠了十几层 SDS，还配一层位运算编码的 QA。这套代码解决三件事：从 HDF4 里精准取层、按 QA 位掩掉云和雪、转成投影和 NoData 都正确的 GeoTIFF 并批量跑。它不是 QGIS 插件，是几段能塞进批处理脚本的独立工具，适合要把 AOD 接进时序分析或 GIS 流程的人。

2. 拆开 HDF4：层、QA 与投影三件事决定了后面所有代码

2.1 先把 SDS 清单列出来，别凭记忆写变量名

HDF4 是容器格式，一个.hdf里可以平铺多个 SDS（Scientific Data Set），也可以再用 Vgroup 组织成组。MCD19A2 走的是 EOS 网格产品那一套，层是平铺的，每层自带自己的_FillValue、scale_factor、valid_range属性。新手最常见的动作是打开官方文档抄一个变量名，然后KeyError——因为不同 Collection 之间层的名字改过，有的层还带前后缀。

我一般先不看文档，直接让工具把文件里有什么打印出来。GDAL 对 HDF4 的支持是把每个 SDS 暴露成一个 subdataset，用gdalinfo就能看到全部名字和尺寸：

# 列出 HDF4 文件里的所有 subdataset，grep 只看名字，避免被元数据刷屏 gdalinfo MCD19A2.A2019001.h27v05.061.2020001.hdf | grep "SUBDATASET_.*_NAME"

输出形如HDF4_EOS:EOS_GRID:"文件路径":MODIS_Grid_1km_2D:Optical_Depth_055，冒号最后一段才是层名。这条命令的价值在于：文件名、层名、层数一次对齐，后面写代码不用猜。缺点是gdalinfo对几百个文件逐个跑很慢，批量场景我会换成 pyhdf 直接读：

from pyhdf.SD import SD, SDC path = "MCD19A2.A2019001.h27v05.061.2020001.hdf" hdf = SD(path, SDC.READ) # datasets() 返回 {层名: (维度, 数据类型, 维度列表)}，一次拿全 for name, (ndim, dtype, dims) in hdf.datasets().items(): print(f"{name:28s} ndim={ndim} dtype={dtype} dims={dims}") hdf.end()

SD打开的是 SDS 接口，datasets()返回的字典里，值的第一项是维数、第二项是 HDF 数据类型常量、第三项是各维长度。注意hdf.end()必须调，HDF4 的文件句柄不释放，批量跑几千个文件时会把文件描述符耗光，报Too many open files。这段代码跑一次，你就拿到了这个产品的真实层清单，后面所有变量名从这里抄，不靠记忆。

2.2 层与层是配对的：AOD 不能单独用

MCD19A2 里最常被误用的关系是“只读 AOD”。AOD 层本身是 int16，有效值域只在 0 到 5000（缩放后 0 到 5）之间，其余全是填充值和被算法标记过的无效像元。真正决定某个像元能不能用的是AOD_QA，而 AOD 的不确定度在AOD_Uncertainty里。这三层必须同尺寸、同网格、同一批读出来，任何一层缺失，后面统计出来的均值都是错的。

我通常关注的核心层如下表。列这张表的目的不是让你背，而是让你知道哪些层要成组读、哪些层只在特定分析里才需要，避免每个文件读十几层把 IO 拖死。

“按需”这两个字的实际含义是：如果你只做 AOD 的月均或年际变化，读 055 和 QA 两层就够，一个文件从 40MB 级别的读取量降到十几 MB，几千个文件累积下来是几十分钟和几小时的差别。反过来，要做大气校正或者反演 Angstrom 指数，047、Column_WV、角度层就都跑不掉，这时候别省。

提示：每层的scale_factor和_FillValue一定从文件属性里读，不要硬编码。产品升级时这两个值调整过，硬编码的结果是悄无声息地算出一堆 0。

2.3 网格、投影与瓦片编号：为什么同一个 h27v05 能直接叠

MCD19A2 用的是全球正弦投影（Sinusoidal），球体半径 6371007.181 m，1km 网格，每个瓦片 1200 行 × 1200 列。文件名里的h27v05是瓦片编号，水平方向 0 到 35、垂直方向 0 到 17，同一个编号的瓦片在任何一天都落在完全相同的平面坐标范围里——这就是为什么你可以把同一瓦片的一整年数据直接按数组相加，而不需要先做几何配准。

投影参数不在 SDS 里，藏在文件的StructMetadata.0属性中，是一段带换行的结构化文本，里面有UpperLeftPointMtrs=(x,y)、XDim、YDim。把它正则抓出来，配合固定像元大小，就能拼出 GeoTIFF 需要的六参数：

import re from pyhdf.SD import SD, SDC PIXEL = 926.625433 # MODIS 正弦投影 1km 网格像元边长，单位米 def read_upper_left(path): """从 StructMetadata.0 里取瓦片左上角平面坐标""" hdf = SD(path, SDC.READ) meta = hdf.attributes().get("StructMetadata.0", "") hdf.end() m = re.search(r"UpperLeftPointMtrs=\(([-\d.]+),([-\d.]+)\)", meta) if not m: raise ValueError(f"{path} 里没找到 UpperLeftPointMtrs，检查是否被裁剪过") return float(m.group(1)), float(m.group(2))

这里有两个容易被忽略的点。第一，UpperLeftPointMtrs是左上角像元的外边缘坐标，写 geotransform 时直接作为原点，像元大小取负值给 Y 方向，别自己再加半个像元，加了就整体偏移 463 米，肉眼看不出来，和别的图层叠加时才发现对不上。第二，PIXEL这个值不是整数 1000，任何拿 1000 去写 geotransform 的脚本，在跨瓦片拼接时都会累积出几百米的缝。

3. 从 HDF4 到 GeoTIFF：一条能批跑的转换链路

3.1 环境与依赖：GDAL 和 pyhdf 的安装顺序有讲究

这套流程依赖两个东西：GDAL（负责写 GeoTIFF、重投影、拼接）和 pyhdf（负责读 HDF4 的 SDS 和属性）。两者在 conda 环境下共存最省事：

# 建议单独建环境，pyhdf 对 HDF4 库版本敏感，别和主环境混 conda create -n mcd19 python=3.10 -y conda activate mcd19 conda install -c conda-forge gdal pyhdf numpy -y # 验证：两个都要能 import，且 GDAL 要能看到 HDF4 驱动 python -c "from osgeo import gdal; print(gdal.__version__); print(gdal.GetDriverByName('HDF4') is not None)"

GetDriverByName('HDF4')返回None说明 GDAL 编译时没带 HDF4 支持，这时候 pyhdf 读数组没问题，但后面想用 GDAL 直接打开 subdataset 就会失败。遇到这种情况，要么换 conda-forge 的 gdal 包，要么整条链路都走 pyhdf 读、纯 numpy 算、GDAL 只负责写，本文的脚本就是后者，不依赖 GDAL 的 HDF4 驱动。

pyhdf装不上的场景也常见，尤其在已经装了系统级 HDF4 的机器上，编译时会链接到错误版本的库。我的处理办法是优先用 conda 包而不是 pip 源码编译，pip 那条路在无网环境里失败的次数太多。

3.2 单文件转换脚本：读、掩、写三步走完

下面这段是一个完整的单文件转换函数，输入一个 hdf 路径，输出一个只含有效 AOD 的 GeoTIFF。它把 2.2 的配对、2.3 的投影、QA 掩码全串起来了：

import os import numpy as np from osgeo import gdal, osr from pyhdf.SD import SD, SDC PIXEL = 926.625433 AOD_NAME = "Optical_Depth_055" QA_NAME = "AOD_QA" FILL = -28672 # AOD 层填充值，运行时以文件属性为准 SCALE = 0.001 QA_KEEP = (0, 1) # 云掩码 0=晴空 1=可能晴空，其余丢弃 def hdf_to_gtiff(src, dst, keep_qa=QA_KEEP): hdf = SD(src, SDC.READ) # 1) 读 AOD 并转浮点 sds = hdf.select(AOD_NAME) aod = sds.get().astype(np.float32) attrs = sds.attributes() fill = attrs.get("_FillValue", FILL) scale = attrs.get("scale_factor", SCALE) # 2) 读 QA，取低 3 位做云掩码 qa = hdf.select(QA_NAME).get() cloud = qa & 0b111 valid = np.isin(cloud, keep_qa) & (aod != fill) & (aod > 0) # 3) 缩放 + 无效置 NaN out = np.where(valid, aod * scale, np.nan).astype(np.float32) # 4) 投影与地理变换 meta = hdf.attributes().get("StructMetadata.0", "") hdf.end() import re m = re.search(r"UpperLeftPointMtrs=\(([-\d.]+),([-\d.]+)\)", meta) x0, y0 = float(m.group(1)), float(m.group(2)) srs = osr.SpatialReference() srs.ImportFromProj4("+proj=sinu +R=6371007.181 +units=m +no_defs") driver = gdal.GetDriverByName("GTiff") ds = driver.Create(dst, out.shape[1], out.shape[0], 1, gdal.GDT_Float32, options=["COMPRESS=DEFLATE", "PREDICTOR=3", "TILED=YES"]) ds.SetProjection(srs.ExportToWkt()) ds.SetGeoTransform((x0, PIXEL, 0, y0, 0, -PIXEL)) band = ds.GetRasterBand(1) band.SetNoDataValue(np.nan) band.WriteArray(out) band.FlushCache() ds = None return int(valid.sum()), out.size

逻辑说明按四步看。第一步sds.get()拿到的是原始 int16 数组，不缩放，缩放放在掩码之后做，避免把填充值 -28672 乘成 -28.672 再判断，那样填充值和真实的负值就分不开了。第二步是关键：AOD_QA的低 3 位是云掩码，0 和 1 是晴空与可能晴空，2 和 3 是有云与可能有云，& 0b111只保留这三位，剩下的位段不参与判断。第三步用np.where一次性把无效像元写成NaN，而不是写 0——写 0 会让后续求均值时把无效像元算成“极清洁天气”，AOD 均值被严重低估，这是血泪经验里最常见的一条。第四步的PREDICTOR=3是浮点预测器，对 AOD 这种平滑场压缩率提升明显，代价是编码稍慢。

参数说明：keep_qa默认(0, 1)是标准档；FILL和SCALE只是兜底，运行时优先取文件属性；SetNoDataValue(np.nan)在部分老版本 GDAL 里会被写成-3.4e38，读回来判断时要用np.isnan而不是等值比较。

3.3 参数怎么调：数据类型、NoData 与压缩的取舍

这一段专门说参数，因为同样的脚本，参数换一换结果就完全不同。第一是数据类型。AOD 原始是 int16、缩放 0.001，转成 Float32 后每个像元 4 字节，1200×1200 一个文件约 5.7MB，压缩后通常 1 到 2MB。如果你的应用只关心 0.01 的精度，可以用 Int16 存“缩放后乘 10000”的整数，体积减半，但下游每个读取环节都得记得除以 10000，一旦有人忘了，量级直接错。我一般保持 Float32，省下的是沟通成本。

第二是 NoData。用NaN最省心，所有 numpy 的nanmean、nanpercentile直接可用；代价是个别老版本 GIS 软件对 NaN 的 NoData 支持不好，显示成全黑或者全白。要兼容这些软件，就写-9999.0并在SetNoDataValue里声明，代码里判断有效像元时把np.isnan换成out != -9999。

第三是压缩。COMPRESS=DEFLATE加TILED=YES是遥感栅格的常规组合，BLOCKXSIZE=256、BLOCKYSIZE=256让分块读取时序数据时只解压用到的部分。不建议用COMPRESS=LZW处理浮点 AOD，压缩率明显低于 DEFLATE。

3.4 批量跑：文件名解析、并发与断点续跑

单文件跑通之后，批量的问题立刻变成三个：怎么从文件名拿到日期和瓦片、怎么并行、跑一半断了怎么接上。文件名格式是MCD19A2.AYYYYDDD.hXXvYY.CCC.YYYYDDDHHMMSS.hdf，按点号切分后第 2 段是年积日，第 3 段是瓦片：

import os, glob from concurrent.futures import ProcessPoolExecutor from datetime import datetime, timedelta def parse_name(path): parts = os.path.basename(path).split(".") doy = parts[1][1:] # A2019001 -> 2019001 year, day = int(doy[:4]), int(doy[4:]) date = datetime(year, 1, 1) + timedelta(days=day - 1) return date.strftime("%Y%m%d"), parts[2] # ('20190101', 'h27v05') def run_one(src, outdir): date, tile = parse_name(src) dst = os.path.join(outdir, f"AOD_{tile}_{date}.tif") if os.path.exists(dst) and os.path.getsize(dst) > 0: return dst, "skip" # 断点续跑：成品存在就跳过 valid, total = hdf_to_gtiff(src, dst) return dst, f"{valid}/{total}" if __name__ == "__main__": files = sorted(glob.glob("/data/mcd19a2/*.hdf")) with ProcessPoolExecutor(max_workers=8) as ex: for dst, msg in ex.map(run_one, files, ["/data/out"] * len(files)): print(dst, msg)

parse_name里把年积日转成公历日期，是为了后面做时序合成时能直接按天排序，而不是按字符串排。run_one开头的存在性检查就是断点续跑的实现，几千个文件的批量任务一定会中断，没有这行就得从头再来。max_workers=8不是越多越好：pyhdf读文件是 IO 密集，但写 GeoTIFF 时 DEFLATE 压缩吃 CPU，实测 8 到 12 个进程在普通工作站上比较均衡，再往上加会因为内存带宽和磁盘争抢反而变慢。

4. AOD_QA 位掩码：筛得松和筛得严，结论能差一倍

4.1 QA 是位运算，不是数值大小

AOD_QA是 uint8，一个字节里塞了多个字段，每个字段占若干 bit。云掩码在低 3 位，邻域、冰雪、地表类型、气候态 AOD 等字段分布在更高的位上。最常见的误用是把它当普通数值比较，比如qa < 3或者qa == 0——前者恰好能过滤掉一部分云，看起来“有效果”，实际把高位的字段一起卷进来了；后者只保留完全干净的像元，样本量掉得厉害。

正确的读法是先移位再取掩码：

import numpy as np # 位定义以产品 User Guide 为准，这里用常量表隔离，换版本只改这张表 CLOUD_MASK = 0b00000111 # bit 0-2 ADJ_MASK = 0b00001000 # bit 3 邻域像元 SNOW_MASK = 0b00110000 # bit 4-5 冰雪（示意，按实际文档核对） def decode_qa(qa): cloud = qa & CLOUD_MASK adj = (qa & ADJ_MASK) >> 3 snow = (qa & SNOW_MASK) >> 4 return cloud, adj, snow

把位段写成常量表而不是散落在代码里，是因为产品升级时位定义调整过，一旦散落，改起来必然漏。>> 3这种右移是为了把字段挪到最低位，后续做== 1之类的判断才直观。

注意：不同 Collection 的AOD_QA位定义确实改过，拿到新数据先对照随产品附带的 User Guide 核对一遍，别直接套老脚本。

4.2 一个可复用的掩码函数

把判断逻辑收敛到一个函数里，好处是三个档位共用一套代码，只传参：

def aod_mask(aod, qa, level="standard", fill=-28672): """返回布尔掩码，True 表示该像元可用""" cloud, adj, snow = decode_qa(qa) base = (aod != fill) & (aod > 0) if level == "loose": # 只去有云，保留可能晴空和邻域像元 return base & (cloud <= 1) elif level == "standard": # 晴空 + 可能晴空，且非邻域污染 return base & (cloud <= 1) & (adj == 0) elif level == "strict": # 只要最干净的晴空，且无雪 return base & (cloud == 0) & (adj == 0) & (snow == 0) raise ValueError(f"未知档位: {level}")

base先排除填充值和小于等于 0 的像元——AOD 为负在物理上不成立，产品里出现负值是反演残差，参与统计会拉低均值。三个档位的差别集中在云、邻域、雪三个字段上。loose适合做长时序趋势，样本量优先；standard是我做区域年均值时的默认档；strict用于单日个例分析或者和地基观测做点对点比对，宁可样本少也要干净。

4.3 三档策略的样本量与偏差权衡

换档位不是“越严越好”，它直接改变样本量和偏差方向。用同一批数据跑一次统计，把有效像元比例打出来，比任何理论说明都直观：

最后一行要注意：strict不是“更准”，而是“更干净”。被剔除的往往是有云边缘和污染过程，留下的全是洁净晴空，年均值反而偏高。这个反直觉的结论我在第一次做区域对比时吃过亏——换了严格掩码之后 AOD 均值不降反升，查了半天才发现是样本选择偏了，不是代码写错。

5. 排查：转换链路上最容易翻车的五个地方

5.1 读写阶段的问题

现象一：读出来的数组全是 -28672，缩放后变成 -28.672。原因是没有判断填充值就做缩放，或者把_FillValue硬编码成了别的数。解决：先读属性拿_FillValue，在缩放之前用aod == fill建掩码，缩放放最后。顺便检查valid_range属性，有的产品会把无效区间的上下界也写进去，双重校验更稳。

现象二：批量跑到几百个文件时报Too many open files。原因是pyhdf.SD打开后没调end()，句柄一直挂着，多进程下更快耗尽。解决：所有读文件的函数用try/finally保证end()执行，或者像 3.2 的脚本那样在读完后立刻hdf.end()再进写流程。系统层的ulimit -n可以临时调大，但那是掩盖问题不是解决。

现象三：GeoTIFF 在 QGIS 里位置没错，和另一套投影数据叠一起整体偏了几百米。原因是UpperLeftPointMtrs是像元外边缘坐标，脚本里又给 X、Y 各加了半个像元。解决：原点直接用，不偏移；Y 方向像元大小取负。验证办法是把同一个瓦片的两个不同产品叠一起，如果 AOD 的云纹理和地表反射率的云纹理对不齐，就是几何错了，肉眼一眼能看出来。

5.2 质量控制与几何阶段的问题

现象四：月均 AOD 出现大片 0 值，像被抠掉了一块。原因是无效像元被写成了 0 而不是 NaN，做均值时 0 参与了计算。解决：写入前统一置 NaN，统计时用np.nanmean；如果已经写成 0，用np.where(arr == 0, np.nan, arr)补救，但要注意真实的 AOD=0 极少，不会造成明显误伤。

现象五：跨瓦片拼接后中间有一条几十米宽的缝。原因有两个，一是像元大小用了 1000 而不是 926.625433，二是拼接时按文件边界而不是按平面坐标对齐。解决：统一用正弦投影的真实像元尺寸，拼接走 GDAL 的BuildVRT或Warp按坐标对齐，不要用数组切片硬拼。判断方法很简单，把拼接结果转成 WGS84 看经纬度网格，缝会变成一条线性的 NaN 带。

6. 重投影、拼接与时序合成：把日产品变成能用的时间序列

单日 GeoTIFF 只是中间产物，真正的分析对象是重投影后的区域序列。我一般用 GDAL 的 Python 绑定做重投影，不自己写坐标变换：

from osgeo import gdal def reproject(src, dst, epsg=4326, res=0.01): opts = gdal.WarpOptions( dstSRS=f"EPSG:{epsg}", xRes=res, yRes=res, resampleAlg="bilinear", # AOD 是连续场，双线性比最近邻合理 srcNodata=float("nan"), dstNodata=float("nan"), multithread=True, warpMemoryLimit=512, ) gdal.Warp(dst, src, options=opts)

resampleAlg选双线性而不是最近邻，是因为 AOD 是连续物理量，最近邻会在重采样后留下块状伪影；但 QA 层必须用最近邻，类别型数据做插值毫无意义。res=0.01度大约对应 1km，和原始分辨率匹配，设得太细只会引入冗余像元。

区域拼接用 VRT 比直接拼文件省事，虚拟文件不占空间，后面做合成时 GDAL 会按需读取：

# 把一天内覆盖研究区的所有瓦片做成一个虚拟镶嵌 gdalbuildvrt -srcnodata nan -vrtnodata nan merged_20190101.vrt AOD_h27v05_20190101.tif AOD_h28v05_20190101.tif

时序合成的坑在于“按天合成”和“按像元合成”的区别。正确做法是先把一年 365 个重投影后的栅格堆成三维数组（或者用gdal_calc分块算），再沿时间轴做nanmean：

import glob import numpy as np from osgeo import gdal def stack_mean(pattern, out_path): files = sorted(glob.glob(pattern)) arrs = [] for f in files: ds = gdal.Open(f) arrs.append(ds.ReadAsArray().astype(np.float32)) ds = None cube = np.stack(arrs, axis=0) # (T, H, W) mean = np.nanmean(cube, axis=0) n = np.sum(~np.isnan(cube), axis=0) # 每像元有效观测天数 ref = gdal.Open(files[0]) drv = gdal.GetDriverByName("GTiff") out = drv.Create(out_path, ref.RasterXSize, ref.RasterYSize, 2, gdal.GDT_Float32) out.SetGeoTransform(ref.GetGeoTransform()) out.SetProjection(ref.GetProjection()) out.GetRasterBand(1).WriteArray(mean) out.GetRasterBand(2).WriteArray(n.astype(np.float32)) out = None

第二个波段存有效观测天数，作用比均值本身还大：任何一个像元的年均值，如果背后只有 3 天观测，它就没有解释力。我习惯在下游分析前先按天数过滤，少于某个阈值（比如月 5 天、年 30 天）的像元直接置 NaN，这一步能挡掉大部分“看起来有值其实不可信”的格子。验证方法也简单：随机抽几十个像元，把原始 hdf 里对应位置的 AOD 和 QA 打出来，手工推一遍，和成品比对，数值对不上就回头查缩放和掩码。

从那以后我每批数据落地前都强制跑一遍三个自检——数组尺寸是否等于 1200×1200、值域是否落在 0 到 5、有效像元比例是否在合理区间，任何一项不对就整批停下来查，不再往下游传。这套习惯救过我好几次，希望帮到你。

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