MCD19A2 AOD批处理：HDF4读取、QA掩膜与投影合成

简介面向遥感、生态监测与气候变化方向的学习者和科研人员这套代码围绕 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(rOptical_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.001MCD19A2 的 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.01fill不要硬编码成全局常量最好从子数据集元数据里读读不到再用常见值兜底。检查时先看 AOD 的 5%、50%、95% 分位数合理范围通常在 -0.1 到 5 之间整层超过 5 或大面积 -0.1 以下先回头查填充值和尺度因子不要急着做合成。3. 投影与网格对齐正弦投影转目标 CRS 的三种落地方式MCD19A2 原生网格是正弦投影全球等积做全球统计很舒服但一旦要叠矢量、裁区域、和别的栅格做像元级运算直接拿正弦投影往 GIS 里一扔错位是常事。投影这一步不是“转一下就行”目标 CRS、分辨率、范围、重采样方式四个参数要一起定否则后面合成时网格对不上前面读得再对也白搭。3.1 正弦投影不是“错”但直接给 GIS 用会挨骂先定目标网格正弦投影本身没问题问题在于它和常见 GIS 图层的坐标网格不一致。区域分析我一般先定目标网格四件事目标 CRS、分辨率、输出范围、重采样方式。做区域统计优先选等面积投影做叠加展示选经纬度或本地投影分辨率不要一上来就放大到 30 米源数据是 1 公里放大只会让像元看起来更细信息量不会变多。下面这张表是我批处理前会先填的配置填完再写代码比边跑边改省事。参数建议值说明目标 CRSEPSG:4326 或等面积投影区域统计优先等面积叠加展示用经纬度分辨率0.01° 或 1 km与源数据匹配不要无脑放大输出范围AOI 外扩半像元避免边缘截断AOD 重采样bilinear连续量双线性比最近邻稳QA 重采样nearest分类量插值会造出小数编码表格里最容易被忽略的是“输出范围外扩半像元”。裁剪时如果范围卡得太死边缘像元会被切掉一半统计时有效像元比例偏低看起来像数据缺失其实是网格没对齐。另一个是 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_srsEPSG:4326, xres0.01, yres0.01): opts gdal.WarpOptions( formatGTiff, dstSRSdst_srs, xResxres, yResyres, resampleAlgbilinear, srcNodata-28672, dstNodata-9999, outputTypegdal.GDT_Float32, multithreadTrue, warpMemoryLimit512, creationOptions[COMPRESSLZW, TILEDYES, BIGTIFFIF_SAFER] ) ds gdal.Warp(dst_path, src_path, optionsopts) ds None warp_aod(aod_047_sub.tif, aod_047_wgs84.tif)逻辑上gdal.Warp先按srcNodata把源填充值排除再按resampleAlg重采样最后用dstNodata填输出空洞。参数上dstSRS填目标 CRSxRes/yRes填目标分辨率srcNodata必须和源数据一致dstNodata用 -9999 这类明显偏离真实范围的值别用 0。warpMemoryLimit512是 MB批量跑大文件时别设太小否则会频繁落临时文件multithreadTrue对多核机器有用但内存吃紧时先关掉。需要固定输出范围时再加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, cropTrue, all_touchedFalse, filledTrue, 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按矢量几何裁栅格cropTrue会缩小输出范围filledTrue用nodata填洞。参数上all_touchedFalse适合均值统计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(rMCD19A2\.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(daysdoy - 1) files sorted(Path(data).glob(MCD19A2.A*.hdf), keyparse_date) for f in files[:3]: print(parse_date(f), f.name)逻辑上正则抓A后面的 4 位年份和 3 位年积日再用datetime(year, 1, 1) timedelta(daysdoy-1)转成日期。参数上doy-1不能少年积日从 1 开始keyparse_date让排序按真实日期走不按文件名字典序。多 tile 场景下glob出来会混着不同 tile合成前按文件名里的 tile 字段再分组否则同一网格会叠进不同区域的数据。4.2 QA 加权合成晴天多给权重云边缘别硬留日合成不是把所有有效像元一视同仁。晴空像元权重高可能晴空权重低云边缘要么给低权重要么直接剔除。下面这段把 QA 位段转成权重再做加权平均权重为 0 的像元不参与分母避免除零。import numpy as np def qa_weight(qa, cloud_offset0, cloud_width3): cloud (qa cloud_offset) ((1 cloud_width) - 1) w np.zeros_like(qa, dtypenp.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, dtypenp.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, axis0) den np.nansum(w, axis0) 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_valid5): # cube: (time, y, x) s pd.Series(range(len(dates)), indexpd.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(axis0) mean np.nanmean(sub, axis0) mean[valid min_valid] np.nan out.append((period, mean)) return out逻辑上pd.Series把时间索引和数组位置绑在一起groupby按年月拿到位置索引再切cube。参数上min_valid5是最低有效天数区域小、云多的时候可以降到 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 单独走resampleAlgnearest重投影到和目标 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 位段取错了或者权重没写进合成。检查项异常信号优先排查有效像元比例低于 10%云掩膜、裁剪范围、填充值AOD 分位数P95 超过 5 或 P50 贴近 0尺度因子、NoData 是否转 NaNQA 分层均值晴空与云边缘差异很小QA 位偏移、权重数组我现在的习惯是每批结果先跑quick_check再看 QA 分层表不通过就不回写月均值。以前图省事跳过这一步后面返工把投影、裁剪、合成全重跑了一遍那个后悔药真不好吃。希望帮到你。本文还有配套的精品资源点击获取