简介MCD19A2数据处理代码包面向遥感与地球科学领域的开发者、研究生及科研人员聚焦MODIS气溶胶光学厚度与地表反射率产品的批量处理流程帮助读者跳过重复手工操作、快速搭建可复用的处理链路。压缩包共12个文件、约11KB以py脚本、pro程序和xml配置为主Python脚本负责重投影、栅格裁剪和AOD值提取到点IDL的pro文件实现均值统计、批量镶嵌以及基于MCTK的融合与重采样xml、iml等用于工程配置适合在IDL与Python混合环境中直接运行并二次修改。已有1248人学习下载说明这套处理流程在同类遥感任务中认可度较高可作为入门与进阶的参考模板。代码覆盖融合、镶嵌、重采样、裁剪、栅格取值等核心环节读者可对照自身数据组织方式理解NDVI等指数的计算与空间匹配思路并据此扩展自己的批量处理脚本。1. 拿到 MCD19A2 的代码之后先搞清楚它替你省掉了哪几件事一个瓦片攒够一年是 365 个 HDF4 文件每个文件里躺着 1200×1200 的整型数组直接拖进桌面 GIS 大概率卡在第一步——软件根本列不出里面的子数据集。MCD19A2 数据处理代码要解决的正是这一段把每天一份的 HDF4 拆成带地理参考的 AOD 栅格顺手把缩放因子、填充值、QA 位掩膜和正弦投影一起处理掉。它是 MAIAC 算法产出的 1km 气溶胶光学厚度产品Terra 与 Aqua 双星观测合成一天一景、按瓦片切分。适合做区域气溶胶时序分析、卫星-AERONET 类地基比对、以及要把 AOD 当作协变量喂进 atmospheric correction 流程的人。如果你只是偶尔打开一两个文件看看手动点两下就够一旦上了年尺度、多瓦片代码就是必需品。2. 拆开 HDF4SDS 层、缩放因子与文件自带的正弦投影HDF4 这层壳劝退过不少人。它不像 GeoTIFF 那样“打开就是一张图”一个 MCD19A2 文件里同时塞了几十个 SDSScientific Data Set每个 SDS 有自己的数据类型、缩放因子和填充值而地理参考不在 SDS 上藏在文件级的StructMetadata.0字符串里。想让后面的批处理不出岔子这一章的三件事必须先落地读哪几层、怎么还原物理量、坐标从哪来。2.1 一份文件里真正要读的只有五六层一个 MCD19A2 文件通常带 1200×1200 的网格SDS 数量在二十上下。全读一遍既慢又没必要常用的就这几层SDS 名称含义数据类型典型用途Optical_Depth_0470.47 μm 通道 AODint16与 055 联合做 Angstrom 推算Optical_Depth_055550 nm 通道 AODint16主产品绝大多数分析用它AOD_QA质量位掩膜uint16云、雪、邻云、置信度全在这AOD_Uncertainty反演不确定度int16做加权平均或精度筛选AOD_Model气溶胶模型编号uint8区分人为/沙尘主导Column_WV柱水汽int16辅助分析、做水汽校正前四层基本是标配。Optical_Depth_047和Optical_Depth_055之间不是简单换算关系别拿 0.47 的值乘个系数当 550 nm 用要做波长插值用两层一起拟合 Angstrom 指数更靠谱这也是代码包里单独提供 047 层读取的原因。AOD_QA是整份数据里最容易读错的一层。它是位打包的 uint16不同 bit 段存不同信息值域不是“越大越好”那种连续质量分。很多人图省事直接用AOD_QA 0当筛选条件结果得到一片近乎全空的图——因为 bit 段里还编码了地表类型、雪盖这些跟质量无关的标记。2.2 用 pyhdf 把 SDS 和它的属性一起掏出来读取 HDF4 的 Python 侧方案里pyhdf是最直接的一个不经过 GDAL不依赖驱动编译纯读数组。代价是它只给你数字坐标得自己算。from pyhdf.SD import SD, SDC import numpy as np f SD(MCD19A2.A2019001.h25v04.061.2019002143456.hdf, SDC.READ) # 1) 取原始整型数组不要在这一步做任何数学运算 sds f.select(Optical_Depth_055) raw sds.get() # shape (1200, 1200), dtype int16 # 2) 属性必须跟着数据一起读缩放和填充都在里面 attrs sds.attributes() scale attrs.get(scale_factor, 1.0) # 常见为 0.001 fill attrs.get(_FillValue, -28672) # 不要硬编码不同版本会变 valid_min attrs.get(valid_range, [-100, 5000])[0] valid_max attrs.get(valid_range, [-100, 5000])[1] # 3) 先掩膜再缩放顺序反了会把填充值也乘进去 aod raw.astype(float32) aod[raw fill] np.nan aod aod * scale sds.endaccess() f.end()这段的核心是第 3 步的顺序。raw fill这个判断必须在乘scale之前做因为填充值是整型域里的哨兵值比如 -28672乘完 0.001 之后变成 -28.672你再去和_FillValue比对就对不上了。第二处是valid_range它不是用来掩膜的是用来做物理范围检查的AOD 理论上可以取到负值反演噪声导致所以下限一般是 -0.1 这个量级而不是 0。把下限直接卡成 0 会系统性抬高区域均值这个偏差在污染背景下的站点比对里能看出来。scale_factor用attrs.get带默认值是防止某些版本没写这个属性。真遇到没写的说明该层本身就是物理量直接跳过缩放即可不要无条件乘 0.001。2.3 地理参考藏在 StructMetadata 里别按标准瓦片硬算这是 MCD19A2 最容易翻车的地方。MAIAC 的 1km 网格虽然也是正弦投影但它的原点定义和 MOD09、MOD04 那套标准正弦瓦片不是同一个直接套“h25v04 左上角坐标 某公式”去算行列号结果会整体偏移一段距离。偏移量不大但足够让你的采样点落错像元。正确做法是读文件自带的网格定义import re meta f.attributes().get(StructMetadata.0, ) # 网格尺寸与四个角点单位是投影米 ul re.search(rUpperLeftPointMtrs\(([-\d\.]),([-\d\.])\), meta) lr re.search(rLowerRightPointMtrs\(([-\d\.]),([-\d\.])\), meta) dim re.search(rXDim(\d)\s*\n\s*YDim(\d), meta) x0, y0 float(ul.group(1)), float(ul.group(2)) # 左上角 x1, y1 float(lr.group(1)), float(lr.group(2)) # 右下角 nx, ny int(dim.group(1)), int(dim.group(2)) res_x (x1 - x0) / nx # 约 1000 米但不一定正好是 1000 res_y (y0 - y1) / ny def rowcol_to_xy(row, col): 行列号 - 正弦投影坐标米 x x0 (col 0.5) * res_x y y0 - (row 0.5) * res_y return x, y0.5是取像元中心这个细节看着小但在 1km 分辨率下做点位匹配时半个像元就是 500 米。ProjParams里还能拿到球体半径正弦投影用的是球体近似半径约 6371007.181 米写 proj4 字符串时要把这个值填进去用 WGS84 椭球的参数去套会带来米级到十米级的系统性误差。提示StructMetadata.0是一段超长字符串不同版本换行符和空格数不一样。正则里用\s*容错别写死成单空格。到这里一个文件已经被拆成“数组 属性 坐标”三件套了。剩下的是怎么把这一步放大到几千个文件以及怎么把正弦投影的坐标变成别人能用的经纬度。3. 批量读取与重投影GDAL 子数据集和 pyhdf 两条路怎么选单个文件怎么读前面已经解决了。真正拉开效率差距的是批处理这一层路径怎么拼、投影怎么转、中断了怎么续。这一段没有哪个现成插件能一把梭哈HDF4 的子数据集命名规则、瓦片的日期编码、重投影时的无效值处理每一环都得自己接上。3.1 GDAL 子数据集路径的拼法与 gdalwarp 参数GDAL 走的是另一条路不把 HDF4 当“文件”看而是当“容器”用HDF4_EOS:EOS_GRID:前缀直接指到某个网格里的某个 SDS。这条路的好处是投影、地理变换、无效值一次到位输出的就是标准 GeoTIFF。第一步永远是先列出子数据集名字别凭记忆拼gdalinfo MCD19A2.A2019001.h25v04.061.2019002143456.hdf | grep SUBDATASET输出里会给出完整字符串形如HDF4_EOS:EOS_GRID:文件路径:MCD19A2:MODIS_Grid_1km_2D:Optical_Depth_055。注意不同版本的前缀层数可能不同有的是产品名:网格名:SDS名有的直接网格名:SDS名所以复制粘贴比手写靠谱。拿到全名之后重投影gdalwarp -of GTiff -overwrite \ -t_srs EPSG:4326 \ -tr 0.01 0.01 \ -r near \ -srcnodata -28672 -dstnodata -9999 \ -te 100 20 120 40 \ HDF4_EOS:EOS_GRID:\MCD19A2.A2019001.h25v04.061.2019002143456.hdf\:MCD19A2:MODIS_Grid_1km_2D:Optical_Depth_055 \ out_2019001.tif参数逐个说清楚-tr 0.01 0.01是输出像元大小0.01 度约 1.1km和源分辨率大致匹配想保细节可以缩到 0.005代价是文件体积四倍。-r near是最近邻用在 AOD 上是有意为之——双线性会把无效区云、填充的值插进有效区边缘一圈全是假数据。-srcnodata必须和文件里的_FillValue对上写错了等于没设。-te是输出范围限定到自己需要的区域能省掉大量空白计算尤其做全国或多瓦片拼接时。如果不想每次手写GDAL 的 Python 绑定更灵活from osgeo import gdal gdal.UseExceptions() sub (HDF4_EOS:EOS_GRID:MCD19A2.A2019001.h25v04.061.hdf: MCD19A2:MODIS_Grid_1km_2D:Optical_Depth_055) src gdal.Open(sub, gdal.GA_ReadOnly) band src.GetRasterBand(1) arr band.ReadAsArray() # int16未缩放 nodata band.GetNoDataValue() # 直接从元数据取不猜gdal.UseExceptions()这一行建议加上。默认情况下 GDAL 出错只返回None批量跑的时候一个静默失败会让人排查到怀疑人生打开异常之后路径错了、驱动缺失了当场抛出来。3.2 pyhdf 读数组 自己算坐标的写法GDAL 省事但有前提你装的 GDAL 必须编译了 HDF4 驱动。自编译或某些 pip 轮子是不带的gdal.Open会直接告诉你“不是支持的文件格式”。这种情况下pyhdf是保底方案——它只要能连上 HDF4 的库就行和 GDAL 无关。pyhdf 路线的代价是重投影要自己做。思路是从StructMetadata.0拿到角点和分辨率用 2.3 里的函数生成一张经纬度查找表再用scipy.ndimage.map_coordinates或rasterio的reproject做重采样。import numpy as np from rasterio.transform import from_bounds from rasterio.warp import reproject, Resampling # 源正弦投影用文件里解析出的 x0/y0/res src_transform from_bounds(x0, y1, x1, y0, nx, ny) src_crs projsinu lon_00 x_00 y_00 R6371007.181 unitsm no_defs # 目标经纬度 0.01 度 dst np.full((2000, 2000), np.nan, dtypefloat32) dst_transform from_bounds(100, 20, 120, 40, 2000, 2000) reproject( sourceaod, destinationdst, src_transformsrc_transform, src_crssrc_crs, dst_transformdst_transform, dst_crsEPSG:4326, resamplingResampling.nearest, src_nodatanp.nan, dst_nodatanp.nan, )关键点是src_transform的构造。from_bounds(left, bottom, right, top, width, height)的参数顺序是左下、右上而文件里给的是左上和右下所以传参时要把 y 换过来写成from_bounds(x0, y1, x1, y0, nx, ny)。这个顺序搞反的表现是图像上下颠倒而数值上看不出来很容易被当成“数据本身有问题”。Resampling.nearest同样是刻意的。AOD 是连续量理论上双线性更平滑但只要源数据里有nan双线性会把nan扩散到邻近像元边缘出现一圈空洞。想在重投影后再平滑那是另一个步骤的事。3.3 批处理脚本目录遍历、日期解析、断点续跑一个瓦片一年 365 个文件多瓦片就是几千个。脚本要解决三件事从文件名里正确解析日期、失败不中断、重跑不重复。文件名格式是MCD19A2.AYYYYDDD.hXXvYY.CCC.时间戳.hdf其中AYYYYDDD是年 年积日DOY不是年月日。这里错过一次就够记一辈子。from pathlib import Path import re, datetime as dt PAT re.compile(rMCD19A2\.A(\d{4})(\d{3})\.(h\d{2}v\d{2})\.(\d{3})\.) def parse_name(p: Path): m PAT.search(p.name) if not m: return None year, doy, tile, ver m.groups() day dt.date(int(year), 1, 1) dt.timedelta(daysint(doy) - 1) return {date: day, tile: tile, version: ver, path: p} def run_batch(root: Path, outdir: Path): outdir.mkdir(parentsTrue, exist_okTrue) for hdf in sorted(root.glob(MCD19A2.*.hdf)): info parse_name(hdf) if info is None: continue out outdir / f{info[tile]}_{info[date]:%Y%m%d}.tif if out.exists() and out.stat().st_size 4096: continue # 断点续跑 try: process_one(hdf, out) # 读取 掩膜 重投影 except Exception as e: (outdir / failed.log).open(a).write(f{hdf.name}\t{e}\n)dt.date(year, 1, 1) timedelta(daysdoy-1)这个写法比strptime更稳因为 Python 的%j在某些版本上和年份组合会有坑。out.stat().st_size 4096这个判断是为了防半截文件——程序被 kill 掉时 GeoTIFF 可能已经创建但没写完只判exists()会把坏文件当成已完成后面统计全是错的。失败写日志而不是直接抛出是为了让一个坏文件不至于毁掉整晚的任务。CC那一段是集合版本号6.1 版本写作061。不同集合的scale_factor和 QA 位定义可能有调整混着处理之前先确认版本一致别把 006 和 061 的文件放进同一个目录。4. QA 位掩膜与缺测把 AOD 从“能读”做到“能用”到这一步AOD 已经是带经纬度的浮点栅格了但直接拿去算均值结论一定是错的。原因有三层云和雪没去掉、填充值没掩干净、极端负值没剔除。这一章就是把这三层剥完让剩下的像元真的能代表大气状态。4.1 用位运算把云掩膜和置信度拆出来AOD_QA是 uint16按 bit 段打包。最常被用到的一组是低 8 位import numpy as np qa f.select(AOD_QA).get().astype(uint16) cloud qa 0b0000000000000111 # bit 0-2云状态 adj (qa 3) 0b1 # bit 3是否紧邻云 surf (qa 4) 0b11 # bit 4-5地表类型 snow (qa 6) 0b11 # bit 6-7雪/冰 clear (cloud 0) (adj 0) (snow 0)位运算的关键是掩码和移位要配对。取 bit 0-2掩码就是0b111取完不用移位取 bit 4-5先右移 4 位再和0b11与顺序反了得到的是原位置的位值完全不对。这类错误不会报错只会让筛选结果看起来“有点奇怪”。具体哪个值代表“晴空”要以文件配套的 ATBD 为准别照抄别人的代码。经验上cloud 0是最严格的一档会用掉相当比例的像元如果能接受稍宽松的条件cloud 1能多留一些样本代价是混入薄云。做趋势分析时筛选条件必须全程一致中途换阈值会让时间序列出现台阶。注意不要用AOD_QA 0做筛选。这个条件实际含义是“所有位段全为 0”等于同时要求地表类型是第一种、无雪、无邻云过于苛刻多数区域会筛成空图。4.2 填充值、物理范围、不确定度三重过滤QA 过了还要过数值关。三层过滤按顺序叠加# 第一层填充值与 NaN raw f.select(Optical_Depth_055).get() fill f.select(Optical_Depth_055).attributes()[_FillValue] aod raw.astype(float32) aod[raw fill] np.nan aod * 0.001 # 第二层物理范围负值不能一刀切 aod[(aod -0.05) | (aod 5.0)] np.nan # 第三层不确定度超过阈值视为不可信 unc f.select(AOD_Uncertainty).get().astype(float32) * 0.001 unc[unc 0] np.nan aod[unc 0.1] np.nan good clear np.isfinite(aod)第二层的下限-0.05是留了余量的。晴空条件下 AOD 真实值可以接近 0反演噪声会让它略微为负把负值全砍掉会引入正偏差——尤其在清洁背景区这个偏差能占到均值的百分之几。上限 5.0 是极端值兜底超过这个数的多半是云边缘残留或地表反射率估计失败。第三层的AOD_Uncertainty是常被忽略的一层。它的量纲同样是 0.001超过 0.1 意味着反演相对误差已经很大。做定量分析时加上这一层样本量会掉一截但和地基观测的相关性能明显改善。4.3 有效像元统计与阈值设定掩膜做完下一个问题是“这一天的这张图到底能不能用”。做法是统计有效像元占比valid_mask good np.isfinite(aod) valid_ratio valid_mask.sum() / valid_mask.size if valid_ratio 0.05: quality_flag insufficient # 有效率太低建议整景丢弃 elif valid_ratio 0.2: quality_flag marginal else: quality_flag ok mean_aod np.nanmean(np.where(valid_mask, aod, np.nan))阈值的取法没有标准答案但有个原则单景阈值和时间合成阈值要分开设。单景有效率低于 5% 基本没法用直接丢做月合成时要求每个月至少有 8 到 10 天有效才输出否则那个格点上的“月均值”其实只代表两三天的天气不具代表性。另一个容易忽略的点是区域统计时的加权。在正弦投影下每个像元的实际面积相等直接求算术平均是对的一旦重投影到等经纬度高纬像元被拉伸不加cos(lat)权重会把高纬度的值算重。这个坑在做跨纬度区域比如整省、整流域统计时特别明显低纬地区还好过了 40 度就必须加权。5. 避坑与排查五条让批处理静默失败的细节下面这五条都是跑批量任务时真会遇到的共同点是不报错或者报错信息指不到根因。按现象、原因、解决三段来记比看文档快。5.1 环境与驱动装得上不等于跑得动现象一gdal.Open返回None或者报not recognized as a supported file format但文件明明存在、路径也没写错用别的工具能打开。原因GDAL 没编译 HDF4 驱动。这不是文件问题是安装包问题。很多通过 pip 安装的 GDAL 轮子为了减小体积只带了 GeoTIFF、JPEG 这类常用驱动HDF4 这种依赖外部库的格式被裁掉了。判断方法很简单gdalinfo --formats | grep -i hdf4没输出就是没带。解决换成带 HDF4 的发行版或者在 conda 环境里装conda install -c conda-forge gdal一般会连同 HDF4 一起解决。如果必须用 pip就得自己编译编译前确认系统里装了 HDF4 的开发库和hdf4-config。这个排查顺序很重要——先查驱动再查代码能省掉半天。现象二import pyhdf报ImportError: libdfalt.so: cannot open shared object file。原因pyhdf 是 HDF4 的 Python 封装运行时要去加载 HDF4 的动态库。pip 装的时候只装 Python 包不会帮你把底层库拉下来所以 import 阶段就断了。解决同样走 conda-forge 安装它会把hdf4作为依赖一起装上。装完用python -c from pyhdf.SD import SD; print(ok)验一下别等到跑批处理才发现。另外如果系统里同时存在多个 HDF4 版本可能出现版本冲突显式指定LD_LIBRARY_PATH指向 conda 环境的 lib 目录能解决大部分情况。5.2 数据与数值结果不对但程序不报错现象三批量跑完区域平均 AOD 比预期低一大截而且低得很有规律像是被什么拉下去了。原因填充值参与了平均。要么没掩_FillValue要么在gdalwarp时没设-srcnodata导致 -28672 这种哨兵值被当成真实数值参与运算。还有一种变体设了-srcnodata但值写错了比如写成了 65535等于没设。解决填充值一律从文件属性里读不硬编码。用 GDAL 时读band.GetNoDataValue()用 pyhdf 时读_FillValue。写死数值的代码在版本升级后一定会出问题因为不同集合的填充值定义可能调整。现象四处理完的 AOD 栅格和 MOD09 的地表反射率叠在一起整体错开一段距离几百米到一两公里。原因拿标准正弦瓦片的行列号公式去套 MCD19A2。MAIAC 用自己的一套 1km 网格定义原点、分辨率都在文件的StructMetadata.0里和 MODIS 标准瓦片网格不是同一个。按标准公式硬算得到的是“看起来对但差一点”的结果最难排查。解决坐标一律从文件自带的元数据里取或者干脆交给 GDAL——它的 HDF4 驱动会读StructMetadata并生成正确的地理变换。自己算坐标只作为备用路径且必须用 2.3 里那套解析逻辑。验证方法挑一个已知地面站位置看它落在哪个像元再和目视解译的湖岸线、海岸线比对。现象五重投影后图像四周出现一圈宽窄不一的无效值带或者整幅图有锯齿状空洞。原因-r bilinear在源数据含nan时会把它扩散到邻域边缘一圈全是插值出来的无效值另一个原因是-tr设得比源分辨率细太多上采样过程中有效像元和无效像元混在一起。解决AOD 这类含大量缺测的连续量重投影阶段用-r near需要平滑的话放到重投影之后单独做。输出分辨率不要盲目缩小0.01 度对应 1.1km 左右和源分辨率匹配就够了要更细的分辨率得接受插值带来的信息损失并且用-dstnodata明确标记。6. 月度合成与交叉校验一个能跑通的最小验证闭环单日图跑通了不代表整条链路是对的。我会用一个月的数据做一次端到端验证合成、出图、和地基观测对一遍三步都过才算这条管线能用。import numpy as np, rasterio, glob files sorted(glob.glob(out/h25v04_2019*.tif)) stack [] for fp in files: with rasterio.open(fp) as ds: stack.append(ds.read(1, maskedTrue).filled(np.nan)) arr np.stack(stack) # (nday, h, w) cnt np.isfinite(arr).sum(axis0) # 每格有效天数 with np.errstate(invalidignore): mean np.nanmean(arr, axis0) mean np.where(cnt 8, mean, np.nan) # 有效天数不足的不输出 with rasterio.open(files[0]) as ds: meta ds.meta.copy() meta.update(dtypefloat32, count1, nodatanp.nan) with rasterio.open(monthly_mean.tif, w, **meta) as dst: dst.write(mean.astype(float32), 1)cnt 8这个阈值不是拍脑袋定的一个月 30 天要求至少 8 天有效意味着允许云覆盖掉七成同时保证均值不至于被两三天代表。阈值调低会让覆盖变好、代表性变差做趋势分析时宁可留空洞也不要虚假填充。合成完要校验。最省事的办法是拿同区域的太阳光度计观测做散点看斜率和相关系数没有地基数据时退一步用同一天的 MOD04 产品做交叉比对主要看空间分布形态是否一致别太在意绝对值的几个百分点差异——两者反演假设不同。几个实操上会卡住的点掩膜后的rasterio读取一定要带maskedTrue否则nodata会被读成 0而 0 在 AOD 里是合法值混进去会把均值拉低np.nanmean在整列全nan时会返回nan并抛 RuntimeWarning用errstate压掉但别用np.mean代替输出时nodata设成nan后续读取能直接识别。血泪经验是重投影那一步的-srcnodata和-dstnodata两个参数——我见过不止一次因为只设了前者、后者默认成 0结果一整月的合成图里出现大块“AOD 恰好等于 0”的区域看着像清洁大气其实是空洞被填了零。从那以后我每次跑完批处理都会先统计一遍零值像元占比超过 5% 就停下来查参数不往下走。希望帮到你。本文还有配套的精品资源点击获取
简介MCD19A2数据处理代码包面向遥感与地球科学领域的开发者、研究生及科研人员聚焦MODIS气溶胶光学厚度与地表反射率产品的批量处理流程帮助读者跳过重复手工操作、快速搭建可复用的处理链路。压缩包共12个文件、约11KB以py脚本、pro程序和xml配置为主Python脚本负责重投影、栅格裁剪和AOD值提取到点IDL的pro文件实现均值统计、批量镶嵌以及基于MCTK的融合与重采样xml、iml等用于工程配置适合在IDL与Python混合环境中直接运行并二次修改。已有1248人学习下载说明这套处理流程在同类遥感任务中认可度较高可作为入门与进阶的参考模板。代码覆盖融合、镶嵌、重采样、裁剪、栅格取值等核心环节读者可对照自身数据组织方式理解NDVI等指数的计算与空间匹配思路并据此扩展自己的批量处理脚本。1. 拿到 MCD19A2 的代码之后先搞清楚它替你省掉了哪几件事一个瓦片攒够一年是 365 个 HDF4 文件每个文件里躺着 1200×1200 的整型数组直接拖进桌面 GIS 大概率卡在第一步——软件根本列不出里面的子数据集。MCD19A2 数据处理代码要解决的正是这一段把每天一份的 HDF4 拆成带地理参考的 AOD 栅格顺手把缩放因子、填充值、QA 位掩膜和正弦投影一起处理掉。它是 MAIAC 算法产出的 1km 气溶胶光学厚度产品Terra 与 Aqua 双星观测合成一天一景、按瓦片切分。适合做区域气溶胶时序分析、卫星-AERONET 类地基比对、以及要把 AOD 当作协变量喂进 atmospheric correction 流程的人。如果你只是偶尔打开一两个文件看看手动点两下就够一旦上了年尺度、多瓦片代码就是必需品。2. 拆开 HDF4SDS 层、缩放因子与文件自带的正弦投影HDF4 这层壳劝退过不少人。它不像 GeoTIFF 那样“打开就是一张图”一个 MCD19A2 文件里同时塞了几十个 SDSScientific Data Set每个 SDS 有自己的数据类型、缩放因子和填充值而地理参考不在 SDS 上藏在文件级的StructMetadata.0字符串里。想让后面的批处理不出岔子这一章的三件事必须先落地读哪几层、怎么还原物理量、坐标从哪来。2.1 一份文件里真正要读的只有五六层一个 MCD19A2 文件通常带 1200×1200 的网格SDS 数量在二十上下。全读一遍既慢又没必要常用的就这几层SDS 名称含义数据类型典型用途Optical_Depth_0470.47 μm 通道 AODint16与 055 联合做 Angstrom 推算Optical_Depth_055550 nm 通道 AODint16主产品绝大多数分析用它AOD_QA质量位掩膜uint16云、雪、邻云、置信度全在这AOD_Uncertainty反演不确定度int16做加权平均或精度筛选AOD_Model气溶胶模型编号uint8区分人为/沙尘主导Column_WV柱水汽int16辅助分析、做水汽校正前四层基本是标配。Optical_Depth_047和Optical_Depth_055之间不是简单换算关系别拿 0.47 的值乘个系数当 550 nm 用要做波长插值用两层一起拟合 Angstrom 指数更靠谱这也是代码包里单独提供 047 层读取的原因。AOD_QA是整份数据里最容易读错的一层。它是位打包的 uint16不同 bit 段存不同信息值域不是“越大越好”那种连续质量分。很多人图省事直接用AOD_QA 0当筛选条件结果得到一片近乎全空的图——因为 bit 段里还编码了地表类型、雪盖这些跟质量无关的标记。2.2 用 pyhdf 把 SDS 和它的属性一起掏出来读取 HDF4 的 Python 侧方案里pyhdf是最直接的一个不经过 GDAL不依赖驱动编译纯读数组。代价是它只给你数字坐标得自己算。from pyhdf.SD import SD, SDC import numpy as np f SD(MCD19A2.A2019001.h25v04.061.2019002143456.hdf, SDC.READ) # 1) 取原始整型数组不要在这一步做任何数学运算 sds f.select(Optical_Depth_055) raw sds.get() # shape (1200, 1200), dtype int16 # 2) 属性必须跟着数据一起读缩放和填充都在里面 attrs sds.attributes() scale attrs.get(scale_factor, 1.0) # 常见为 0.001 fill attrs.get(_FillValue, -28672) # 不要硬编码不同版本会变 valid_min attrs.get(valid_range, [-100, 5000])[0] valid_max attrs.get(valid_range, [-100, 5000])[1] # 3) 先掩膜再缩放顺序反了会把填充值也乘进去 aod raw.astype(float32) aod[raw fill] np.nan aod aod * scale sds.endaccess() f.end()这段的核心是第 3 步的顺序。raw fill这个判断必须在乘scale之前做因为填充值是整型域里的哨兵值比如 -28672乘完 0.001 之后变成 -28.672你再去和_FillValue比对就对不上了。第二处是valid_range它不是用来掩膜的是用来做物理范围检查的AOD 理论上可以取到负值反演噪声导致所以下限一般是 -0.1 这个量级而不是 0。把下限直接卡成 0 会系统性抬高区域均值这个偏差在污染背景下的站点比对里能看出来。scale_factor用attrs.get带默认值是防止某些版本没写这个属性。真遇到没写的说明该层本身就是物理量直接跳过缩放即可不要无条件乘 0.001。2.3 地理参考藏在 StructMetadata 里别按标准瓦片硬算这是 MCD19A2 最容易翻车的地方。MAIAC 的 1km 网格虽然也是正弦投影但它的原点定义和 MOD09、MOD04 那套标准正弦瓦片不是同一个直接套“h25v04 左上角坐标 某公式”去算行列号结果会整体偏移一段距离。偏移量不大但足够让你的采样点落错像元。正确做法是读文件自带的网格定义import re meta f.attributes().get(StructMetadata.0, ) # 网格尺寸与四个角点单位是投影米 ul re.search(rUpperLeftPointMtrs\(([-\d\.]),([-\d\.])\), meta) lr re.search(rLowerRightPointMtrs\(([-\d\.]),([-\d\.])\), meta) dim re.search(rXDim(\d)\s*\n\s*YDim(\d), meta) x0, y0 float(ul.group(1)), float(ul.group(2)) # 左上角 x1, y1 float(lr.group(1)), float(lr.group(2)) # 右下角 nx, ny int(dim.group(1)), int(dim.group(2)) res_x (x1 - x0) / nx # 约 1000 米但不一定正好是 1000 res_y (y0 - y1) / ny def rowcol_to_xy(row, col): 行列号 - 正弦投影坐标米 x x0 (col 0.5) * res_x y y0 - (row 0.5) * res_y return x, y0.5是取像元中心这个细节看着小但在 1km 分辨率下做点位匹配时半个像元就是 500 米。ProjParams里还能拿到球体半径正弦投影用的是球体近似半径约 6371007.181 米写 proj4 字符串时要把这个值填进去用 WGS84 椭球的参数去套会带来米级到十米级的系统性误差。提示StructMetadata.0是一段超长字符串不同版本换行符和空格数不一样。正则里用\s*容错别写死成单空格。到这里一个文件已经被拆成“数组 属性 坐标”三件套了。剩下的是怎么把这一步放大到几千个文件以及怎么把正弦投影的坐标变成别人能用的经纬度。3. 批量读取与重投影GDAL 子数据集和 pyhdf 两条路怎么选单个文件怎么读前面已经解决了。真正拉开效率差距的是批处理这一层路径怎么拼、投影怎么转、中断了怎么续。这一段没有哪个现成插件能一把梭哈HDF4 的子数据集命名规则、瓦片的日期编码、重投影时的无效值处理每一环都得自己接上。3.1 GDAL 子数据集路径的拼法与 gdalwarp 参数GDAL 走的是另一条路不把 HDF4 当“文件”看而是当“容器”用HDF4_EOS:EOS_GRID:前缀直接指到某个网格里的某个 SDS。这条路的好处是投影、地理变换、无效值一次到位输出的就是标准 GeoTIFF。第一步永远是先列出子数据集名字别凭记忆拼gdalinfo MCD19A2.A2019001.h25v04.061.2019002143456.hdf | grep SUBDATASET输出里会给出完整字符串形如HDF4_EOS:EOS_GRID:文件路径:MCD19A2:MODIS_Grid_1km_2D:Optical_Depth_055。注意不同版本的前缀层数可能不同有的是产品名:网格名:SDS名有的直接网格名:SDS名所以复制粘贴比手写靠谱。拿到全名之后重投影gdalwarp -of GTiff -overwrite \ -t_srs EPSG:4326 \ -tr 0.01 0.01 \ -r near \ -srcnodata -28672 -dstnodata -9999 \ -te 100 20 120 40 \ HDF4_EOS:EOS_GRID:\MCD19A2.A2019001.h25v04.061.2019002143456.hdf\:MCD19A2:MODIS_Grid_1km_2D:Optical_Depth_055 \ out_2019001.tif参数逐个说清楚-tr 0.01 0.01是输出像元大小0.01 度约 1.1km和源分辨率大致匹配想保细节可以缩到 0.005代价是文件体积四倍。-r near是最近邻用在 AOD 上是有意为之——双线性会把无效区云、填充的值插进有效区边缘一圈全是假数据。-srcnodata必须和文件里的_FillValue对上写错了等于没设。-te是输出范围限定到自己需要的区域能省掉大量空白计算尤其做全国或多瓦片拼接时。如果不想每次手写GDAL 的 Python 绑定更灵活from osgeo import gdal gdal.UseExceptions() sub (HDF4_EOS:EOS_GRID:MCD19A2.A2019001.h25v04.061.hdf: MCD19A2:MODIS_Grid_1km_2D:Optical_Depth_055) src gdal.Open(sub, gdal.GA_ReadOnly) band src.GetRasterBand(1) arr band.ReadAsArray() # int16未缩放 nodata band.GetNoDataValue() # 直接从元数据取不猜gdal.UseExceptions()这一行建议加上。默认情况下 GDAL 出错只返回None批量跑的时候一个静默失败会让人排查到怀疑人生打开异常之后路径错了、驱动缺失了当场抛出来。3.2 pyhdf 读数组 自己算坐标的写法GDAL 省事但有前提你装的 GDAL 必须编译了 HDF4 驱动。自编译或某些 pip 轮子是不带的gdal.Open会直接告诉你“不是支持的文件格式”。这种情况下pyhdf是保底方案——它只要能连上 HDF4 的库就行和 GDAL 无关。pyhdf 路线的代价是重投影要自己做。思路是从StructMetadata.0拿到角点和分辨率用 2.3 里的函数生成一张经纬度查找表再用scipy.ndimage.map_coordinates或rasterio的reproject做重采样。import numpy as np from rasterio.transform import from_bounds from rasterio.warp import reproject, Resampling # 源正弦投影用文件里解析出的 x0/y0/res src_transform from_bounds(x0, y1, x1, y0, nx, ny) src_crs projsinu lon_00 x_00 y_00 R6371007.181 unitsm no_defs # 目标经纬度 0.01 度 dst np.full((2000, 2000), np.nan, dtypefloat32) dst_transform from_bounds(100, 20, 120, 40, 2000, 2000) reproject( sourceaod, destinationdst, src_transformsrc_transform, src_crssrc_crs, dst_transformdst_transform, dst_crsEPSG:4326, resamplingResampling.nearest, src_nodatanp.nan, dst_nodatanp.nan, )关键点是src_transform的构造。from_bounds(left, bottom, right, top, width, height)的参数顺序是左下、右上而文件里给的是左上和右下所以传参时要把 y 换过来写成from_bounds(x0, y1, x1, y0, nx, ny)。这个顺序搞反的表现是图像上下颠倒而数值上看不出来很容易被当成“数据本身有问题”。Resampling.nearest同样是刻意的。AOD 是连续量理论上双线性更平滑但只要源数据里有nan双线性会把nan扩散到邻近像元边缘出现一圈空洞。想在重投影后再平滑那是另一个步骤的事。3.3 批处理脚本目录遍历、日期解析、断点续跑一个瓦片一年 365 个文件多瓦片就是几千个。脚本要解决三件事从文件名里正确解析日期、失败不中断、重跑不重复。文件名格式是MCD19A2.AYYYYDDD.hXXvYY.CCC.时间戳.hdf其中AYYYYDDD是年 年积日DOY不是年月日。这里错过一次就够记一辈子。from pathlib import Path import re, datetime as dt PAT re.compile(rMCD19A2\.A(\d{4})(\d{3})\.(h\d{2}v\d{2})\.(\d{3})\.) def parse_name(p: Path): m PAT.search(p.name) if not m: return None year, doy, tile, ver m.groups() day dt.date(int(year), 1, 1) dt.timedelta(daysint(doy) - 1) return {date: day, tile: tile, version: ver, path: p} def run_batch(root: Path, outdir: Path): outdir.mkdir(parentsTrue, exist_okTrue) for hdf in sorted(root.glob(MCD19A2.*.hdf)): info parse_name(hdf) if info is None: continue out outdir / f{info[tile]}_{info[date]:%Y%m%d}.tif if out.exists() and out.stat().st_size 4096: continue # 断点续跑 try: process_one(hdf, out) # 读取 掩膜 重投影 except Exception as e: (outdir / failed.log).open(a).write(f{hdf.name}\t{e}\n)dt.date(year, 1, 1) timedelta(daysdoy-1)这个写法比strptime更稳因为 Python 的%j在某些版本上和年份组合会有坑。out.stat().st_size 4096这个判断是为了防半截文件——程序被 kill 掉时 GeoTIFF 可能已经创建但没写完只判exists()会把坏文件当成已完成后面统计全是错的。失败写日志而不是直接抛出是为了让一个坏文件不至于毁掉整晚的任务。CC那一段是集合版本号6.1 版本写作061。不同集合的scale_factor和 QA 位定义可能有调整混着处理之前先确认版本一致别把 006 和 061 的文件放进同一个目录。4. QA 位掩膜与缺测把 AOD 从“能读”做到“能用”到这一步AOD 已经是带经纬度的浮点栅格了但直接拿去算均值结论一定是错的。原因有三层云和雪没去掉、填充值没掩干净、极端负值没剔除。这一章就是把这三层剥完让剩下的像元真的能代表大气状态。4.1 用位运算把云掩膜和置信度拆出来AOD_QA是 uint16按 bit 段打包。最常被用到的一组是低 8 位import numpy as np qa f.select(AOD_QA).get().astype(uint16) cloud qa 0b0000000000000111 # bit 0-2云状态 adj (qa 3) 0b1 # bit 3是否紧邻云 surf (qa 4) 0b11 # bit 4-5地表类型 snow (qa 6) 0b11 # bit 6-7雪/冰 clear (cloud 0) (adj 0) (snow 0)位运算的关键是掩码和移位要配对。取 bit 0-2掩码就是0b111取完不用移位取 bit 4-5先右移 4 位再和0b11与顺序反了得到的是原位置的位值完全不对。这类错误不会报错只会让筛选结果看起来“有点奇怪”。具体哪个值代表“晴空”要以文件配套的 ATBD 为准别照抄别人的代码。经验上cloud 0是最严格的一档会用掉相当比例的像元如果能接受稍宽松的条件cloud 1能多留一些样本代价是混入薄云。做趋势分析时筛选条件必须全程一致中途换阈值会让时间序列出现台阶。注意不要用AOD_QA 0做筛选。这个条件实际含义是“所有位段全为 0”等于同时要求地表类型是第一种、无雪、无邻云过于苛刻多数区域会筛成空图。4.2 填充值、物理范围、不确定度三重过滤QA 过了还要过数值关。三层过滤按顺序叠加# 第一层填充值与 NaN raw f.select(Optical_Depth_055).get() fill f.select(Optical_Depth_055).attributes()[_FillValue] aod raw.astype(float32) aod[raw fill] np.nan aod * 0.001 # 第二层物理范围负值不能一刀切 aod[(aod -0.05) | (aod 5.0)] np.nan # 第三层不确定度超过阈值视为不可信 unc f.select(AOD_Uncertainty).get().astype(float32) * 0.001 unc[unc 0] np.nan aod[unc 0.1] np.nan good clear np.isfinite(aod)第二层的下限-0.05是留了余量的。晴空条件下 AOD 真实值可以接近 0反演噪声会让它略微为负把负值全砍掉会引入正偏差——尤其在清洁背景区这个偏差能占到均值的百分之几。上限 5.0 是极端值兜底超过这个数的多半是云边缘残留或地表反射率估计失败。第三层的AOD_Uncertainty是常被忽略的一层。它的量纲同样是 0.001超过 0.1 意味着反演相对误差已经很大。做定量分析时加上这一层样本量会掉一截但和地基观测的相关性能明显改善。4.3 有效像元统计与阈值设定掩膜做完下一个问题是“这一天的这张图到底能不能用”。做法是统计有效像元占比valid_mask good np.isfinite(aod) valid_ratio valid_mask.sum() / valid_mask.size if valid_ratio 0.05: quality_flag insufficient # 有效率太低建议整景丢弃 elif valid_ratio 0.2: quality_flag marginal else: quality_flag ok mean_aod np.nanmean(np.where(valid_mask, aod, np.nan))阈值的取法没有标准答案但有个原则单景阈值和时间合成阈值要分开设。单景有效率低于 5% 基本没法用直接丢做月合成时要求每个月至少有 8 到 10 天有效才输出否则那个格点上的“月均值”其实只代表两三天的天气不具代表性。另一个容易忽略的点是区域统计时的加权。在正弦投影下每个像元的实际面积相等直接求算术平均是对的一旦重投影到等经纬度高纬像元被拉伸不加cos(lat)权重会把高纬度的值算重。这个坑在做跨纬度区域比如整省、整流域统计时特别明显低纬地区还好过了 40 度就必须加权。5. 避坑与排查五条让批处理静默失败的细节下面这五条都是跑批量任务时真会遇到的共同点是不报错或者报错信息指不到根因。按现象、原因、解决三段来记比看文档快。5.1 环境与驱动装得上不等于跑得动现象一gdal.Open返回None或者报not recognized as a supported file format但文件明明存在、路径也没写错用别的工具能打开。原因GDAL 没编译 HDF4 驱动。这不是文件问题是安装包问题。很多通过 pip 安装的 GDAL 轮子为了减小体积只带了 GeoTIFF、JPEG 这类常用驱动HDF4 这种依赖外部库的格式被裁掉了。判断方法很简单gdalinfo --formats | grep -i hdf4没输出就是没带。解决换成带 HDF4 的发行版或者在 conda 环境里装conda install -c conda-forge gdal一般会连同 HDF4 一起解决。如果必须用 pip就得自己编译编译前确认系统里装了 HDF4 的开发库和hdf4-config。这个排查顺序很重要——先查驱动再查代码能省掉半天。现象二import pyhdf报ImportError: libdfalt.so: cannot open shared object file。原因pyhdf 是 HDF4 的 Python 封装运行时要去加载 HDF4 的动态库。pip 装的时候只装 Python 包不会帮你把底层库拉下来所以 import 阶段就断了。解决同样走 conda-forge 安装它会把hdf4作为依赖一起装上。装完用python -c from pyhdf.SD import SD; print(ok)验一下别等到跑批处理才发现。另外如果系统里同时存在多个 HDF4 版本可能出现版本冲突显式指定LD_LIBRARY_PATH指向 conda 环境的 lib 目录能解决大部分情况。5.2 数据与数值结果不对但程序不报错现象三批量跑完区域平均 AOD 比预期低一大截而且低得很有规律像是被什么拉下去了。原因填充值参与了平均。要么没掩_FillValue要么在gdalwarp时没设-srcnodata导致 -28672 这种哨兵值被当成真实数值参与运算。还有一种变体设了-srcnodata但值写错了比如写成了 65535等于没设。解决填充值一律从文件属性里读不硬编码。用 GDAL 时读band.GetNoDataValue()用 pyhdf 时读_FillValue。写死数值的代码在版本升级后一定会出问题因为不同集合的填充值定义可能调整。现象四处理完的 AOD 栅格和 MOD09 的地表反射率叠在一起整体错开一段距离几百米到一两公里。原因拿标准正弦瓦片的行列号公式去套 MCD19A2。MAIAC 用自己的一套 1km 网格定义原点、分辨率都在文件的StructMetadata.0里和 MODIS 标准瓦片网格不是同一个。按标准公式硬算得到的是“看起来对但差一点”的结果最难排查。解决坐标一律从文件自带的元数据里取或者干脆交给 GDAL——它的 HDF4 驱动会读StructMetadata并生成正确的地理变换。自己算坐标只作为备用路径且必须用 2.3 里那套解析逻辑。验证方法挑一个已知地面站位置看它落在哪个像元再和目视解译的湖岸线、海岸线比对。现象五重投影后图像四周出现一圈宽窄不一的无效值带或者整幅图有锯齿状空洞。原因-r bilinear在源数据含nan时会把它扩散到邻域边缘一圈全是插值出来的无效值另一个原因是-tr设得比源分辨率细太多上采样过程中有效像元和无效像元混在一起。解决AOD 这类含大量缺测的连续量重投影阶段用-r near需要平滑的话放到重投影之后单独做。输出分辨率不要盲目缩小0.01 度对应 1.1km 左右和源分辨率匹配就够了要更细的分辨率得接受插值带来的信息损失并且用-dstnodata明确标记。6. 月度合成与交叉校验一个能跑通的最小验证闭环单日图跑通了不代表整条链路是对的。我会用一个月的数据做一次端到端验证合成、出图、和地基观测对一遍三步都过才算这条管线能用。import numpy as np, rasterio, glob files sorted(glob.glob(out/h25v04_2019*.tif)) stack [] for fp in files: with rasterio.open(fp) as ds: stack.append(ds.read(1, maskedTrue).filled(np.nan)) arr np.stack(stack) # (nday, h, w) cnt np.isfinite(arr).sum(axis0) # 每格有效天数 with np.errstate(invalidignore): mean np.nanmean(arr, axis0) mean np.where(cnt 8, mean, np.nan) # 有效天数不足的不输出 with rasterio.open(files[0]) as ds: meta ds.meta.copy() meta.update(dtypefloat32, count1, nodatanp.nan) with rasterio.open(monthly_mean.tif, w, **meta) as dst: dst.write(mean.astype(float32), 1)cnt 8这个阈值不是拍脑袋定的一个月 30 天要求至少 8 天有效意味着允许云覆盖掉七成同时保证均值不至于被两三天代表。阈值调低会让覆盖变好、代表性变差做趋势分析时宁可留空洞也不要虚假填充。合成完要校验。最省事的办法是拿同区域的太阳光度计观测做散点看斜率和相关系数没有地基数据时退一步用同一天的 MOD04 产品做交叉比对主要看空间分布形态是否一致别太在意绝对值的几个百分点差异——两者反演假设不同。几个实操上会卡住的点掩膜后的rasterio读取一定要带maskedTrue否则nodata会被读成 0而 0 在 AOD 里是合法值混进去会把均值拉低np.nanmean在整列全nan时会返回nan并抛 RuntimeWarning用errstate压掉但别用np.mean代替输出时nodata设成nan后续读取能直接识别。血泪经验是重投影那一步的-srcnodata和-dstnodata两个参数——我见过不止一次因为只设了前者、后者默认成 0结果一整月的合成图里出现大块“AOD 恰好等于 0”的区域看着像清洁大气其实是空洞被填了零。从那以后我每次跑完批处理都会先统计一遍零值像元占比超过 5% 就停下来查参数不往下走。希望帮到你。本文还有配套的精品资源点击获取
