简介这份代码包围绕NASA MCD19A2地表反射率与植被指数产品整理了一套可复用的遥感批处理流程面向具备IDL与Python基础、从事地球科学、环境监测与气候变化研究的数据处理人员主要解决逐景手动处理时效率低的问题。内容涉及重投影、批量融合与镶嵌、重采样、栅格裁剪以及把AOD等栅格值提取到地面点等环节。压缩包共12个文件约11KB以3个pro基于MCTK/IDL的批处理主程序、3个py重投影与裁剪脚本为主辅以4个xml等工程配置体量轻便、按功能分文件组织便于对照修改。已有1248人学习下载在实际科研场景中具备参考价值。读者可从中了解MCTK批量调用的写法、各环节参数设置思路与项目目录结构并在此基础上做二次开发。1. MCD19A2 原始 HDF 的读取难点第一份 MCD19A2 打开时容易误判一个瓦片一个 HDF里面十几个科学数据集看着像能直接拖进 GIS 出图。真转成 TIFF 叠加到影像上才发现气溶胶反演结果铺在正弦投影的独立网格里行列号和经纬度对不上0.001 的尺度因子没乘-28672 的填充值被当成真实 AOD 参与统计一天下来整幅图全是诡异条纹。这份 MCD19A2 数据处理代码解决的就是这段从原始 HDF 到可分析栅格的脏活读层、解尺度、拆 QA 位、重投影、跨轨道拼接、按月合成。它不是插件形态的东西没有安装包也没有界面就是一组能改参数、能挂进批处理调度里的脚本适合做大气校正、颗粒物估算、长时序 AOD 趋势的从业者拿去当骨架改。2. 认清网格与波段先把行列号换成坐标再谈像素值动手写转换脚本之前必须把两件事确认清楚数据落在什么投影网格上每个 SDS 的数值语义是什么。这两件事搞错后面所有代码都在放大错误——投影错了拼图永远有缝语义错了AOD 会出现 0.05 到 6553.5 这种离谱量级。2.1 正弦投影网格瓦片原点与像元中心的换算MCD19A2 用的是全球正弦投影分瓦片存储瓦片编号 h水平0–35、v垂直0–17单瓦片边长 1111950 米1 km 产品每瓦片 1200×1200 像元所以像元边长正好是 1111950/1200 926.625 米。左上角原点固定任何一行一列都能反算出投影坐标这是后续写 GeoTIFF 时构造仿射变换的基础。import numpy as np R 6371007.181 # 正弦投影球体半径mMODIS 网格定义值 TILE 1111950.0 # 单瓦片边长m PIX TILE / 1200.0 # 1 km 产品像元边长 926.625 m X_MIN -20015109.354 # h0 瓦片左边界 Y_MAX 10007554.677 # v0 瓦片上边界 def tile_origin(h, v): 返回瓦片左上角的投影坐标m return X_MIN h * TILE, Y_MAX - v * TILE def rowcol_to_xy(h, v, row, col): 行列号 - 像元中心的投影坐标注意 y 向下递减 x0, y0 tile_origin(h, v) x x0 (col 0.5) * PIX y y0 - (row 0.5) * PIX return x, y if __name__ __main__: print(tile_origin(25, 4)) # 华东一带常用瓦片 print(rowcol_to_xy(25, 4, 0, 0)) # 左上角第一个像元中心这段代码里最容易被忽略的是0.5数组索引指的是像元左上角而地理坐标通常按像元中心记差半个像元在 1 km 尺度上就是 460 米做小区域对比时足够让整幅图偏移一条边。Y_MAX - v * TILE的减号也别写反HDF 里第 0 行在北、第 1199 行在南而 GeoTIFF 的仿射变换里 y 方向步长必须是负值两者对应不上就会出现上下颠倒的图。拿到投影坐标后转经纬度交给 pyproj投影串直接写projsinu lon_00 a6371007.181 b6371007.181 unitsm。不要图省事用 WGS84 椭球参数去套球体半径差 20 多公里跨瓦片拼接时接缝会明显错开。2.2 SDS 清点哪些是反演值哪些是质量标记同一个文件里数值型数据集和质量型数据集的处理方式完全相反前者要乘尺度因子、要去填充值、可以插值后者是位域打包的整型只能按位与、按位右移绝不能线性插值。开工前先用属性把每个 SDS 的scale_factor、_FillValue、valid_range打印出来做成一张对照表贴在脚本顶部比到处翻文档快得多。SDS 名称含义类型尺度因子填充值处理要点Optical_Depth_047470 nm 气溶胶光学厚度整型0.001-28672乘 0.001 后再判有效Optical_Depth_055550 nm 气溶胶光学厚度整型0.001-28672大气校正常用这一层AOD_QA质量与云标记整型无依版本按位解析不插值AOD_Uncertainty反演不确定度整型0.001-28672可与 AOD 配对做加权Column_WV柱水汽整型0.001-28672辅助大气校正AngstromExp_470-780埃指数整型0.001-28672判断气溶胶粒径有两点提醒。第一部分 SDS 会多一个前导维度按轨道记录的角度类数据直接按二维数组索引会报维度错误或者读出乱数据读取前先调.dimensions()看清楚再决定取第 0 个轨道还是做跨轨道合成。第二不同 Collection 版本的 SDS 名单和 QA 位序都可能变脚本里不要把位定义硬编码散落在各处统一收进一个配置字典换版本只改一处。2.3 用 pyhdf 读一层先核对数值分布再往下走读原始文件优先用 pyhdf它对 HDF4 的支持最直接。别急着一次性把整文件读进内存先单个 SDS 试读打印维度、类型、属性再统计有效像元占比和分位数。这一步做完你基本能判断这个瓦片当天是晴天还是整片云。from pyhdf.SD import SD, SDC import numpy as np PATH MCD19A2.A2020150.h25v04.061.2020152034857.hdf hdf SD(PATH, SDC.READ) print(sorted(hdf.datasets().keys())) # 先看文件里到底有什么 sds hdf.select(Optical_Depth_047) raw sds.get() # 整型数组未做任何缩放 attrs sds.attributes() print(sds.dimensions(), raw.dtype, attrs.get(scale_factor), attrs.get(_FillValue)) fill attrs.get(_FillValue, -28672) scale attrs.get(scale_factor, 0.001) valid raw ! fill # 先做掩膜再做缩放 aod np.where(valid, raw.astype(float32) * scale, np.nan) print(有效像元占比 %.2f%% % (100 * valid.mean())) print(分位数, np.nanpercentile(aod, [5, 50, 95])) sds.endaccess() hdf.end()顺序非常关键先掩膜再缩放。反过来写-28672 乘 0.001 变成 -28.672既不是填充值也不是物理量后面拿它算均值、算分位数全是脏数据而且这种错误不会报错只会在几个月后做趋势分析时才暴露。用np.nan而不是 0 表示无效也是同样的道理0 是一个合法的 AOD 值混淆之后整幅图的有效像元比例就失真了。分位数那一行是低成本的健康检查正常的陆地 AOD 中位数大概在 0.1–0.5 量级如果打印出来中位数是 -28 或者 6500说明掩膜或缩放写错了不用往下跑。3. HDF 转 GeoTIFF投影元数据、重采样与跨轨道拼接读通一个 SDS 只是开头真正耗时间的是把它变成带正确空间参考、能被 QGIS 或下游模型直接吃的栅格。这一层做扎实后面所有流程都省事这一层偷懒后面每个环节都要还债。3.1 直接写正弦投影 GeoTIFF能不转就不转很多人拿到 HDF 第一反应是立刻转成 WGS84 经纬度理由是下游好用。但如果下游只是做区域统计、时序分析、与同类栅格叠加保留原始正弦投影更划算不重采样就无损像元面积一致跨瓦片拼接时不用考虑重采样边界的半像元错位。重投影只做一次放在流程最后别在中间环节反复转。from affine import Affine from rasterio.crs import CRS import rasterio SINU CRS.from_proj4( projsinu lon_00 x_00 y_00 a6371007.181 b6371007.181 unitsm no_defs ) PIX 926.625 X_MIN, Y_MAX -20015109.354, 10007554.677 TILE 1111950.0 def write_sinu_tif(out_path, arr, h, v, nodata-9999.0): x0 X_MIN h * TILE y0 Y_MAX - v * TILE transform Affine(PIX, 0.0, x0, 0.0, -PIX, y0) # y 步长为负 profile { driver: GTiff, height: arr.shape[0], width: arr.shape[1], count: 1, dtype: float32, crs: SINU, transform: transform, nodata: nodata, tiled: True, blockxsize: 256, blockysize: 256, compress: deflate, zlevel: 6, } data arr.astype(float32) data[~np.isfinite(data)] nodata # NaN 统一转 nodata with rasterio.open(out_path, w, **profile) as dst: dst.write(data, 1) dst.update_tags(TILE_Hh, TILE_Vv, SOURCEMCD19A2)几个参数值得单独说。tiledTrue加 256 的分块是为了后面按窗口读取时不用整幅解压做区域裁剪和时序采样时差别很明显。compressdeflate对 AOD 这种大片 nodata 的浮点栅格压缩比通常能到 5 倍以上代价是写入稍慢磁盘紧张时值。nodata-9999.0选一个物理上不可能出现的值别用 -28672 直接当 nodata那个数字在原始整型文件里是填充标记混进浮点结果里迟早出问题。3.2 gdalwarp 重采样连续量和分类量不能用同一个算法确实需要经纬度网格时交给 gdalwarp 批处理比自己写重投影稳。核心就三个参数目标参考系、目标分辨率、重采样算法。算法选错是这一环节最常见的翻车点。#!/usr/bin/env bash set -euo pipefail IN$1 # 输入3.1 产出的正弦投影 GeoTIFF OUT$2 # 输出WGS84 经纬度栅格 gdalwarp \ -t_srs EPSG:4326 \ -tr 0.01 0.01 \ -r bilinear \ -srcnodata -9999 -dstnodata -9999 \ -multi -wo NUM_THREADS4 \ -co COMPRESSDEFLATE -co TILEDYES \ -overwrite $IN $OUT-r bilinear只适合 AOD、水汽这类连续物理量换成 QA、云掩膜、土地覆盖这些整型分类层必须改成-r near否则插值会造出根本不存在的分类码比如把清晰0和云3插成 1.5掩膜逻辑直接失效。-tr 0.01 0.01是约 1 km 的经纬度分辨率在低纬度略细、高纬度略粗如果做全国尺度统计更稳妥的做法是重投影到等积投影如 Albers再统计避免高纬像元被拉伸导致面积加权错误。-multi与NUM_THREADS一起用时注意别和进程级并行叠加线程开太多反而抢 IO。3.3 一天多轨道拼接重叠区靠 QA 选优不靠平均一颗卫星一天过境多次同一瓦片里可能有 2–3 条轨道的数据叠在同一块区域。直接相加取平均是最糟的处理方式云边缘的坏反演会把晴空的好值拉低而且平均之后无法追溯来源。正确逻辑是给每条轨道算一个质量分重叠区只保留最高分的那一个。# 1) 把当天同一瓦片的各轨道成果列进清单 ls work/2020150/h25v04/*.tif list_2020150_h25v04.txt # 2) 建 VRT 做快速拼接先看整体覆盖不做最终成果 gdalbuildvrt -srcnodata -9999 -vrtnodata -9999 \ -input_file_list list_2020150_h25v04.txt \ mosaic_2020150_h25v04.vrt # 3) 落成实体文件按最近邻保留首个有效值 gdal_translate -of GTiff -co COMPRESSDEFLATE -co TILEDYES \ mosaic_2020150_h25v04.vrt mosaic_2020150_h25v04.tifgdalbuildvrt默认行为是后面的文件覆盖前面的所以清单顺序就是优先级顺序。我一般会把 QA 最好的轨道排在最前面再配合-srcnodata保证只有有效像元参与覆盖这样重叠区自然取到优质值。想更严格就让脚本先读 QA 层算每条轨道的清晰像元占比按占比排序生成清单几行 Python 的事比事后修补划算。拼接完必须做一次覆盖检查统计 nodata 占比再和单轨道的有效像元率对比。如果拼接后有效像元反而变少多半是-srcnodata没设或设错覆盖逻辑把整块区域盖成了空白。4. 批量化任务切分、并行与断点续跑单瓦片单天跑通只是验证真正的工作量是几年、几十个瓦片、上千个文件。批量脚本写不好跑到第 800 个文件崩掉前面全白干这种亏吃一次就够。4.1 按日期 × 瓦片切任务别按文件切原始目录里的文件名带日期、瓦片号、处理时间戳同一个瓦片同一天可能有多个版本文件。按文件切任务会导致重复计算和结果覆盖顺序不确定。正确做法是先扫一遍目录用正则解析出日期和瓦片号按(日期, 瓦片)聚合每个任务只保留时间戳最新的那个文件。import re, sqlite3, pathlib PAT re.compile(rMCD19A2\.A(\d{7})\.(h\d{2}v\d{2})\.(\d{3})\.(\d)\.hdf$) def build_tasks(root, db_path): best {} for p in pathlib.Path(root).rglob(*.hdf): m PAT.search(p.name) if not m: continue date, tile, ver, stamp m.groups() key (date, tile) if key not in best or stamp best[key][1]: best[key] (str(p), stamp) con sqlite3.connect(db_path) con.execute(CREATE TABLE IF NOT EXISTS tasks( key TEXT PRIMARY KEY, date TEXT, tile TEXT, path TEXT, status TEXT DEFAULT pending, tries INTEGER DEFAULT 0)) con.executemany( INSERT OR REPLACE INTO tasks(key,date,tile,path,status) VALUES(?,?,?,?, pending), [(f{d}_{t}, d, t, p) for (d, t), (p, _) in best.items()]) con.commit() print(任务数, con.execute(SELECT COUNT(*) FROM tasks).fetchone()[0])正则里的(\d{3})是 Collection 号(\d)是处理时间戳按字符串比较取最大即可因为时间戳位数固定。这个聚合逻辑还有一个附带好处能立刻看出数据缺口比如某年 7 月整月缺 h25v05查一下SELECT date, COUNT(*) FROM tasks GROUP BY date HAVING COUNT(*) N就知道比事后发现时序断档再回头补数据强。4.2 并行度IO 密集的活别把核数开满这类处理九成时间花在解压 HDF 和读写磁盘上不是 CPU 计算。进程数开到物理核数往往比开一半还慢因为磁盘队列被打满每个进程都在等 IO。我的经验值是物理核数的一半SSD 上可以到 3/4机械盘上老老实实 4 个进程。from multiprocessing import Pool from osgeo import gdal gdal.UseExceptions() def init_worker(): # 子进程里压住 GDAL 内部缓存和线程避免多进程互相抢 gdal.SetConfigOption(GDAL_CACHEMAX, 256) # 单位 MB gdal.SetConfigOption(GDAL_NUM_THREADS, 1) gdal.SetConfigOption(CPL_TMPDIR, /dev/shm) # 中间文件走内存盘 def run_one(task): key, path, h, v task try: arr read_and_scale(path, h, v) # 第 2、3 章的读层逻辑 write_sinu_tif(fout/{key}.tif, arr, h, v) return key, ok except Exception as e: return key, ffail: {e} if __name__ __main__: tasks load_pending(db_path) # 从 tasks 表取 pending with Pool(processes6, initializerinit_worker) as pool: for key, status in pool.imap_unordered(run_one, tasks, chunksize1): update_status(db_path, key, status) # 逐条回写崩溃可续GDAL_CACHEMAX必须按进程设不设的话每个子进程都想吃默认缓存8 个进程能把内存吃穿。chunksize1看似低效但单个任务耗时不均匀有的瓦片全是云几秒就完有的要重投影几十秒粗分块会让快的进程干等慢的。逐条imap_unordered加回写是断点续跑的前提。4.3 断点续跑状态写进库别靠文件名猜用输出文件是否存在来判断任务是否完成是小脚本阶段最省事、批量化阶段最坑的做法。写到一半被 kill 掉会留下一个残缺的 TIFF下次运行看到文件存在就跳过这个瓦片永远缺数据而且没有任何报错。用数据库记状态加上先写临时文件再原子重命名这个动作问题就解决了import os, sqlite3 def update_status(db_path, key, status): con sqlite3.connect(db_path) con.execute(UPDATE tasks SET status?, triestries1, WHERE key?.replace(, WHERE, WHERE), (status, key)) con.commit(); con.close() def safe_write(out_path, writer, *args): tmp out_path .part # 先写临时文件 writer(tmp, *args) os.replace(tmp, out_path) # 原子替换要么完整要么不存在-- 重跑时只挑失败和未完成的最多重试 3 次 SELECT key, path FROM tasks WHERE status IN (pending, fail) AND tries 3 ORDER BY date, tile;os.replace在同一文件系统内是原子操作不会出现半个文件被当成完整成果的情况。tries字段是为了挡住那些必然失败的坏文件——比如下载不完整的 HDF——无限重试只会把队列卡死。跑完之后看一眼失败清单通常能归纳出两三类原因比逐个调试高效。5. 避坑与排查五条踩过的记录5.1 出图全是一片负值或条纹现象转出来的 TIFF 在 QGIS 里一片黑或者规则条纹统计均值是 -28 左右。原因没做填充值掩膜就直接乘了尺度因子-28672 × 0.001 -28.672 被当成真实 AOD。解决读取后第一件事raw ! _FillValue建掩膜再缩放最后把无效位置统一设成 nodata 再落盘。检查手段是打印 5% 和 95% 分位数正常陆地区域应落在 0–2 之间。5.2 云掩膜做完了有效像元反而更少现象加了 QA 过滤明明是大晴天有效像元比例从 60% 掉到 5%。原因AOD_QA 是位域打包的整型不同 Collection 的位序改过直接照搬别人脚本里的位表等于用错的掩膜在筛数据。解决先把 QA 层所有出现过的取值列出来np.unique(qa, return_countsTrue)对照产品文档的位定义表确认每个字段的起止位再写(qa start) mask。位表收进配置字典换版本只改一处。5.3 拼接后出现明显的直角接缝现象两个瓦片拼起来边界处 AOD 差 0.2 以上形成一条直线。原因一边用-r bilinear重采样过另一边还是原始正弦网格或者两次重投影的目标分辨率、对齐原点不一致像元错开半个。解决全流程只重投影一次所有瓦片用同一套-te和-tr参数拼接前统一gdalinfo检查 CRS、transform、尺寸三项是否完全一致。5.4 图上下颠倒或左右镜像现象瓦片位置对但地物整个翻转和海陆边界完全镜像。原因HDF 第 0 行在北写 GeoTIFF 时仿射变换 y 步长写成了正值或者用了arr[::-1]翻转数组却没改 transform两处翻转抵消或叠加。解决只在一处控制方向——数组保持原始顺序仿射变换里 y 步长固定为-PIX。验证方法是取瓦片四角坐标转经纬度和海陆位置比对。5.5 跑到一半内存爆掉现象单文件测试没问题批量跑第 20 个左右进程被 OOM kill。原因默认按 float64 读一次把整文件所有 SDS 读进内存子进程的 GDAL 缓存又没限制几个进程叠加就把内存吃穿。解决只读需要的 SDS立刻转 float32读一个写一个及时释放子进程初始化里设GDAL_CACHEMAX进程数按物理核数一半起步。做时序统计时用窗口分块读别整幅读。6. 验收QA 直方图 地基观测对表把结果钉死流程跑通不代表性我一般在批量出结果前后各做一次验证两次都不省。第一次是落盘时自动跑的体检指标有效像元比例、AOD 中位数、超过 3.0 的异常像元占比、QA 各取值分布。晴天瓦片有效比例通常能到 40% 以上全月平均低于 15% 就要怀疑是不是碰到了持续的云季或者掩膜写错。import numpy as np def health_check(aod, qa): ok np.isfinite(aod) stats { valid_ratio: float(ok.mean()), median: float(np.nanmedian(aod)), p95: float(np.nanpercentile(aod, 95)), outlier_ratio: float(np.nansum(aod 3.0) / max(ok.sum(), 1)), } # QA 位域先看取值分布再决定阈值 vals, cnts np.unique(qa, return_countsTrue) stats[qa_top] sorted(zip(cnts.tolist(), vals.tolist()), reverseTrue)[:5] return stats第二次验证是把月合成结果和地基观测对表。做法是按站点经纬度取周围 3×3 像元的中位数避开单像元噪声和同期日均值配对算相关系数、RMSE 和偏差。经验上经过 QA 过滤的日尺度 AOD 与地基观测的相关系数能到 0.85 以上RMSE 在 0.1 左右如果相关系数只有 0.5八成是时间匹配没做对——卫星过境是瞬时值地基是小时均值直接拿整天平均去比会引入系统偏差。检查项合理区间异常时的首要怀疑有效像元比例月均 15%–60%QA 位表或云掩膜AOD 中位数0.1–0.5陆地尺度因子或填充值异常值占比3.0 2%云边缘未剔除与地基相关系数 0.85时间窗或空间窗匹配有一个提速细节月合成不要按月逐天重读原始 HDF先把日成果落成按瓦片分目录的 GeoTIFF再用 VRT 或按窗口读取做聚合重复利用率高得多。同一批数据后面还要做季、年合成或者换 QA 阈值重算日成果留着就是后悔药。从那以后我每批数据落盘前都强制走一遍这套体检加对表宁可多花十分钟看直方图也不愿意一周后回头查为什么某年 3 月的 AOD 突然低了 30%。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/11/5 4:03:23 · 阅读 1.2k MCD19A2气溶胶AOD处理：HDF转GeoTIFF与批量化实践
简介这份代码包围绕NASA MCD19A2地表反射率与植被指数产品整理了一套可复用的遥感批处理流程面向具备IDL与Python基础、从事地球科学、环境监测与气候变化研究的数据处理人员主要解决逐景手动处理时效率低的问题。内容涉及重投影、批量融合与镶嵌、重采样、栅格裁剪以及把AOD等栅格值提取到地面点等环节。压缩包共12个文件约11KB以3个pro基于MCTK/IDL的批处理主程序、3个py重投影与裁剪脚本为主辅以4个xml等工程配置体量轻便、按功能分文件组织便于对照修改。已有1248人学习下载在实际科研场景中具备参考价值。读者可从中了解MCTK批量调用的写法、各环节参数设置思路与项目目录结构并在此基础上做二次开发。1. MCD19A2 原始 HDF 的读取难点第一份 MCD19A2 打开时容易误判一个瓦片一个 HDF里面十几个科学数据集看着像能直接拖进 GIS 出图。真转成 TIFF 叠加到影像上才发现气溶胶反演结果铺在正弦投影的独立网格里行列号和经纬度对不上0.001 的尺度因子没乘-28672 的填充值被当成真实 AOD 参与统计一天下来整幅图全是诡异条纹。这份 MCD19A2 数据处理代码解决的就是这段从原始 HDF 到可分析栅格的脏活读层、解尺度、拆 QA 位、重投影、跨轨道拼接、按月合成。它不是插件形态的东西没有安装包也没有界面就是一组能改参数、能挂进批处理调度里的脚本适合做大气校正、颗粒物估算、长时序 AOD 趋势的从业者拿去当骨架改。2. 认清网格与波段先把行列号换成坐标再谈像素值动手写转换脚本之前必须把两件事确认清楚数据落在什么投影网格上每个 SDS 的数值语义是什么。这两件事搞错后面所有代码都在放大错误——投影错了拼图永远有缝语义错了AOD 会出现 0.05 到 6553.5 这种离谱量级。2.1 正弦投影网格瓦片原点与像元中心的换算MCD19A2 用的是全球正弦投影分瓦片存储瓦片编号 h水平0–35、v垂直0–17单瓦片边长 1111950 米1 km 产品每瓦片 1200×1200 像元所以像元边长正好是 1111950/1200 926.625 米。左上角原点固定任何一行一列都能反算出投影坐标这是后续写 GeoTIFF 时构造仿射变换的基础。import numpy as np R 6371007.181 # 正弦投影球体半径mMODIS 网格定义值 TILE 1111950.0 # 单瓦片边长m PIX TILE / 1200.0 # 1 km 产品像元边长 926.625 m X_MIN -20015109.354 # h0 瓦片左边界 Y_MAX 10007554.677 # v0 瓦片上边界 def tile_origin(h, v): 返回瓦片左上角的投影坐标m return X_MIN h * TILE, Y_MAX - v * TILE def rowcol_to_xy(h, v, row, col): 行列号 - 像元中心的投影坐标注意 y 向下递减 x0, y0 tile_origin(h, v) x x0 (col 0.5) * PIX y y0 - (row 0.5) * PIX return x, y if __name__ __main__: print(tile_origin(25, 4)) # 华东一带常用瓦片 print(rowcol_to_xy(25, 4, 0, 0)) # 左上角第一个像元中心这段代码里最容易被忽略的是0.5数组索引指的是像元左上角而地理坐标通常按像元中心记差半个像元在 1 km 尺度上就是 460 米做小区域对比时足够让整幅图偏移一条边。Y_MAX - v * TILE的减号也别写反HDF 里第 0 行在北、第 1199 行在南而 GeoTIFF 的仿射变换里 y 方向步长必须是负值两者对应不上就会出现上下颠倒的图。拿到投影坐标后转经纬度交给 pyproj投影串直接写projsinu lon_00 a6371007.181 b6371007.181 unitsm。不要图省事用 WGS84 椭球参数去套球体半径差 20 多公里跨瓦片拼接时接缝会明显错开。2.2 SDS 清点哪些是反演值哪些是质量标记同一个文件里数值型数据集和质量型数据集的处理方式完全相反前者要乘尺度因子、要去填充值、可以插值后者是位域打包的整型只能按位与、按位右移绝不能线性插值。开工前先用属性把每个 SDS 的scale_factor、_FillValue、valid_range打印出来做成一张对照表贴在脚本顶部比到处翻文档快得多。SDS 名称含义类型尺度因子填充值处理要点Optical_Depth_047470 nm 气溶胶光学厚度整型0.001-28672乘 0.001 后再判有效Optical_Depth_055550 nm 气溶胶光学厚度整型0.001-28672大气校正常用这一层AOD_QA质量与云标记整型无依版本按位解析不插值AOD_Uncertainty反演不确定度整型0.001-28672可与 AOD 配对做加权Column_WV柱水汽整型0.001-28672辅助大气校正AngstromExp_470-780埃指数整型0.001-28672判断气溶胶粒径有两点提醒。第一部分 SDS 会多一个前导维度按轨道记录的角度类数据直接按二维数组索引会报维度错误或者读出乱数据读取前先调.dimensions()看清楚再决定取第 0 个轨道还是做跨轨道合成。第二不同 Collection 版本的 SDS 名单和 QA 位序都可能变脚本里不要把位定义硬编码散落在各处统一收进一个配置字典换版本只改一处。2.3 用 pyhdf 读一层先核对数值分布再往下走读原始文件优先用 pyhdf它对 HDF4 的支持最直接。别急着一次性把整文件读进内存先单个 SDS 试读打印维度、类型、属性再统计有效像元占比和分位数。这一步做完你基本能判断这个瓦片当天是晴天还是整片云。from pyhdf.SD import SD, SDC import numpy as np PATH MCD19A2.A2020150.h25v04.061.2020152034857.hdf hdf SD(PATH, SDC.READ) print(sorted(hdf.datasets().keys())) # 先看文件里到底有什么 sds hdf.select(Optical_Depth_047) raw sds.get() # 整型数组未做任何缩放 attrs sds.attributes() print(sds.dimensions(), raw.dtype, attrs.get(scale_factor), attrs.get(_FillValue)) fill attrs.get(_FillValue, -28672) scale attrs.get(scale_factor, 0.001) valid raw ! fill # 先做掩膜再做缩放 aod np.where(valid, raw.astype(float32) * scale, np.nan) print(有效像元占比 %.2f%% % (100 * valid.mean())) print(分位数, np.nanpercentile(aod, [5, 50, 95])) sds.endaccess() hdf.end()顺序非常关键先掩膜再缩放。反过来写-28672 乘 0.001 变成 -28.672既不是填充值也不是物理量后面拿它算均值、算分位数全是脏数据而且这种错误不会报错只会在几个月后做趋势分析时才暴露。用np.nan而不是 0 表示无效也是同样的道理0 是一个合法的 AOD 值混淆之后整幅图的有效像元比例就失真了。分位数那一行是低成本的健康检查正常的陆地 AOD 中位数大概在 0.1–0.5 量级如果打印出来中位数是 -28 或者 6500说明掩膜或缩放写错了不用往下跑。3. HDF 转 GeoTIFF投影元数据、重采样与跨轨道拼接读通一个 SDS 只是开头真正耗时间的是把它变成带正确空间参考、能被 QGIS 或下游模型直接吃的栅格。这一层做扎实后面所有流程都省事这一层偷懒后面每个环节都要还债。3.1 直接写正弦投影 GeoTIFF能不转就不转很多人拿到 HDF 第一反应是立刻转成 WGS84 经纬度理由是下游好用。但如果下游只是做区域统计、时序分析、与同类栅格叠加保留原始正弦投影更划算不重采样就无损像元面积一致跨瓦片拼接时不用考虑重采样边界的半像元错位。重投影只做一次放在流程最后别在中间环节反复转。from affine import Affine from rasterio.crs import CRS import rasterio SINU CRS.from_proj4( projsinu lon_00 x_00 y_00 a6371007.181 b6371007.181 unitsm no_defs ) PIX 926.625 X_MIN, Y_MAX -20015109.354, 10007554.677 TILE 1111950.0 def write_sinu_tif(out_path, arr, h, v, nodata-9999.0): x0 X_MIN h * TILE y0 Y_MAX - v * TILE transform Affine(PIX, 0.0, x0, 0.0, -PIX, y0) # y 步长为负 profile { driver: GTiff, height: arr.shape[0], width: arr.shape[1], count: 1, dtype: float32, crs: SINU, transform: transform, nodata: nodata, tiled: True, blockxsize: 256, blockysize: 256, compress: deflate, zlevel: 6, } data arr.astype(float32) data[~np.isfinite(data)] nodata # NaN 统一转 nodata with rasterio.open(out_path, w, **profile) as dst: dst.write(data, 1) dst.update_tags(TILE_Hh, TILE_Vv, SOURCEMCD19A2)几个参数值得单独说。tiledTrue加 256 的分块是为了后面按窗口读取时不用整幅解压做区域裁剪和时序采样时差别很明显。compressdeflate对 AOD 这种大片 nodata 的浮点栅格压缩比通常能到 5 倍以上代价是写入稍慢磁盘紧张时值。nodata-9999.0选一个物理上不可能出现的值别用 -28672 直接当 nodata那个数字在原始整型文件里是填充标记混进浮点结果里迟早出问题。3.2 gdalwarp 重采样连续量和分类量不能用同一个算法确实需要经纬度网格时交给 gdalwarp 批处理比自己写重投影稳。核心就三个参数目标参考系、目标分辨率、重采样算法。算法选错是这一环节最常见的翻车点。#!/usr/bin/env bash set -euo pipefail IN$1 # 输入3.1 产出的正弦投影 GeoTIFF OUT$2 # 输出WGS84 经纬度栅格 gdalwarp \ -t_srs EPSG:4326 \ -tr 0.01 0.01 \ -r bilinear \ -srcnodata -9999 -dstnodata -9999 \ -multi -wo NUM_THREADS4 \ -co COMPRESSDEFLATE -co TILEDYES \ -overwrite $IN $OUT-r bilinear只适合 AOD、水汽这类连续物理量换成 QA、云掩膜、土地覆盖这些整型分类层必须改成-r near否则插值会造出根本不存在的分类码比如把清晰0和云3插成 1.5掩膜逻辑直接失效。-tr 0.01 0.01是约 1 km 的经纬度分辨率在低纬度略细、高纬度略粗如果做全国尺度统计更稳妥的做法是重投影到等积投影如 Albers再统计避免高纬像元被拉伸导致面积加权错误。-multi与NUM_THREADS一起用时注意别和进程级并行叠加线程开太多反而抢 IO。3.3 一天多轨道拼接重叠区靠 QA 选优不靠平均一颗卫星一天过境多次同一瓦片里可能有 2–3 条轨道的数据叠在同一块区域。直接相加取平均是最糟的处理方式云边缘的坏反演会把晴空的好值拉低而且平均之后无法追溯来源。正确逻辑是给每条轨道算一个质量分重叠区只保留最高分的那一个。# 1) 把当天同一瓦片的各轨道成果列进清单 ls work/2020150/h25v04/*.tif list_2020150_h25v04.txt # 2) 建 VRT 做快速拼接先看整体覆盖不做最终成果 gdalbuildvrt -srcnodata -9999 -vrtnodata -9999 \ -input_file_list list_2020150_h25v04.txt \ mosaic_2020150_h25v04.vrt # 3) 落成实体文件按最近邻保留首个有效值 gdal_translate -of GTiff -co COMPRESSDEFLATE -co TILEDYES \ mosaic_2020150_h25v04.vrt mosaic_2020150_h25v04.tifgdalbuildvrt默认行为是后面的文件覆盖前面的所以清单顺序就是优先级顺序。我一般会把 QA 最好的轨道排在最前面再配合-srcnodata保证只有有效像元参与覆盖这样重叠区自然取到优质值。想更严格就让脚本先读 QA 层算每条轨道的清晰像元占比按占比排序生成清单几行 Python 的事比事后修补划算。拼接完必须做一次覆盖检查统计 nodata 占比再和单轨道的有效像元率对比。如果拼接后有效像元反而变少多半是-srcnodata没设或设错覆盖逻辑把整块区域盖成了空白。4. 批量化任务切分、并行与断点续跑单瓦片单天跑通只是验证真正的工作量是几年、几十个瓦片、上千个文件。批量脚本写不好跑到第 800 个文件崩掉前面全白干这种亏吃一次就够。4.1 按日期 × 瓦片切任务别按文件切原始目录里的文件名带日期、瓦片号、处理时间戳同一个瓦片同一天可能有多个版本文件。按文件切任务会导致重复计算和结果覆盖顺序不确定。正确做法是先扫一遍目录用正则解析出日期和瓦片号按(日期, 瓦片)聚合每个任务只保留时间戳最新的那个文件。import re, sqlite3, pathlib PAT re.compile(rMCD19A2\.A(\d{7})\.(h\d{2}v\d{2})\.(\d{3})\.(\d)\.hdf$) def build_tasks(root, db_path): best {} for p in pathlib.Path(root).rglob(*.hdf): m PAT.search(p.name) if not m: continue date, tile, ver, stamp m.groups() key (date, tile) if key not in best or stamp best[key][1]: best[key] (str(p), stamp) con sqlite3.connect(db_path) con.execute(CREATE TABLE IF NOT EXISTS tasks( key TEXT PRIMARY KEY, date TEXT, tile TEXT, path TEXT, status TEXT DEFAULT pending, tries INTEGER DEFAULT 0)) con.executemany( INSERT OR REPLACE INTO tasks(key,date,tile,path,status) VALUES(?,?,?,?, pending), [(f{d}_{t}, d, t, p) for (d, t), (p, _) in best.items()]) con.commit() print(任务数, con.execute(SELECT COUNT(*) FROM tasks).fetchone()[0])正则里的(\d{3})是 Collection 号(\d)是处理时间戳按字符串比较取最大即可因为时间戳位数固定。这个聚合逻辑还有一个附带好处能立刻看出数据缺口比如某年 7 月整月缺 h25v05查一下SELECT date, COUNT(*) FROM tasks GROUP BY date HAVING COUNT(*) N就知道比事后发现时序断档再回头补数据强。4.2 并行度IO 密集的活别把核数开满这类处理九成时间花在解压 HDF 和读写磁盘上不是 CPU 计算。进程数开到物理核数往往比开一半还慢因为磁盘队列被打满每个进程都在等 IO。我的经验值是物理核数的一半SSD 上可以到 3/4机械盘上老老实实 4 个进程。from multiprocessing import Pool from osgeo import gdal gdal.UseExceptions() def init_worker(): # 子进程里压住 GDAL 内部缓存和线程避免多进程互相抢 gdal.SetConfigOption(GDAL_CACHEMAX, 256) # 单位 MB gdal.SetConfigOption(GDAL_NUM_THREADS, 1) gdal.SetConfigOption(CPL_TMPDIR, /dev/shm) # 中间文件走内存盘 def run_one(task): key, path, h, v task try: arr read_and_scale(path, h, v) # 第 2、3 章的读层逻辑 write_sinu_tif(fout/{key}.tif, arr, h, v) return key, ok except Exception as e: return key, ffail: {e} if __name__ __main__: tasks load_pending(db_path) # 从 tasks 表取 pending with Pool(processes6, initializerinit_worker) as pool: for key, status in pool.imap_unordered(run_one, tasks, chunksize1): update_status(db_path, key, status) # 逐条回写崩溃可续GDAL_CACHEMAX必须按进程设不设的话每个子进程都想吃默认缓存8 个进程能把内存吃穿。chunksize1看似低效但单个任务耗时不均匀有的瓦片全是云几秒就完有的要重投影几十秒粗分块会让快的进程干等慢的。逐条imap_unordered加回写是断点续跑的前提。4.3 断点续跑状态写进库别靠文件名猜用输出文件是否存在来判断任务是否完成是小脚本阶段最省事、批量化阶段最坑的做法。写到一半被 kill 掉会留下一个残缺的 TIFF下次运行看到文件存在就跳过这个瓦片永远缺数据而且没有任何报错。用数据库记状态加上先写临时文件再原子重命名这个动作问题就解决了import os, sqlite3 def update_status(db_path, key, status): con sqlite3.connect(db_path) con.execute(UPDATE tasks SET status?, triestries1, WHERE key?.replace(, WHERE, WHERE), (status, key)) con.commit(); con.close() def safe_write(out_path, writer, *args): tmp out_path .part # 先写临时文件 writer(tmp, *args) os.replace(tmp, out_path) # 原子替换要么完整要么不存在-- 重跑时只挑失败和未完成的最多重试 3 次 SELECT key, path FROM tasks WHERE status IN (pending, fail) AND tries 3 ORDER BY date, tile;os.replace在同一文件系统内是原子操作不会出现半个文件被当成完整成果的情况。tries字段是为了挡住那些必然失败的坏文件——比如下载不完整的 HDF——无限重试只会把队列卡死。跑完之后看一眼失败清单通常能归纳出两三类原因比逐个调试高效。5. 避坑与排查五条踩过的记录5.1 出图全是一片负值或条纹现象转出来的 TIFF 在 QGIS 里一片黑或者规则条纹统计均值是 -28 左右。原因没做填充值掩膜就直接乘了尺度因子-28672 × 0.001 -28.672 被当成真实 AOD。解决读取后第一件事raw ! _FillValue建掩膜再缩放最后把无效位置统一设成 nodata 再落盘。检查手段是打印 5% 和 95% 分位数正常陆地区域应落在 0–2 之间。5.2 云掩膜做完了有效像元反而更少现象加了 QA 过滤明明是大晴天有效像元比例从 60% 掉到 5%。原因AOD_QA 是位域打包的整型不同 Collection 的位序改过直接照搬别人脚本里的位表等于用错的掩膜在筛数据。解决先把 QA 层所有出现过的取值列出来np.unique(qa, return_countsTrue)对照产品文档的位定义表确认每个字段的起止位再写(qa start) mask。位表收进配置字典换版本只改一处。5.3 拼接后出现明显的直角接缝现象两个瓦片拼起来边界处 AOD 差 0.2 以上形成一条直线。原因一边用-r bilinear重采样过另一边还是原始正弦网格或者两次重投影的目标分辨率、对齐原点不一致像元错开半个。解决全流程只重投影一次所有瓦片用同一套-te和-tr参数拼接前统一gdalinfo检查 CRS、transform、尺寸三项是否完全一致。5.4 图上下颠倒或左右镜像现象瓦片位置对但地物整个翻转和海陆边界完全镜像。原因HDF 第 0 行在北写 GeoTIFF 时仿射变换 y 步长写成了正值或者用了arr[::-1]翻转数组却没改 transform两处翻转抵消或叠加。解决只在一处控制方向——数组保持原始顺序仿射变换里 y 步长固定为-PIX。验证方法是取瓦片四角坐标转经纬度和海陆位置比对。5.5 跑到一半内存爆掉现象单文件测试没问题批量跑第 20 个左右进程被 OOM kill。原因默认按 float64 读一次把整文件所有 SDS 读进内存子进程的 GDAL 缓存又没限制几个进程叠加就把内存吃穿。解决只读需要的 SDS立刻转 float32读一个写一个及时释放子进程初始化里设GDAL_CACHEMAX进程数按物理核数一半起步。做时序统计时用窗口分块读别整幅读。6. 验收QA 直方图 地基观测对表把结果钉死流程跑通不代表性我一般在批量出结果前后各做一次验证两次都不省。第一次是落盘时自动跑的体检指标有效像元比例、AOD 中位数、超过 3.0 的异常像元占比、QA 各取值分布。晴天瓦片有效比例通常能到 40% 以上全月平均低于 15% 就要怀疑是不是碰到了持续的云季或者掩膜写错。import numpy as np def health_check(aod, qa): ok np.isfinite(aod) stats { valid_ratio: float(ok.mean()), median: float(np.nanmedian(aod)), p95: float(np.nanpercentile(aod, 95)), outlier_ratio: float(np.nansum(aod 3.0) / max(ok.sum(), 1)), } # QA 位域先看取值分布再决定阈值 vals, cnts np.unique(qa, return_countsTrue) stats[qa_top] sorted(zip(cnts.tolist(), vals.tolist()), reverseTrue)[:5] return stats第二次验证是把月合成结果和地基观测对表。做法是按站点经纬度取周围 3×3 像元的中位数避开单像元噪声和同期日均值配对算相关系数、RMSE 和偏差。经验上经过 QA 过滤的日尺度 AOD 与地基观测的相关系数能到 0.85 以上RMSE 在 0.1 左右如果相关系数只有 0.5八成是时间匹配没做对——卫星过境是瞬时值地基是小时均值直接拿整天平均去比会引入系统偏差。检查项合理区间异常时的首要怀疑有效像元比例月均 15%–60%QA 位表或云掩膜AOD 中位数0.1–0.5陆地尺度因子或填充值异常值占比3.0 2%云边缘未剔除与地基相关系数 0.85时间窗或空间窗匹配有一个提速细节月合成不要按月逐天重读原始 HDF先把日成果落成按瓦片分目录的 GeoTIFF再用 VRT 或按窗口读取做聚合重复利用率高得多。同一批数据后面还要做季、年合成或者换 QA 阈值重算日成果留着就是后悔药。从那以后我每批数据落盘前都强制走一遍这套体检加对表宁可多花十分钟看直方图也不愿意一周后回头查为什么某年 3 月的 AOD 突然低了 30%。希望帮到你。本文还有配套的精品资源点击获取