简介面向遥感与地理信息领域的初学者及科研应用人员这份文档围绕MODIS数据的背景、下载与处理展开系统梳理可帮助缺乏卫星遥感经验的读者快速建立从数据获取到分析应用的完整认知。包内仅含1个doc文档压缩包约503KB篇幅紧凑却覆盖了ESE与EOS计划脉络、Terra与Aqua双星互补观测机制、36个波段的光谱范围与空间分辨率、2330公里视场及705公里太阳同步极轨等关键技术指标并给出波段设计目标与典型应用领域的对照。文档同时说明NASA官网的产品与时间范围选取流程、HDF格式文件的读取方式以及借助ENVI、ArcGIS或开源工具完成预处理、辐射与几何校正、时序合成和植被指数、地表温度提取的思路。已有304人学习适合作为课程作业、课题入门或环境变化监测项目的前期参考资料。1. MODIS数据1到2天覆盖全球的36波段到底适合谁如果你要把一个区域几年的植被长势、地表温度或者云量变化跑一遍又不想在数据采购上耗预算MODIS数据大概率是最先该试的那一套。Terra 和 Aqua 两颗太阳同步极轨卫星互相补位上午星和下午星各自过境1 到 2 天能把全球表面扫一遍星上搭载的中分辨率成像光谱仪有 36 个波段光谱范围从 0.405 微米一直铺到 14.385 微米空间分辨率分 250 米、500 米、1000 米三档幅宽 2330 公里。它做不了米级精细制图也替代不了高光谱专题产品但做区域趋势、长时序对比、大范围快速摸底这套数据的性价比很难绕开。下面按产品选型、下载、HDF4 读取、MRT 投影、排错、合成验证的顺序把这份资料里能直接落地的部分拆开讲。2. 产品分级与文件名MOD、MYD、MCD三个前缀背后的选型账2.1 0级到4级先分清你要的是原始观测还是算好的产品MODIS 标准数据产品按处理深度分成 0 到 4 级这个分级不是学术分类它直接决定你下载的文件里装的是什么、后面要不要自己写算法。0 级是地面站直接收到的原始数据一般用户碰不到1A 级做了数据重建和时间配准仍然是没有辐射校正的原始计数值1B 级在 1A 基础上做了辐射校正和定标输出的是辐亮度MOD02 就是这一级。真正开始“能用”的是 2 级以后2 级产品和 1B 具有相同空间分辨率、覆盖同一地理区域比如 MOD09 表面反射率、MOD11 地表温度3 级是把多天数据放进统一时空栅格做合成比如 MOD13 的 16 天 NDVI4 级则是跑完模型之后的结果比如 MOD17 的初级生产力。选型的判断依据很简单如果你要自己做指数计算、自己控制云掩膜就从 2 级甚至 1B 起步如果你只想要一条能直接分析的时间序列3 级产品省掉的预处理时间远超它损失的那点灵活性。我一般会先确认研究问题的时间尺度——日尺度看 2 级旬月尺度看 3 级年度碳收支看 4 级别一上来就下 1B几十景原始辐亮度堆在硬盘里光定标和几何定位就能拖掉一周。2.2 44种标准产品怎么倒着选标准产品按内容分成大气、陆地、海洋三大类合计 44 种。新手最容易犯的错是照着产品编号从 MOD01 往下翻翻到哪算哪。更靠谱的做法是从应用倒推先写清楚你要的输出是什么再找对应的短名。下面这张表是日常用得最多的十几个覆盖了土地覆被、植被、温度、火情、云和气溶胶几个方向。短名内容空间分辨率时间频率典型用途MOD021B 辐亮度250/500/1000 m每轨自定义算法输入MOD03地理定位1 km每轨1B 几何校正必备MOD04气溶胶1 km10 km 产品日大气颗粒物研究MOD06云产品1 km日云量、云顶高度MOD09表面反射率250/500 m日指数计算、分类MOD11地表温度与发射率1 km日/8 天城市热岛、干旱MOD13NDVI/EVI250 m16 天植被长势MOD14热异常与火情1 km日火点监测MOD15LAI/FPAR1 km8 天生态模型输入MOD17GPP/NPP1 km8 天/年碳循环MOD35云掩膜250 m/1 km日筛云MCD12Q1土地覆被类型500 m年下垫面分类表里没有列的去产品手册里按“短名 版本号”检索不要凭记忆猜。版本号很关键同一产品不同 Collection 的处理基线不同定标系数和 QA 位定义都可能变混用会得到一条斜率对不上的时间序列。2.3 文件名拆成七段每一段都影响后面脚本MODIS 的文件名不是随手拼的拿MOD01.A2002032.1015.003.2002034184356.hdf举例七段各有含义写批处理脚本时经常要按段提取。段位示例含义脚本里的用法aMOD01产品短名MOD 为 TerraMYD 为 AquaMCD 为两者合成判断平台决定后续用哪套参数bA2002数据采集年份按年建目录c032儒略日第 32 天即 2 月 1 日转成日期和时间序列对齐d1015UTC 时分10 点 15 分判断是上午星还是下午星过境e003Collection 版本号同一批次必须统一f2002034184356生产时间年儒略日时分秒排查“同一天两个版本”时用ghdf扩展名全部为 HDF 格式决定用 HDF4 还是 HDF5 库打开最容易翻车的是 c 段和 d 段。儒略日不是月日032是 2 月 1 日不是 3 月 2 日1015是 UTC 时间换算到北京时间要加 8 小时所以一景标着A2020032.1015的数据实际对应北京时间 2 月 1 日 18 点 15 分。做日尺度统计时如果直接拿文件名里的日期当本地日期跨年那几天必错。2.4 正弦投影与h17v04瓦片编号不是随便起的大多数标准产品用正弦投影Sinusoidal赤道处按 10 度乘 10 度切网格行号从左上角 (0,0) 起算到右下角 (35,17)全球共 36 列 18 行。文件里的h17v04就是水平第 17 列、垂直第 4 行一景大致覆盖中纬度一个省级尺度的范围。瓦片编号的价值在于同一h17v04的多年数据可以直接按像元对齐做时间序列不需要先做空间重采样跨瓦片研究才需要拼接而拼接一定伴随投影和像元网格的变化。这里有个常见误解以为同一瓦片号在所有产品里范围完全一致。250 米、500 米、1 公里产品的瓦片覆盖范围是嵌套的但像元数不同MOD09GQ250 米和 MOD09GA500 米/1 公里的阵列尺寸差一倍混着读之前先看清楚SDS的维度。投影参数、瓦片号、分辨率这三个信息凑齐才能判断两景数据能不能直接逐像元运算。3. 下载从经纬度检索到wget批量拉取的四个步骤3.1 检索条件写成“经纬度框UTC时间段”数据检索页通常提供按经纬度框选、按瓦片号选、按日期范围选三种方式。区域研究用经纬度框最直接输入左上角和右下角坐标再填时间范围。时间一定按 UTC 填要取北京时间 1 月 16 日全天的数据UTC 范围是 1 月 15 日 16:00 到 1 月 16 日 16:00直接填 1 月 16 日 00:00 到 24:00 会漏掉前 8 小时。检索结果里会列出每景的文件名、大小、平台、版本先别急着全选用文件名里的 d 段确认过境时间是不是你要的上午星或下午星。检索结果一般支持导出成文本清单这个清单就是后面批量下载的输入。清单里通常带完整下载链接但链接可能带会话参数过几小时就失效所以脚本里更稳的做法是只保留文件名把数据目录作为变量拼前缀。L1B 数据记得把 MOD03 一起勾上它是 1 公里、500 米、250 米产品做几何纠正的定位文件漏下这一景后面打开 MOD02 会发现没有经纬度返工重下最费时间。3.2 账号与token把登录这一步提前做掉数据平台现在基本都走统一登录匿名 FTP 的老路子已经不好用。注册账号后在个人中心生成一个下载 token之后所有请求带Authorization: Bearer头即可不用每次输密码也不会因为会话过期在半夜断掉。token 等同于密码别写进公开仓库脚本里用环境变量传入。下面这段是批量下载的骨架把清单文件、数据目录、token 三个变量换掉就能跑。#!/usr/bin/env bash # 批量下载 MODIS 数据断点续传、失败重试、限速 set -u export ED_TOKEN你的下载token BASE_URLhttps://检索页复制下来的数据目录 OUT_DIR./modis_raw mkdir -p $OUT_DIR while read -r name; do [ -z $name ] continue # 跳过空行 [ -f $OUT_DIR/$name ] continue # 已存在则跳过 wget -c -nc --tries5 --timeout60 --wait2 \ --headerAuthorization: Bearer ${ED_TOKEN} \ -P $OUT_DIR ${BASE_URL}/${name} || echo FAILED $name failed.txt done urls.txt-c是断点续传网络中断后重跑不会从头拉-nc防止覆盖已完成的文件--tries5和--timeout60对付偶发的连接超时--wait2每次请求间隔两秒批量几百景时能明显降低被限流的概率。失败的文件名单独记到failed.txt第二天换时段重跑比盯着终端一屏屏翻要省心。数据目录不要手拼直接从检索结果里复制一层不同产品、不同版本的目录结构不完全一样。3.3 下载完先校验别等读的时候才发现是HTML下载中断最隐蔽的结果不是文件缺失而是文件存在但只有几 MB打开一看是登录失效后返回的 HTML 错误页。HDF4 文件有固定的文件头正常文件开头是 HDF 标识用几行 Python 就能批量筛一遍顺便统计每个文件的 SDS 数量数量为零的直接标红。from pathlib import Path from pyhdf.SD import SD, SDC bad [] for p in Path(modis_raw).glob(*.hdf): try: h SD(str(p), SDC.READ) n len(h.datasets()) # SDS 数量 h.end() if n 0: bad.append((p.name, no SDS)) except Exception as e: bad.append((p.name, repr(e))) for name, reason in bad: print(name, reason) print(待重下:, len(bad))SD打开失败基本就是文件截断或根本不是 HDF能打开但datasets()为空多半是下到了错误页。这个检查跑一遍几十秒比后面在重投影环节报一个看不懂的错要划算。顺手把文件按“产品短名/年/儒略日”建目录后面 MRT 批处理和质量筛选都按这个层级遍历省得再写一套路径解析。3.4 清单、日志、版本三样东西留好批量下载最怕两周后回头找“当时到底下了哪几景”。我的习惯是每批数据在目录下留三个文件原始清单urls.txt、下载日志failed.txt、以及一个manifest.csv记录文件名、大小、下载时间、Collection 版本。看着繁琐等审稿人或合作者问“你这批数据版本是不是统一”的时候这份清单就是后悔药。尤其是 MOD 和 MYD 混用、或者跨 Collection 拼接时版本不一致带来的反射率偏差不会报错只会让趋势线的斜率悄悄偏掉。4. HDF4结构与MRT重投影SDS读取到Albers输出的完整链路4.1 HDF4的层级SDS、Vdata、Vgroup与属性MODIS 的.hdf是 HDF4 格式不是 HDF5两者不兼容。HDF4 内部是层次结构最常用的是 SDSScientific Dataset也就是一维到多维的数组反射率、温度、QA 都是 SDSVdata 存表格型数据Vgroup 负责把相关数据集组织在一起每个 SDS 自带属性定标系数、填充值、有效范围都藏在属性里。用pyhdf.SD打开文件后先列一遍数据集名和维度能避免“凭记忆 select 一个不存在的名字”这种低级错误。from pyhdf.SD import SD, SDC import numpy as np f MOD09GA.A2020032.h17v04.006.2020040034512.hdf hdf SD(f, SDC.READ) for name, info in hdf.datasets().items(): print(name, info[0], info[1:]) # 名称、类型、维度输出里会看到sur_refl_b01到sur_refl_b07、state_1km、solar_zenith、satellite_zenith等。info[0]是数据类型info[1:]是各维长度。MOD09GA 的 500 米反射率维度大约是 1200×12001 公里 QA 是 600×600两者不能直接逐像元运算要么先重采样要么用 500 米波段自己的 QA。读之前先看清楚维度比读完之后发现形状对不上再回头查要快得多。4.2 读反射率scale_factor和_FillValue必须一起处理MODIS 的反射率 SDS 存的是整型不是物理值直接读出来是几千上万的整数。每个 SDS 的属性里有scale_factor、add_offset、_FillValue、valid_range正确顺序是先按填充值掩膜再乘比例加偏移。顺序反了填充值也会被定标成一个看似合理的负值后面统计最小值时就会莫名其妙出现 -0.1 这种数。sds hdf.select(sur_refl_b01) raw sds.get().astype(float32) attr sds.attributes() scale attr.get(scale_factor, 1.0) offset attr.get(add_offset, 0.0) fill attr.get(_FillValue, -28672) refl raw * scale offset # 定标 refl[raw fill] np.nan # 先掩膜再统计 print(refl.shape, np.nanmin(refl), np.nanmax(refl)) sds.end()MOD09 的scale_factor通常是 0.0001_FillValue是 -28672valid_range在 -100 到 16000 之间。算完的反射率理论值域在 0 到 1 附近出现大于 1 的值先别慌看看是不是云边缘或观测几何极端的情况再检查有没有把valid_range之外的像元算进来。对于 MOD11 的地表温度scale_factor是 0.02填充值是 0单位是开尔文要转摄氏度还得减 273.15这些都在产品手册的属性表里别靠猜。4.3 geo_field与经纬度查找表几何定位的两种来源1B 数据本身不带完整的逐像元经纬度MOD03 才是定位文件里面有Latitude、Longitude、SolarZenith、SensorZenith等数据集维度与 1 公里波段一致。做几何校正时可以把 MOD03 的经纬度读出来按行列号对应到 MOD02 的像元上也可以在 MRT 重投影时直接把 MOD03 作为定位输入。两种方式的区别在于前者适合自己写校正算法、做精确的像元级配准后者适合批量出图MRT 内部会处理定位速度快但可控参数少。读 MOD03 的经纬度时注意维度顺序HDF4 的数组是行优先Python 读出来是[row, col]和影像的上下左右方向一致但有些产品会做翻转或偏移拿一景有明显地物特征的数据比如海岸线叠一下就能确认。经纬度是浮点数组文件体积不小批量处理时不要一次性把几十景的经纬度都读进内存按瓦片循环、读完即释放。4.4 MRT先拼接再投影Albers参数一次配好MRTMODIS Reprojection Tool负责重投影、拼接、重采样、裁剪和波段子集。它的处理顺序是先把同一投影下的多幅影像拼接再做投影转换这个顺序不能反如果先投影再拼接接缝处的像元网格对不齐会出现一条明显的错位带。参数文件.prm里最关键的几项是输出投影类型、标准纬线、中央经线、基准面和重采样方式。国内区域制图常用 Albers 等积圆锥投影双标准纬线 25 度和 47 度、中央经线 105 度是一套通用参数省级范围可以按轮廓再调比如东西跨度大的区域把中央经线往中间挪。参数常用值说明OUTPUT_PROJECTION_TYPEALBERS等积圆锥适合面积统计STDPR125.0第一标准纬线STDPR247.0第二标准纬线CENTRAL_MERIDIAN105.0中央经线按研究区调整ORIGIN_LATITUDE0.0原点纬度DATUMWGS84基准面RESAMPLING_TYPEBILINEAR连续量用双线性分类产品用最近邻OUTPUT_PIXEL_SIZE500输出像元大小单位米重采样方式要按数据类型选反射率、温度这类连续量用双线性土地覆被、QA 这类分类或位标志必须用最近邻否则一个 3 类像元被插成 2.5 类后面统计类别面积全是错的。输出像元大小不要小于原始分辨率500 米产品输出 250 米不会增加信息只会让文件大四倍。参数文件写好后用-R传入文件清单、-O指定输出目录批量跑。MRT_HOME/opt/mrt/bin $MRT_HOME/resample -P albers.prm -R filelist.txt -O ./out_albers -V-V打开详细日志每景处理成功还是失败都会打印批量几百景时把日志重定向到文件失败景单独重跑。MRT 对文件名有要求输入清单里写绝对路径更稳输出文件名会自动加波段后缀后期按规则批量改名即可。跑完先抽一景和原始正弦投影数据叠一下看边界是不是吻合再接下去做时间序列。5. 避坑排查下载、读取、投影五个翻车现场5.1 日期差一天UTC和儒略日叠在一起现象按“1 月 16 日”检索下回来的数据文件名是A2003015做日统计时这一天的值对不上。原因文件名里的日期是数据采集的 UTC 儒略日检索页填的时间范围也是 UTC。北京时间 1 月 16 日 00:00 到 08:00 属于 UTC 1 月 15 日这 8 小时的数据会落在A2003015里。如果只按文件名日期分组就会把本地日期的一整天拆成两半。解决统一在脚本里做时区转换把文件名解析成 UTC 时间后再加 8 小时得到本地日期按本地日期分组或者干脆全程用 UTC 做时间轴出图时再标注本地时间。两种都行关键是别混。5.2 只下MOD02打开没有地理定位现象MOD02 能打开但找不到经纬度数据集做裁剪和重投影时无从下手。原因MOD02 是 1B 辐亮度几何定位信息在配套的 MOD03 文件里两者是分开下载的。检索时只勾了 MOD02MOD03 就漏了。解决检索清单里同时勾选同一天、同一轨的 MOD03下载后按文件名里的时间戳配对。MRT 重投影时把 MOD03 作为定位输入或者在 Python 里读 MOD03 的Latitude、Longitude做自定义校正。配对规则写进脚本别靠人工比对。5.3 用h5py读HDF4报错还找不到原因现象h5py.File(MOD09GA...hdf)直接抛异常提示不是 HDF5 文件。原因MODIS 标准产品是 HDF4h5py只认 HDF5。两个格式同名不同代API 完全不兼容。网上不少示例代码混着写复制过来就报错。解决HDF4 用pyhdfHDF5 用h5py先看文件头再选库。GDAL 两类都能读gdalinfo一跑就知道格式。批量处理前统一检查一遍文件头避免一半文件用一种库、另一半用另一种库。5.4 反射率值域不对scale_factor没乘现象读出来的反射率在 0 到 10000 之间做指数计算全是异常值。原因SDS 存的是整型物理值要乘scale_factor再加add_offset。只读了get()没读属性或者顺序反了先定标再掩膜填充值被算成了有效值。解决固定写成“读原始数组 → 读属性 → 掩膜填充值 → 定标”四步封装成函数复用。每次换产品先打印attributes()看一遍MOD09、MOD11、MOD13 的 scale 和 fill 都不一样不能一套代码硬套所有产品。5.5 MRT输出接缝、错位与内存爆掉现象多景拼接后中间有一条错位带或者跑到一半进程被系统杀掉。原因错位多半是先投影再拼接或者输入清单里混了不同 Collection、不同分辨率的文件内存爆掉通常是一次性把几十景全丢给 MRT或者输出像元设得过细。解决确认 MRT 的处理顺序是先拼接后投影输入清单按产品、版本、分辨率分组不同组分开跑批量时按瓦片或按年分批每批 10 到 20 景跑完检查接缝再继续。输出像元不小于原始分辨率连续量用双线性、分类量用最近邻。跑之前先用两景做小样本测试参数确认无误再放开批量。6. 从单日瓦片到月度合成QA筛云与结果验证的一个笨办法单日反射率直接平均是最容易出问题的做法云、云影、观测几何差异都会混进月度值里。我的习惯是先读state_1km做云筛再按像元统计有效观测数有效数低于阈值比如月内少于 5 次的像元直接置空不参与后续分析。MOD09GA 的state_1km是位标志低两位表示云状态0 为晴空1 为有云2 为混合3 为未定义先掩膜再定标。import numpy as np from pyhdf.SD import SD, SDC def read_clear_refl(path, bandsur_refl_b01): hdf SD(path, SDC.READ) raw hdf.select(band).get().astype(float32) attr hdf.select(band).attributes() qa hdf.select(state_1km).get() hdf.end() refl raw * attr.get(scale_factor, 1.0) attr.get(add_offset, 0.0) refl[raw attr.get(_FillValue, -28672)] np.nan cloud qa 0b11 # 低两位云状态 refl[(cloud ! 0) | (qa 1)] np.nan # 非晴空和填充值都剔除 return refl这段代码只做了一件事把非晴空像元变成NaN后续用np.nanmean做月度合成时自然跳过。跑完一个月后别急着出图先看两个数有效像元比例和月度均值的直方图。有效比例低于 30% 的区域说明这个月云太多合成结果的可信度要打折扣直方图出现双峰或者异常拖尾多半是 QA 位定义用错了版本或者把云影漏掉了。拿一景已知水体或城市核心区的像元值做对照水体反射率在可见光波段应该低而平城市建成区在短波红外明显偏高对不上就回去查定标和掩膜。从那以后我每次下完一批数据都强制先跑一遍“文件头检查 SDS 列表 QA 有效比例”这三步再谈分析和出图。这套流程不高级但能把大部分返工挡在前期。文中这些参数、脚本和产品对照表我整理成了同一份笔记和那份 MODIS 资料放在一起需要的按文末方式自取。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/30 5:09:46 · 阅读 1.2k MODIS数据选型、下载、HDF4读取、MRT重投影与QA筛云
简介面向遥感与地理信息领域的初学者及科研应用人员这份文档围绕MODIS数据的背景、下载与处理展开系统梳理可帮助缺乏卫星遥感经验的读者快速建立从数据获取到分析应用的完整认知。包内仅含1个doc文档压缩包约503KB篇幅紧凑却覆盖了ESE与EOS计划脉络、Terra与Aqua双星互补观测机制、36个波段的光谱范围与空间分辨率、2330公里视场及705公里太阳同步极轨等关键技术指标并给出波段设计目标与典型应用领域的对照。文档同时说明NASA官网的产品与时间范围选取流程、HDF格式文件的读取方式以及借助ENVI、ArcGIS或开源工具完成预处理、辐射与几何校正、时序合成和植被指数、地表温度提取的思路。已有304人学习适合作为课程作业、课题入门或环境变化监测项目的前期参考资料。1. MODIS数据1到2天覆盖全球的36波段到底适合谁如果你要把一个区域几年的植被长势、地表温度或者云量变化跑一遍又不想在数据采购上耗预算MODIS数据大概率是最先该试的那一套。Terra 和 Aqua 两颗太阳同步极轨卫星互相补位上午星和下午星各自过境1 到 2 天能把全球表面扫一遍星上搭载的中分辨率成像光谱仪有 36 个波段光谱范围从 0.405 微米一直铺到 14.385 微米空间分辨率分 250 米、500 米、1000 米三档幅宽 2330 公里。它做不了米级精细制图也替代不了高光谱专题产品但做区域趋势、长时序对比、大范围快速摸底这套数据的性价比很难绕开。下面按产品选型、下载、HDF4 读取、MRT 投影、排错、合成验证的顺序把这份资料里能直接落地的部分拆开讲。2. 产品分级与文件名MOD、MYD、MCD三个前缀背后的选型账2.1 0级到4级先分清你要的是原始观测还是算好的产品MODIS 标准数据产品按处理深度分成 0 到 4 级这个分级不是学术分类它直接决定你下载的文件里装的是什么、后面要不要自己写算法。0 级是地面站直接收到的原始数据一般用户碰不到1A 级做了数据重建和时间配准仍然是没有辐射校正的原始计数值1B 级在 1A 基础上做了辐射校正和定标输出的是辐亮度MOD02 就是这一级。真正开始“能用”的是 2 级以后2 级产品和 1B 具有相同空间分辨率、覆盖同一地理区域比如 MOD09 表面反射率、MOD11 地表温度3 级是把多天数据放进统一时空栅格做合成比如 MOD13 的 16 天 NDVI4 级则是跑完模型之后的结果比如 MOD17 的初级生产力。选型的判断依据很简单如果你要自己做指数计算、自己控制云掩膜就从 2 级甚至 1B 起步如果你只想要一条能直接分析的时间序列3 级产品省掉的预处理时间远超它损失的那点灵活性。我一般会先确认研究问题的时间尺度——日尺度看 2 级旬月尺度看 3 级年度碳收支看 4 级别一上来就下 1B几十景原始辐亮度堆在硬盘里光定标和几何定位就能拖掉一周。2.2 44种标准产品怎么倒着选标准产品按内容分成大气、陆地、海洋三大类合计 44 种。新手最容易犯的错是照着产品编号从 MOD01 往下翻翻到哪算哪。更靠谱的做法是从应用倒推先写清楚你要的输出是什么再找对应的短名。下面这张表是日常用得最多的十几个覆盖了土地覆被、植被、温度、火情、云和气溶胶几个方向。短名内容空间分辨率时间频率典型用途MOD021B 辐亮度250/500/1000 m每轨自定义算法输入MOD03地理定位1 km每轨1B 几何校正必备MOD04气溶胶1 km10 km 产品日大气颗粒物研究MOD06云产品1 km日云量、云顶高度MOD09表面反射率250/500 m日指数计算、分类MOD11地表温度与发射率1 km日/8 天城市热岛、干旱MOD13NDVI/EVI250 m16 天植被长势MOD14热异常与火情1 km日火点监测MOD15LAI/FPAR1 km8 天生态模型输入MOD17GPP/NPP1 km8 天/年碳循环MOD35云掩膜250 m/1 km日筛云MCD12Q1土地覆被类型500 m年下垫面分类表里没有列的去产品手册里按“短名 版本号”检索不要凭记忆猜。版本号很关键同一产品不同 Collection 的处理基线不同定标系数和 QA 位定义都可能变混用会得到一条斜率对不上的时间序列。2.3 文件名拆成七段每一段都影响后面脚本MODIS 的文件名不是随手拼的拿MOD01.A2002032.1015.003.2002034184356.hdf举例七段各有含义写批处理脚本时经常要按段提取。段位示例含义脚本里的用法aMOD01产品短名MOD 为 TerraMYD 为 AquaMCD 为两者合成判断平台决定后续用哪套参数bA2002数据采集年份按年建目录c032儒略日第 32 天即 2 月 1 日转成日期和时间序列对齐d1015UTC 时分10 点 15 分判断是上午星还是下午星过境e003Collection 版本号同一批次必须统一f2002034184356生产时间年儒略日时分秒排查“同一天两个版本”时用ghdf扩展名全部为 HDF 格式决定用 HDF4 还是 HDF5 库打开最容易翻车的是 c 段和 d 段。儒略日不是月日032是 2 月 1 日不是 3 月 2 日1015是 UTC 时间换算到北京时间要加 8 小时所以一景标着A2020032.1015的数据实际对应北京时间 2 月 1 日 18 点 15 分。做日尺度统计时如果直接拿文件名里的日期当本地日期跨年那几天必错。2.4 正弦投影与h17v04瓦片编号不是随便起的大多数标准产品用正弦投影Sinusoidal赤道处按 10 度乘 10 度切网格行号从左上角 (0,0) 起算到右下角 (35,17)全球共 36 列 18 行。文件里的h17v04就是水平第 17 列、垂直第 4 行一景大致覆盖中纬度一个省级尺度的范围。瓦片编号的价值在于同一h17v04的多年数据可以直接按像元对齐做时间序列不需要先做空间重采样跨瓦片研究才需要拼接而拼接一定伴随投影和像元网格的变化。这里有个常见误解以为同一瓦片号在所有产品里范围完全一致。250 米、500 米、1 公里产品的瓦片覆盖范围是嵌套的但像元数不同MOD09GQ250 米和 MOD09GA500 米/1 公里的阵列尺寸差一倍混着读之前先看清楚SDS的维度。投影参数、瓦片号、分辨率这三个信息凑齐才能判断两景数据能不能直接逐像元运算。3. 下载从经纬度检索到wget批量拉取的四个步骤3.1 检索条件写成“经纬度框UTC时间段”数据检索页通常提供按经纬度框选、按瓦片号选、按日期范围选三种方式。区域研究用经纬度框最直接输入左上角和右下角坐标再填时间范围。时间一定按 UTC 填要取北京时间 1 月 16 日全天的数据UTC 范围是 1 月 15 日 16:00 到 1 月 16 日 16:00直接填 1 月 16 日 00:00 到 24:00 会漏掉前 8 小时。检索结果里会列出每景的文件名、大小、平台、版本先别急着全选用文件名里的 d 段确认过境时间是不是你要的上午星或下午星。检索结果一般支持导出成文本清单这个清单就是后面批量下载的输入。清单里通常带完整下载链接但链接可能带会话参数过几小时就失效所以脚本里更稳的做法是只保留文件名把数据目录作为变量拼前缀。L1B 数据记得把 MOD03 一起勾上它是 1 公里、500 米、250 米产品做几何纠正的定位文件漏下这一景后面打开 MOD02 会发现没有经纬度返工重下最费时间。3.2 账号与token把登录这一步提前做掉数据平台现在基本都走统一登录匿名 FTP 的老路子已经不好用。注册账号后在个人中心生成一个下载 token之后所有请求带Authorization: Bearer头即可不用每次输密码也不会因为会话过期在半夜断掉。token 等同于密码别写进公开仓库脚本里用环境变量传入。下面这段是批量下载的骨架把清单文件、数据目录、token 三个变量换掉就能跑。#!/usr/bin/env bash # 批量下载 MODIS 数据断点续传、失败重试、限速 set -u export ED_TOKEN你的下载token BASE_URLhttps://检索页复制下来的数据目录 OUT_DIR./modis_raw mkdir -p $OUT_DIR while read -r name; do [ -z $name ] continue # 跳过空行 [ -f $OUT_DIR/$name ] continue # 已存在则跳过 wget -c -nc --tries5 --timeout60 --wait2 \ --headerAuthorization: Bearer ${ED_TOKEN} \ -P $OUT_DIR ${BASE_URL}/${name} || echo FAILED $name failed.txt done urls.txt-c是断点续传网络中断后重跑不会从头拉-nc防止覆盖已完成的文件--tries5和--timeout60对付偶发的连接超时--wait2每次请求间隔两秒批量几百景时能明显降低被限流的概率。失败的文件名单独记到failed.txt第二天换时段重跑比盯着终端一屏屏翻要省心。数据目录不要手拼直接从检索结果里复制一层不同产品、不同版本的目录结构不完全一样。3.3 下载完先校验别等读的时候才发现是HTML下载中断最隐蔽的结果不是文件缺失而是文件存在但只有几 MB打开一看是登录失效后返回的 HTML 错误页。HDF4 文件有固定的文件头正常文件开头是 HDF 标识用几行 Python 就能批量筛一遍顺便统计每个文件的 SDS 数量数量为零的直接标红。from pathlib import Path from pyhdf.SD import SD, SDC bad [] for p in Path(modis_raw).glob(*.hdf): try: h SD(str(p), SDC.READ) n len(h.datasets()) # SDS 数量 h.end() if n 0: bad.append((p.name, no SDS)) except Exception as e: bad.append((p.name, repr(e))) for name, reason in bad: print(name, reason) print(待重下:, len(bad))SD打开失败基本就是文件截断或根本不是 HDF能打开但datasets()为空多半是下到了错误页。这个检查跑一遍几十秒比后面在重投影环节报一个看不懂的错要划算。顺手把文件按“产品短名/年/儒略日”建目录后面 MRT 批处理和质量筛选都按这个层级遍历省得再写一套路径解析。3.4 清单、日志、版本三样东西留好批量下载最怕两周后回头找“当时到底下了哪几景”。我的习惯是每批数据在目录下留三个文件原始清单urls.txt、下载日志failed.txt、以及一个manifest.csv记录文件名、大小、下载时间、Collection 版本。看着繁琐等审稿人或合作者问“你这批数据版本是不是统一”的时候这份清单就是后悔药。尤其是 MOD 和 MYD 混用、或者跨 Collection 拼接时版本不一致带来的反射率偏差不会报错只会让趋势线的斜率悄悄偏掉。4. HDF4结构与MRT重投影SDS读取到Albers输出的完整链路4.1 HDF4的层级SDS、Vdata、Vgroup与属性MODIS 的.hdf是 HDF4 格式不是 HDF5两者不兼容。HDF4 内部是层次结构最常用的是 SDSScientific Dataset也就是一维到多维的数组反射率、温度、QA 都是 SDSVdata 存表格型数据Vgroup 负责把相关数据集组织在一起每个 SDS 自带属性定标系数、填充值、有效范围都藏在属性里。用pyhdf.SD打开文件后先列一遍数据集名和维度能避免“凭记忆 select 一个不存在的名字”这种低级错误。from pyhdf.SD import SD, SDC import numpy as np f MOD09GA.A2020032.h17v04.006.2020040034512.hdf hdf SD(f, SDC.READ) for name, info in hdf.datasets().items(): print(name, info[0], info[1:]) # 名称、类型、维度输出里会看到sur_refl_b01到sur_refl_b07、state_1km、solar_zenith、satellite_zenith等。info[0]是数据类型info[1:]是各维长度。MOD09GA 的 500 米反射率维度大约是 1200×12001 公里 QA 是 600×600两者不能直接逐像元运算要么先重采样要么用 500 米波段自己的 QA。读之前先看清楚维度比读完之后发现形状对不上再回头查要快得多。4.2 读反射率scale_factor和_FillValue必须一起处理MODIS 的反射率 SDS 存的是整型不是物理值直接读出来是几千上万的整数。每个 SDS 的属性里有scale_factor、add_offset、_FillValue、valid_range正确顺序是先按填充值掩膜再乘比例加偏移。顺序反了填充值也会被定标成一个看似合理的负值后面统计最小值时就会莫名其妙出现 -0.1 这种数。sds hdf.select(sur_refl_b01) raw sds.get().astype(float32) attr sds.attributes() scale attr.get(scale_factor, 1.0) offset attr.get(add_offset, 0.0) fill attr.get(_FillValue, -28672) refl raw * scale offset # 定标 refl[raw fill] np.nan # 先掩膜再统计 print(refl.shape, np.nanmin(refl), np.nanmax(refl)) sds.end()MOD09 的scale_factor通常是 0.0001_FillValue是 -28672valid_range在 -100 到 16000 之间。算完的反射率理论值域在 0 到 1 附近出现大于 1 的值先别慌看看是不是云边缘或观测几何极端的情况再检查有没有把valid_range之外的像元算进来。对于 MOD11 的地表温度scale_factor是 0.02填充值是 0单位是开尔文要转摄氏度还得减 273.15这些都在产品手册的属性表里别靠猜。4.3 geo_field与经纬度查找表几何定位的两种来源1B 数据本身不带完整的逐像元经纬度MOD03 才是定位文件里面有Latitude、Longitude、SolarZenith、SensorZenith等数据集维度与 1 公里波段一致。做几何校正时可以把 MOD03 的经纬度读出来按行列号对应到 MOD02 的像元上也可以在 MRT 重投影时直接把 MOD03 作为定位输入。两种方式的区别在于前者适合自己写校正算法、做精确的像元级配准后者适合批量出图MRT 内部会处理定位速度快但可控参数少。读 MOD03 的经纬度时注意维度顺序HDF4 的数组是行优先Python 读出来是[row, col]和影像的上下左右方向一致但有些产品会做翻转或偏移拿一景有明显地物特征的数据比如海岸线叠一下就能确认。经纬度是浮点数组文件体积不小批量处理时不要一次性把几十景的经纬度都读进内存按瓦片循环、读完即释放。4.4 MRT先拼接再投影Albers参数一次配好MRTMODIS Reprojection Tool负责重投影、拼接、重采样、裁剪和波段子集。它的处理顺序是先把同一投影下的多幅影像拼接再做投影转换这个顺序不能反如果先投影再拼接接缝处的像元网格对不齐会出现一条明显的错位带。参数文件.prm里最关键的几项是输出投影类型、标准纬线、中央经线、基准面和重采样方式。国内区域制图常用 Albers 等积圆锥投影双标准纬线 25 度和 47 度、中央经线 105 度是一套通用参数省级范围可以按轮廓再调比如东西跨度大的区域把中央经线往中间挪。参数常用值说明OUTPUT_PROJECTION_TYPEALBERS等积圆锥适合面积统计STDPR125.0第一标准纬线STDPR247.0第二标准纬线CENTRAL_MERIDIAN105.0中央经线按研究区调整ORIGIN_LATITUDE0.0原点纬度DATUMWGS84基准面RESAMPLING_TYPEBILINEAR连续量用双线性分类产品用最近邻OUTPUT_PIXEL_SIZE500输出像元大小单位米重采样方式要按数据类型选反射率、温度这类连续量用双线性土地覆被、QA 这类分类或位标志必须用最近邻否则一个 3 类像元被插成 2.5 类后面统计类别面积全是错的。输出像元大小不要小于原始分辨率500 米产品输出 250 米不会增加信息只会让文件大四倍。参数文件写好后用-R传入文件清单、-O指定输出目录批量跑。MRT_HOME/opt/mrt/bin $MRT_HOME/resample -P albers.prm -R filelist.txt -O ./out_albers -V-V打开详细日志每景处理成功还是失败都会打印批量几百景时把日志重定向到文件失败景单独重跑。MRT 对文件名有要求输入清单里写绝对路径更稳输出文件名会自动加波段后缀后期按规则批量改名即可。跑完先抽一景和原始正弦投影数据叠一下看边界是不是吻合再接下去做时间序列。5. 避坑排查下载、读取、投影五个翻车现场5.1 日期差一天UTC和儒略日叠在一起现象按“1 月 16 日”检索下回来的数据文件名是A2003015做日统计时这一天的值对不上。原因文件名里的日期是数据采集的 UTC 儒略日检索页填的时间范围也是 UTC。北京时间 1 月 16 日 00:00 到 08:00 属于 UTC 1 月 15 日这 8 小时的数据会落在A2003015里。如果只按文件名日期分组就会把本地日期的一整天拆成两半。解决统一在脚本里做时区转换把文件名解析成 UTC 时间后再加 8 小时得到本地日期按本地日期分组或者干脆全程用 UTC 做时间轴出图时再标注本地时间。两种都行关键是别混。5.2 只下MOD02打开没有地理定位现象MOD02 能打开但找不到经纬度数据集做裁剪和重投影时无从下手。原因MOD02 是 1B 辐亮度几何定位信息在配套的 MOD03 文件里两者是分开下载的。检索时只勾了 MOD02MOD03 就漏了。解决检索清单里同时勾选同一天、同一轨的 MOD03下载后按文件名里的时间戳配对。MRT 重投影时把 MOD03 作为定位输入或者在 Python 里读 MOD03 的Latitude、Longitude做自定义校正。配对规则写进脚本别靠人工比对。5.3 用h5py读HDF4报错还找不到原因现象h5py.File(MOD09GA...hdf)直接抛异常提示不是 HDF5 文件。原因MODIS 标准产品是 HDF4h5py只认 HDF5。两个格式同名不同代API 完全不兼容。网上不少示例代码混着写复制过来就报错。解决HDF4 用pyhdfHDF5 用h5py先看文件头再选库。GDAL 两类都能读gdalinfo一跑就知道格式。批量处理前统一检查一遍文件头避免一半文件用一种库、另一半用另一种库。5.4 反射率值域不对scale_factor没乘现象读出来的反射率在 0 到 10000 之间做指数计算全是异常值。原因SDS 存的是整型物理值要乘scale_factor再加add_offset。只读了get()没读属性或者顺序反了先定标再掩膜填充值被算成了有效值。解决固定写成“读原始数组 → 读属性 → 掩膜填充值 → 定标”四步封装成函数复用。每次换产品先打印attributes()看一遍MOD09、MOD11、MOD13 的 scale 和 fill 都不一样不能一套代码硬套所有产品。5.5 MRT输出接缝、错位与内存爆掉现象多景拼接后中间有一条错位带或者跑到一半进程被系统杀掉。原因错位多半是先投影再拼接或者输入清单里混了不同 Collection、不同分辨率的文件内存爆掉通常是一次性把几十景全丢给 MRT或者输出像元设得过细。解决确认 MRT 的处理顺序是先拼接后投影输入清单按产品、版本、分辨率分组不同组分开跑批量时按瓦片或按年分批每批 10 到 20 景跑完检查接缝再继续。输出像元不小于原始分辨率连续量用双线性、分类量用最近邻。跑之前先用两景做小样本测试参数确认无误再放开批量。6. 从单日瓦片到月度合成QA筛云与结果验证的一个笨办法单日反射率直接平均是最容易出问题的做法云、云影、观测几何差异都会混进月度值里。我的习惯是先读state_1km做云筛再按像元统计有效观测数有效数低于阈值比如月内少于 5 次的像元直接置空不参与后续分析。MOD09GA 的state_1km是位标志低两位表示云状态0 为晴空1 为有云2 为混合3 为未定义先掩膜再定标。import numpy as np from pyhdf.SD import SD, SDC def read_clear_refl(path, bandsur_refl_b01): hdf SD(path, SDC.READ) raw hdf.select(band).get().astype(float32) attr hdf.select(band).attributes() qa hdf.select(state_1km).get() hdf.end() refl raw * attr.get(scale_factor, 1.0) attr.get(add_offset, 0.0) refl[raw attr.get(_FillValue, -28672)] np.nan cloud qa 0b11 # 低两位云状态 refl[(cloud ! 0) | (qa 1)] np.nan # 非晴空和填充值都剔除 return refl这段代码只做了一件事把非晴空像元变成NaN后续用np.nanmean做月度合成时自然跳过。跑完一个月后别急着出图先看两个数有效像元比例和月度均值的直方图。有效比例低于 30% 的区域说明这个月云太多合成结果的可信度要打折扣直方图出现双峰或者异常拖尾多半是 QA 位定义用错了版本或者把云影漏掉了。拿一景已知水体或城市核心区的像元值做对照水体反射率在可见光波段应该低而平城市建成区在短波红外明显偏高对不上就回去查定标和掩膜。从那以后我每次下完一批数据都强制先跑一遍“文件头检查 SDS 列表 QA 有效比例”这三步再谈分析和出图。这套流程不高级但能把大部分返工挡在前期。文中这些参数、脚本和产品对照表我整理成了同一份笔记和那份 MODIS 资料放在一起需要的按文末方式自取。希望帮到你。本文还有配套的精品资源点击获取