简介一份围绕MODIS数据入门与应用的文档面向遥感、地理信息、环境监测等方向的学生和初级从业者帮助梳理从传感器背景到数据获取与处理的关键知识。压缩包内仅含1个doc文件约503KB体量轻便适合快速浏览和按主题查阅。内容涵盖NASA地球科学事业与EOS计划、Terra/Aqua卫星分工、MODIS 36个波段及空间分辨率等技术指标并延伸到产品下载流程、HDF格式读取以及ENVI、ArcGIS等工具的处理思路。文中还涉及预处理、校正、合成与分析等环节可辅助理解植被指数、地表温度等常见应用。已有304人学习适合用作MODIS入门参考和课堂/项目前的基础资料。1. 600MB 的 HDF 拖进软件只显示一片灰MODIS数据到底该从哪儿下手第一次接触 MODIS数据的人几乎都会卡在同一个画面上文件下下来了双击打不开拖进 GIS 软件只弹出一个灰蒙蒙的空框图层列表里全是EV_1KM_RefSB、sur_refl_b01、250m_16_days_NDVI这类没头没尾的名字。这不是文件坏了而是 MODIS 的组织方式和普通影像完全不同——一个.hdf里塞着几十个科学数据集、属性表和定标参数波段被拆成一个个“抽屉”你不点名去取它就什么都不给你看。这份资料覆盖的范围比一个 HDF 文件大得多从 Terra/Aqua 双星的轨道设计、36 个波段的物理含义到 L0~L4 五级产品体系、44 种标准产品的命名规则再到检索下单、批量下载、MRT 重投影与 Albers 参数配置最后落到定标和有效值处理。适合第一次做植被指数、地表温度或土地覆被时间序列的人也适合被 HDF4 折磨过一轮、想把下载与批处理固化成脚本的老手。往下读之前先记住一句话MODIS 的坑八成不在算法而在“你下错产品”和“你没做定标”。2. 36 个波段与 44 种产品先定产品谱系再谈下载绝大多数翻车都不是从代码开始的是从检索页开始的。选了MOD021KM想做 NDVI或者下了MOD13Q1想反演气溶胶后面所有工作都是白费。这一章要做的事只有一件在看到下载按钮之前先能说清自己要的是哪一级、哪个短名、哪个分辨率桶。2.1 三个分辨率桶与两类波段为什么不能“全波段一把梭”MODIS 的 36 个波段不是平均分布的它被硬切成了三个分辨率桶250m 只有波段 1、2500m 是波段 3~7剩下 8~36 全部是 1km。这个设计直接决定了很多产品为什么长成现在这样——MOD13Q1的 NDVI 是 250m因为它只依赖波段 1、2MOD11A2的地表温度是 1km因为热红外波段本来就只有 1km。波段还分两类1~19 是反射光谱段探测的是地表和云对太阳短波辐射的反射物理量是反射率20~36 是发射光谱段探测的是地物自身的热辐射物理量是亮温。这两类数据的定标方式、有效值范围、甚至可视化方式都不一样。把波段 31 当反射率去算归一化指数结果会是一张全黑或者全白的图很多人到这一步还以为是自己算错了。常用波段的用途可以对着下面这张表看不用背全 36 个记住七个就够日常用波段光谱范围(nm)星下点分辨率典型用途1620–670250m叶绿素吸收NDVI 红波段2841–876250m近红外植被与云检测3459–479500m土壤与植被差异4545–565500m绿峰水体浑浊度51230–1250500m叶片与冠层水分差异61628–1652500m雪与云判别72105–2155500m土地覆盖与云特征20–233.66–4.08 μm1km海面温度、火点与火山298.40–8.70 μm1km表面温度3110.78–11.28 μm1km地表温度主波段3211.77–12.27 μm1km分裂窗云顶高度表里能读出一个规律凡是做植被、土地覆盖、积雪的基本都落在 1~7凡是做温度、火点、云的基本都落在 20 以后。这也解释了为什么陆地产品多是 250m/500m大气和海洋产品多是 1km。提示波段 13 有 13h 和 13l 两个高增益与低增益波段 14 同理。做叶绿素荧光时选错增益数值会差一个量级且不会有任何报错。2.2 L0 到 L444 种标准产品里真正会用到的十来个MODIS 标准产品按处理深度分五级。L0 是地面站直接收到的原始比特流L1A 是重建、时间配准、加了辅助数据之后的产物L1B 是做过辐射定标的标准产品也就是常说的MOD02L2 以上是在 L1B 基础上做应用处理得到的。这里有个容易被忽略的关键差别L1B 和 L2 是沿轨条带swath每一景跟着卫星走是斜的L3 和 L4 是投影到正弦格网上的瓦片tile是正的、可按瓦片号索引的。这个差别决定了要不要做几何校正。拿 L1B 直接出图你会得到一条歪着的、带弯曲边缘的影像必须先配合MOD03地理定位文件做几何纠正而拿MOD13Q1这类 L3 瓦片它已经在正弦投影下规整好了直接重投影就行。新手最常走错的路径就是想做时间序列却从 L1B 开始结果光是几何校正和拼接就耗掉一周。下面这张表是实际项目里出现频率最高的十来个短名覆盖了陆地、大气、海洋三个方向短名内容分辨率时间粒度平台前缀MOD021KML1B 辐射亮度250m/500m/1km5 分钟段TerraMOD03地理定位文件1km同 L1BTerraMOD04_L2气溶胶10km/1km日TerraMOD06_L2云产品1km日TerraMOD09GA地表反射率500m/1km日TerraMOD11A2地表温度与发射率1km8 天合成TerraMOD13Q1NDVI/EVI250m16 天合成TerraMOD15A2HLAI/FPAR500m8 天合成TerraMOD17A3H净初级生产力500m年TerraMCD12Q1土地覆盖类型500m年TerraAquaMCD15A2LAI/FPAR 融合1km8 天TerraAquaMOD10A2积雪覆盖500m8 天Terra前缀是本篇最该刻进脑子的一条规则MOD是 Terra上午星MYD是 Aqua下午星MCD是两颗星融合。同一个产品名换成MYD前缀过境时间就从上午 10:30 变成下午 13:30。做物候或者能量平衡时上午星和下午星的地表温度能差好几度混着用会得到一条根本解释不通的曲线。选型逻辑可以归纳成一句话做算法反演、需要原始辐射亮度或自定义指数走 L1B MOD03做区域制图、时间序列、统计分析走 L3/L4 瓦片产品。前者自由度高但工作量大后者开箱即用但只有官方定好的那几十个变量。2.3 文件名逐段拆解从 MOD13Q1.A2020032.h25v04.006 读出全部元数据MODIS 的文件名本身就是一份元数据只是被压缩成了一串点分字段。以MOD13Q1.A2020032.h25v04.006.2020052184512.hdf为例MOD13Q1产品短名MOD 是 Terra13 表示植被指数Q 表示 250m 分辨率档1 表示产品系列内的编号。A2020032A加四位年份加三位年积日。2020 年第 32 天也就是 2 月 1 日。h25v04正弦投影下的瓦片行列号水平方向 h00~h35垂直方向 v00~v17每块在赤道处约 10°×10°。006Collection 版本号。不同版本之间的定标系数和算法可能改过跨版本拼时间序列会引入系统性偏差。2020052184512生产时间格式是yyyydddhhmmss即 2020 年第 052 天 18 时 45 分 12 秒。.hdf文件格式。L1B 的文件名结构略有不同会多出一段采集时刻比如MOD021KM.A2020032.0345.061.2020032213456.hdf中间的0345就是 UTC 03:45 的过境时刻这个字段在做多星匹配时非常有用。手工数位太容易出错我一般直接写个函数把文件名解析成结构化字典后面做批量索引、按日期筛选、按版本分组都靠它import re from datetime import datetime, timedelta # 兼容 L1B带时刻、无瓦片和 L3带瓦片、无时刻两种命名 PATTERN re.compile( r^(?Pproduct[A-Z]{3}\d{2}[A-Z0-9]{0,3}) r\.A(?Pyear\d{4})(?Pdoy\d{3}) r(?:\.(?Putc_hhmm\d{4}))? r(?:\.(?Ptileh\d{2}v\d{2}))? r\.(?Pversion\d{3}) r\.(?Pprod_time\d{13}) r\.(?Pexthdf)$ ) def parse_modis_name(name): m PATTERN.match(name) if not m: return None d m.groupdict() # 年积日转真实日期1 月 1 日是第 1 天所以要减 1 d[acq_date] datetime(int(d[year]), 1, 1) timedelta(daysint(d[doy]) - 1) # 生产时间 yyyydddhhmmss p d[prod_time] d[prod_datetime] ( datetime(int(p[0:4]), 1, 1) timedelta(daysint(p[4:7]) - 1) ).replace(hourint(p[7:9]), minuteint(p[9:11]), secondint(p[11:13])) return d print(parse_modis_name(MOD13Q1.A2020032.h25v04.006.2020052184512.hdf)) print(parse_modis_name(MOD021KM.A2020032.0345.061.2020032213456.hdf))这段代码的关键在两个地方。一是正则里的两组可选字段都用(?:...)?包起来因为 L1B 有时间无瓦片、L3 有瓦片无时间如果不加可选两类文件必然有一类匹配不上。二是年积日的换算要减 1datetime(year, 1, 1) timedelta(daysn-1)才等于第 n 天写成daysn会整体偏移一天这种错在按日聚合时特别隐蔽——数据一天不缺但每一天都对不齐。有了这个解析结果按年份、按瓦片、按版本分组或者筛出某段日期内的文件都只是普通的列表过滤不用再靠文件名切片数位置。2.4 用 gdalinfo 验货下载完先看 SDS 清单再写脚本文件到手之后不要急着写批处理先花两分钟验一下货。GDAL 对 HDF4 的支持是开箱即用的一条gdalinfo就能把文件里所有子数据集的完整路径打印出来# 列出 HDF4 文件内所有子数据集SUBDATASET gdalinfo MOD13Q1.A2020032.h25v04.006.2020052184512.hdf | grep -A2 SUBDATASET | head -40 # 只看某个具体子数据集的元信息尺寸、投影、类型、NODATA gdalinfo HDF4_EOS:EOS_GRID:MOD13Q1.A2020032.h25v04.006.2020052184512.hdf:MODIS_Grid_16DAY_250m_500m:250m_16_days_NDVI输出里最该关注四项Size is给出行列数正弦瓦片通常是 4800×4800 这个量级Coordinate System is给的是正弦投影带一个中央经线参数Type一般是Int16说明它存的是缩放后的整数而不是浮点物理量NoData Value决定后续掩膜怎么写。子数据集路径是有固定语法的HDF4_EOS:EOS_GRID:文件路径:网格名:数据集名。网格名和数据集名在不同产品、不同版本里会变所以脚本里不要硬编码正确的做法是先用GetSubDatasets()把清单取出来再用关键字过滤from osgeo import gdal fp MOD13Q1.A2020032.h25v04.006.2020052184512.hdf ds gdal.Open(fp) for path, desc in ds.GetSubDatasets(): print(desc)desc一般是[4800x4800] 250m_16_days_NDVI MODIS_Grid_16DAY_250m_500m这种格式用if NDVI in desc就能精准定位比手写完整路径稳得多。注意同一个产品在不同 Collection 版本下网格名可能不同。脚本里写死MODIS_Grid_16DAY_250m_500m在换版本后会静默失败返回 None而不是报错。所有硬编码路径前都建议先跑一遍清单打印。3. 从检索页到硬盘1B 与 L3 产品的两条下单路线和批量脚本搞清楚要什么之后下一个问题是去哪里拿。MODIS 的下载入口不是统一的而是按产品类型分在两个体系里。走错门不是拿不到数据而是根本搜不出结果很多人到这一步会怀疑是自己网络或者账号有问题。3.1 两类门户的分工swath 产品走订单tile 产品走检索L1B 和 L2 级产品MOD02、MOD03、MOD04、MOD06、MOD35等沿轨条带组织检索维度是时间加空间范围走的是大气产品那条线。这类门户的检索界面通常要求你填起止时间、经纬度框然后列出覆盖这个范围的若干条轨道段勾选后提交订单。L3/L4 产品MOD09、MOD11、MOD13、MOD15、MOD17、MCD12Q1等按正弦瓦片组织检索维度是时间加瓦片号。这类门户一般提供地图点选或者经纬度反查瓦片号的功能选中瓦片和时间范围后直接给出一批文件链接不需要走订单流程。两条路线的差别不只是界面。swath 产品的单景覆盖范围是一条窄带一天之内同一地区可能被多条轨道部分覆盖要做完整覆盖就得把相邻轨道段都下下来再拼接瓦片产品是规整的 10°×10° 网格同一瓦片同一天只有一景取数逻辑简单得多。这也是为什么做区域时间序列的人都优先选瓦片产品。实际检索时还有一个绕不开的问题时区。门户上显示和填写的时间全是 UTC北京时间比 UTC 快 8 小时。要拿到北京时间 1 月 16 日全天的数据UTC 检索范围应该是 1 月 15 日 16:00 到 1 月 16 日 16:00而不是 1 月 16 日 00:00 到 24:00。按后者填会丢掉北京时间当天前 8 小时的数据而且丢得毫无痕迹——文件数看着正常只是早上那几景没了。提示做长时间序列时我习惯把检索范围前后各多放一天下载完再按解析出的采集时刻裁掉多余的。宁可多下几景也不要在边界上少一天。3.2 空间范围、瓦片与覆盖完整性经纬度框看着简单坑却不少。第一个坑是跨 180° 经线如果研究区横跨日界线用minLon170, maxLon-170这种写法大部分检索接口会直接返回空结果。正确做法是拆成两个框分别检索或者干脆改用瓦片号检索。第二个坑是研究区跨瓦片。一个 10°×10° 的瓦片对中等尺度的研究区来说往往不够比如经度 111°–120°、纬度 35°–43° 这样的范围至少要覆盖 2×2 共四块瓦片边缘还可能再搭上一块。少下一块的表现是拼接结果边缘缺一条而缺的那条往往正好是研究区边界。第三个坑是瓦片号手算。正弦投影下瓦片的经度跨度随纬度变化赤道处是 10°高纬度会变宽所以用简单的等经纬度除法算出来的瓦片号经常差一格。可靠的做法是在检索页的地图上直接框选让系统把覆盖的瓦片列出来把这几个瓦片号记下来后面所有检索复用同一组。3.3 批量下载脚本令牌认证与断点续传门户上的“一键下载”适合拿几景试试水真正做项目必须走脚本。现在的数据门户普遍用 Earthdata 账号登录脚本里用账号密码硬编码既不安全也不稳定标准做法是在账号设置里生成一个令牌token通过 HTTP 头或者 cookie 传给下载工具。#!/usr/bin/env bash set -euo pipefail # 数据门户的归档根路径具体值以检索页给出的下载链接前缀为准 BASEhttps://data-portal-host/archive/allData PRODUCTMOD13Q1 VERSION006 YEAR2020 DOY032 # 令牌放在环境变量里不要写进脚本提交到仓库 TOKEN${EARTHDATA_TOKEN:?请先导出 EARTHDATA_TOKEN} # 断点续传 失败重试 限速避免被服务端掐断 wget --headerAuthorization: Bearer ${TOKEN} \ --continue \ --tries5 \ --waitretry10 \ --timeout60 \ --limit-rate8m \ --no-check-certificate \ -P ./raw \ ${BASE}/${PRODUCT}/${VERSION}/${YEAR}/${DOY}/ # 只保留 hdf 主文件清掉同目录下的 .xml 元数据如果不需要溯源 find ./raw -name *.hdf.xml -delete echo 已下载 $(ls -1 ./raw/*.hdf | wc -l) 个 HDF 文件脚本里有几个参数值得单独说。--continue是断点续传的开关MODIS 单景动辄几百兆一次下完的概率并不高--tries5 --waitretry10让它在遇到 5xx 时退避重试而不是直接放弃--limit-rate8m是血泪经验并发拉满容易被服务端限流反而更慢压到 8MB/s 左右反而稳定。如果要下的是按订单生成的 swath 产品脚本思路一样只是 URL 列表要先从订单页导出成文本文件再用wget -i list.txt批量拉。# 从订单页导出的链接清单逐行下载每行一个 URL while IFS read -r url; do [ -z $url ] continue wget --headerAuthorization: Bearer ${EARTHDATA_TOKEN} \ --continue --tries5 -P ./raw $url done order_links.txtwhile read这种写法看起来笨但比xargs -P的好处是失败可感知某一行断了日志里能定位到具体是哪个文件重跑时--continue会接着下。并行下载在文件数上千时确实快但一旦出现半截文件排查成本远高于省下的那点时间。3.4 下载完先做完整性检查脚本跑完不算结束还要验一遍。三个检查动作文件数量对不对、文件大小有没有异常的、GDAL 能不能打开。# 1) 逐个文件做 HDF4 结构校验打不开的会报错 for f in ./raw/*.hdf; do gdalinfo $f /dev/null 21 || echo 损坏: $f done # 2) 找出明显偏小的文件半截下载 find ./raw -name *.hdf -size -1M -print # 3) 统计文件数和检索页给出的数量对一下 ls -1 ./raw/*.hdf | wc -l半截文件是批量下载里最阴的一种它存在、它有大小、文件名也对但gdalinfo打不开或者子数据集缺一半。如果不在入库前筛出来等跑到批处理第 300 个文件时再崩重跑成本会翻好几倍。我的习惯是下载完立刻跑一遍上面的循环把坏文件挪到./bad目录重新下一遍。4. HDF4 拆箱子数据集、定标系数与 MRT 重投影文件完好、路径也对接下来才是真正的处理。这一步要解决三件事怎么把数据从 HDF 里取出来、怎么把整数 DN 变成有物理意义的反射率或亮温、怎么把正弦投影转成能和其他数据叠加的投影。4.1 HDF4 与 HDF-EOS一个文件里的三层结构MODIS 用的是 HDF4 格式部分新产品线已转向 HDF5 或 NetCDF但经典产品仍是 HDF4。HDF4 本身是个通用容器NASA 在它上面加了一层约定叫 HDF-EOS用来规范地理数据的组织方式。所以一个 MODIS 文件里其实有三层东西第一层是科学数据集SDS就是真正的栅格数据每个波段、每个指数都是独立的 SDS有自己的尺寸、类型、填充值。第二层是属性attributes分全局属性和 SDS 属性定标系数scale_factor、add_offset、有效范围valid_range、填充值_FillValue都藏在这里不读属性直接拿数组算结果一定是错的。第三层是 Vdata 和 Vgroup存储地理定位查找表、质量标志的位定义等辅助信息元数据浏览工具能看到普通处理一般碰不到。这解释了为什么不能把.hdf当.tif用。TIFF 一个文件基本对应一张图HDF 一个文件是一个小型数据库你必须先“查表”再“取数”。4.2 三条读取路径GDAL、pyhdf、h5py 的适用边界读取 MODIS 的 HDF4常见的有三条路各有各的适用场景选错了会浪费很多时间。第一条是 GDAL也是我最推荐的默认路径。它把 HDF4 的子数据集抽象成了虚拟文件路径读出来直接就是带投影和地理变换的栅格后面的重投影、裁剪、统计全都能用 GDAL 生态一条龙做掉。缺点是必须先知道子数据集全名靠GetSubDatasets()枚举。from osgeo import gdal import numpy as np fp MOD13Q1.A2020032.h25v04.006.2020052184512.hdf ds gdal.Open(fp) # 从清单里按关键字找目标子数据集避免硬编码网格名 target [p for p, d in ds.GetSubDatasets() if 250m_16_days_NDVI in d][0] band gdal.Open(target).GetRasterBand(1) arr band.ReadAsArray().astype(float32) scale band.GetScale() or 1.0 # 常见为 1e-4 offset band.GetOffset() or 0.0 nodata band.GetNoDataValue() # 先掩膜再定标顺序反了会把填充值也乘上缩放系数 arr[arr nodata] np.nan ndvi arr * scale offset print(范围:, np.nanmin(ndvi), np.nanmax(ndvi)) # 正常应落在 -0.2 ~ 1.0这段代码有两个细节。一是band.GetScale()在某些版本或某些子数据集上返回None所以要给默认值 1.0否则整个数组会变成 NaN。二是掩膜必须在定标之前做如果先乘 1e-4 再比较填充值-3000已经变成了-0.3和nodata对不上等于没掩。第二条是 pyhdfHDF4 的原生 Python 绑定。它的优势是能拿到 GDAL 抽象掉的东西比如 Vdata 里的质量标志位定义、完整的属性字典。做质量控制、精度验证、或者要读一些冷门产品时绕不开它。from pyhdf.SD import SD, SDC hdf SD(MOD11A2.A2020032.h25v04.006.2020052184512.hdf, SDC.READ) print(hdf.datasets().keys()) # 所有 SDS 名称 sds hdf.select(LST_Day_1km) data sds[:] # 读出为 numpy 数组 attrs sds.attributes() # 定标系数、有效范围都在这里 print(attrs.get(scale_factor), attrs.get(valid_range))第三条是 h5py这条路走不通必须说清楚。h5py 只能读 HDF5遇到 HDF4 会直接抛OSError: Unable to open file (file signature not found)。网上不少教程把两者混为一谈新手照着敲完报错还以为是文件损坏。判断方法很简单gdalinfo能打开且显示HDF4_EOS就是 HDF4用 GDAL 或 pyhdf显示HDF5才是 h5py 的战场。4.3 MRT 命令行先拼接后重投影的固定动作MODIS 的正弦投影和常见的等经纬度、Albers、UTM 都对不上必须重投影才能和矢量边界、其他遥感数据叠加。官方工具 MRTMODIS Reprojection Tool是这条链路上最稳的一个虽然有图形界面但批量处理必须用命令行。MRT 的核心逻辑是“先拼接、后投影”。多幅同一天、相邻瓦片的影像mrtmosaic会先按原正弦投影拼成一整块resample再做投影转换。这个顺序不能反反过来先投影再拼接接缝处会因为重采样误差出现明显的条带。# 第一步把同一天的四块瓦片列进清单做拼接 cat tile_list.txt EOF MOD13Q1.A2020032.h25v04.006.2020052184512.hdf MOD13Q1.A2020032.h26v04.006.2020052184512.hdf MOD13Q1.A2020032.h25v05.006.2020052184512.hdf MOD13Q1.A2020032.h26v05.006.2020052184512.hdf EOF # 光谱子集按需选这里保留 3 个 SDS顺序与文件内 SDS 顺序一致 mrtmosaic -i tile_list.txt -o mosaic_2020032.hdf -s 1 1 1 # 第二步用参数文件做重投影 resample -p albers.prm -i mosaic_2020032.hdf -o mosaic_2020032_albers.tif-s 1 1 1这个参数容易被忽略它的意思是“每个输入文件的第 1 个、第 1 个、第 1 个 SDS 参与拼接”位置对应输出波段。如果拼的是MOD13Q1它的 SDS 顺序里 NDVI 是第 1 个、EVI 是第 2 个想要两个都保留就得写成-s 1 1 2 2四个文件各取前两个。写错的后果不是报错而是输出波段数和你预期的不一样。albers.prm是纯文本的投影参数文件内容如下INPUT_FILENAME mosaic_2020032.hdf SPECTRAL_SUBSET ( 1 1 1 1 ) SPATIAL_SUBSET_TYPE INPUT_LAT_LONG SPATIAL_SUBSET_UL_CORNER ( 43.0 111.0 ) SPATIAL_SUBSET_LR_CORNER ( 35.0 121.0 ) OUTPUT_FILENAME mosaic_2020032_albers.tif RESAMPLING_TYPE BILINEAR OUTPUT_PROJECTION_TYPE ALBERS OUTPUT_PROJECTION_PARAMETERS ( 0.0 0.0 25.0 47.0 105.0 0.0 0.0 0.0 ) DATUM WGS84 OUTPUT_PIXEL_SIZE 250参数串的顺序是(SMajor SMinor STDPR1 STDPR2 CenMer OriginLat FE FN)。前两个填 0.0 时MRT 会自动套用 WGS84 的长短半轴6378137.0 和 6356752.314245 米不用手打那一长串。STDPR125、STDPR247、CenMer105是国内中低纬度区域制图最常用的一组双标准纬线选在 25° 和 47°中央经线 105°这样全国范围内的面积变形能控制得比较均匀。如果是单个省份标准纬线可以按该省南北跨度收窄中央经线取省域中心经度变形会更小。RESAMPLING_TYPE选BILINEAR是连续量NDVI、反射率、温度的常规做法如果是土地覆盖这类分类产品必须换成NEAREST_NEIGHBOR双线性会把类别码插值成中间值MCD12Q1的类别 1 和 2 插出来变成 1.5等于凭空造了个不存在的类型。注意SPATIAL_SUBSET_TYPE INPUT_LAT_LONG时UL 是左上角、LR 是右下角顺序写反先 LR 后 UL不会报错但输出可能是一张空图或者尺寸为 0 的文件。5. 五个翻车现场MODIS 下载与处理里的避坑清单这一章记录的都是真实踩过、且反复有人踩的坑。每一条按“现象—原因—解决”写遇到对应症状可以直接对号入座。5.1 现象L1B 打开是一条歪的带子边缘发黑原因漏下 MOD03解决配对下载并做几何校正MOD021KM是沿轨条带数据它的每个像元在文件里是按扫描顺序排列的没有地理坐标只有仪器坐标。想让它变成正射影像必须配合MOD03地理定位文件里面存着每个像元的经纬度查找表。漏下MOD03的表现就是能打开、能出图、但影像斜着弯着边缘一圈黑。解决办法是在检索时把MOD03和MOD02按同一时刻成对下载文件名里的时间字段A2020032.0345必须完全一致。做几何校正时MRT 的 swath 模式或者 ENVI 的 GLT 校正都能用前者更快后者可控性更强。如果是做区域制图其实更建议直接绕开 L1B改用MOD09GA这类已经做过定位的 L2G 产品省掉整个几何校正环节。5.2 现象时间序列中间莫名少几天原因按北京时间检索 UTC解决检索窗口整体前移 8 小时并加宽按“2020-01-01 到 2020-12-31”去检索拿到的实际是 UTC 的这一整年对应北京时间 2020-01-01 08:00 到 2021-01-01 08:00。开头少了 8 小时结尾多了 8 小时中间还可能出现某天完全没有数据的情况——因为那颗星那天的过境时刻恰好跨了 UTC 日界而你按 UTC 日聚合时把它归到了另一天。解决办法是检索时把窗口前后各放宽一天下载完用第 2 章那个解析函数把采集时刻提出来统一换算成北京时间再做日聚合。另外L3 产品的时间标签是合成周期的起始日16 天合成的MOD13Q1标的是周期第一天不是中间那天做物候提取时这个偏移必须考虑进去否则整条曲线会平移 8 天。5.3 现象脚本读文件报Unable to open file或者返回None原因HDF4/HDF5 混淆、子数据集路径写错解决先枚举再取三种情况会长得一模一样。第一种是拿 h5py 去读 HDF4 文件报签名不匹配。第二种是gdal.Open()打开了主文件但忘了取子数据集ReadAsArray()返回None因为主文件本身不是栅格。第三种是子数据集路径里的网格名或数据集名拼错了一个字符gdal.Open()静默返回None不抛异常。排查顺序建议固定下来先用gdalinfo 文件名看它是 HDF4 还是 HDF5再用GetSubDatasets()打印全部路径从输出里复制而不是手打最后确认路径里的引号是英文双引号且文件路径部分用绝对路径。Windows 下的路径反斜杠还要转成正斜杠否则 GDAL 的虚拟文件路径解析会失败。5.4 现象拼接后接缝处有条带或者整块是 NoData原因先投影后拼接、或 UL/LR 顺序反了解决固定“先拼后投”顺序并检查角点接缝条带的根因是重采样顺序。每块瓦片单独重投影边缘像元在重采样时会因为邻域超出瓦片范围而被填成 NoData拼起来就是一条缝。正确顺序是先mrtmosaic在正弦投影下拼成完整一块再resample整体投影接缝处的邻域信息是完整的。整块 NoData 则多半是空间子集写错了。SPATIAL_SUBSET_UL_CORNER必须比SPATIAL_SUBSET_LR_CORNER的纬度大、经度小写反了交集为空MRT 会输出一个 0 字节或者全填充的文件而且退出码是 0看着像成功。我的习惯是在 prm 文件里把角点写成注释放在旁边改参数时肉眼对一遍。5.5 现象算出来的 NDVI 全是几百上千的整数原因没做定标解决读 scale_factor 并先掩膜后换算MODIS 的 L3 产品为了省存储把浮点物理量压成了Int16。NDVI 存的是真实值 × 10000地表温度存的是开尔文 × 50MOD11A2的scale_factor是 0.02反射率存的是真实值 × 10000。直接把数组拿去算或者画图得到的是一张数值范围在 -3000 到 10000 之间的图看着“有数据”其实完全不对。正确顺序是三步先按valid_range和_FillValue做掩膜再乘scale_factor加add_offset最后才做指数计算或统计。MOD13Q1的 NDVI 有效范围是 -2000 到 10000超出这个范围的都是填充值或者云污染像元。还有一点MOD11A2的LST_Day_1km在海洋和常年冰雪区是填充值 0如果直接参与均值统计会把结果拉低务必先掩掉。6. Albers 参数与批量流水线让每一次重投影都对得上投影参数这件事写一次和写一百次是两种心态。单次处理时随便本文还有配套的精品资源点击获取