简介《MODIS数据综合处理软件 V1.0》是一套面向遥感、GIS及生态环境研究人员的桌面工具集成MODIS陆面产品批量处理、多核并行计算与可视化分析功能支持NDVI、EVI、ET/PET、LST、LAI、GPP/NPP等常用产品适用于科研数据预处理、资源环境监测与教学演示等场景。软件附详细使用手册安装包与手册打包为一个zip压缩包共2000个文件文件类型以JavaScript、JSON、HTML等前端程序为主兼有Python脚本、PDF文档、CSS样式、XML配置及Markdown笔记能覆盖软件运行、用户查阅与二次开发等多种用途压缩包整体176.28MB目录结构清晰解压即可部署。已有140人学习下载。借助该软件用户无需编写复杂代码即可完成多文件批量导入、一键处理与结果出图多核并行机制可显著缩短大规模数据计算时间可视化工具则便于快速把握地表参数的时空变化趋势显著提升MODIS数据处理与制图效率是相关领域数据分析流程中一款实用型辅助软件。1. 用 MODIS 数据综合处理软件之前先回答“能不能直接用”这个问题把“MODIS数据综合处理软件 V1.0”这串字拆开看它不像某个商业软件那么神秘——它想解决的是每个用过 MODIS 数据的人都撞过的墙从 NASA 下载下来的 HDF 文件用 GIS 打开一片黑或者看到一个波段根本不知道是什么投影。我最早拿到 MOD11A1 地表温度产品时直接把它当普通栅格丢进 ArcGIS出来的图数值上万分我一度以为卫星坏了后来才意识到是没做综合处理。所谓“综合处理软件 V1.0”在我的理解里就是一套把下载、拼接、重投影、裁剪、质量控制和批量导出的流程固化成脚本和参数模板的工具箱。它面向的人很具体做地表温度、植被指数、水体或土地利用变化研究又不想每次都在 GUI 里点半小时鼠标的从业者和研究生。这篇文章就按这个思路从数据格式认知讲到落地代码再讲清楚热词里那个灵魂拷问——MODIS 下载的地表温度数据到底能不能直接用。答案是能但绝对不能直接画图。2. 先认清 MODIS 的“脾气”HDF 格式、产品家族与正弦投影2.1 HDF 文件结构和子数据集为什么打开一个文件有十几个图层MODIS 的 L2/L3 标准产品是 HDF4 格式后缀 .hdf。HDF 本身是一个容器里面可以装多个科学数据集SDS专业说法叫子数据集。你用普通软件打开一个 MOD11A1.hdf看到的是“栅格”和“属性表”但其实里面装了两类东西一类是实际的地理数据例如白天地表温度LST_Day_1km、夜间地表温度、质量控制波段QC_Day、观测时间、视角等另一类是元数据包括投影信息、尺度因子、有效值范围、填充值。这两个东西混在一个文件里是新手最容易困惑的地方。用 GDAL 的命令行工具看一眼就知道里面有什么gdalinfo MOD11A1.A2020001.h23v04.006.2020001030000.hdf输出会列出类似这样的子数据集Subdatasets: SUBDATASET_1_NAMEHDF4_EOS:EOS_GRID:MOD11A1.A2020001.h23v04.006.2020001030000.hdf:MOD_Grid_Daily_1km_LST:LST_Day_1km SUBDATASET_2_NAMEHDF4_EOS:EOS_GRID:MOD11A1.A2020001.h23v04.006.2020001030000.hdf:MOD_Grid_Daily_1km_LST:QC_Day这里每一行是一个独立子数据集。LST_Day_1km 是数值型地表温度灰度值QC_Day 是 8 位无符号整数质量控制位域。我的习惯是先把这两层拆出来做预处理而不是一上来就放到拼接工具里——因为 QC 在你做重采样和裁剪时会被插值污染一旦污染后续质量筛选全是错的。拆层的代码后面会说这里先建立概念所谓“处理一个 MODIS 文件”第一步永远是“从 HDF 里把这个子数据集正确地抠出来”。2.2 MODIS 产品家族你下载的到底是什么东西MODIS 数据有几十种产品命名规则是 MODTerra或 MYDAqua开头后接三位产品号。常见的一些做地表覆盖和气候研究会用到的产品我整理如下产品 ID分辨率时间合成主要内容MOD11A1 / MYD11A11 km每日地表温度与发射率含质量控制波段MOD11A2 / MYD11A21 km8 天地表温度 8 天合成LST_Day_1km 是均值MOD13A1500 m16 天NDVI / EVI 植被指数MOD09GA500 m每日地表反射率分波段MOD141 km每日热异常 / 火点MCD43A4500 m每日Nadir BRDF 调整反射率区分这些产品很重要因为“综合处理”的流程会根据产品不同而改变。比如 MOD13A1 本身是 16 天合成处理时不需要再做时间平均而 MOD11A1 是瞬时值一天里 Terra 和 Aqua 各过境两次处理时如果要把白天和晚上分开就得同时读文件名里的时间字段和产品内部的时间 SDS。我在处理地表温度时一般会把日产品和 8 天合成产品分开两条流水线因为前者适合做极端温度事件分析后者适合做气候态平均。除了产品号文件名里还有一个关键信息分片编号。MODIS 标准产品采用正弦投影的分片网格tile grid整个地球被横竖切成若干块文件名里 h23v04 这样的编号就代表横向第 23 片、纵向第 4 片。中国的陆地基本覆盖在 h23v04、h24v04、h25v04、h26v04、h27v04 等几片上跨片研究区必须把所有覆盖到的片子都下载下来再做拼接不能只下中间一片。2.3 正弦投影是怎么回事为什么 WGS84 坐标的人容易栽在这里MODIS 标准产品的默认投影是正弦投影SinusoidalGDAL 能识别但 ArcGIS 的老版本和一些在线平台不一定认打开后经常出现整个图层定位到海里或者影像拉伸变形的现象。正弦投影本质上是一种等积伪圆柱投影赤道附近变形小高纬度变形大。它在全球尺度做统计时面积保持得好但到了区域尺度想要和矢量边界、气象站点或者其他卫星资料叠加就必须重投影到你自己的目标坐标系。常见做法有两种一是转成 WGS84 地理坐标系EPSG:4326适合做全球或大区域制图二是转成 UTM根据研究区经度选带号适合做局部定量分析。我的经验是做地表温度时间序列用 WGS84 就够了因为你不涉及面积算量做森林覆盖或蒸散发涉及面积统计的必须转 UTM 或者转 Albers 等积投影。这一步看着简单但如果在拼接前不做后面裁剪出来的数据就会出现边界错位。所以综合处理软件 V1.0 的第一版核心应该包含四件事HDF 子数据集提取、拼接、重投影、裁剪。这四个步骤的顺序不能乱。先拼接再重投影比先重投影再拼接稳定得多因为原始 tile 之间是严格无缝的一旦各自重投影两个相邻 tile 的边界会因为重采样产生细缝和重叠后面拼起来要处理羽化和接边线属于给自己找罪受。3. 用 GDAL 搭一条预处理流水线从 HDF 抠数据到拼接与裁剪3.1 下载与目录组织给 V1.0 定一个不会被自己搞乱的目录结构我一般在 NASA Earthdata Search 上下载数据选定时间范围和 tile 号后一次把覆盖研究区的所有 tile 都勾选下载。下载后先做文件归档否则处理到一半你会发现某个日期的文件缺失蹲在电脑前抓瞎。推荐的目录结构是MODIS/ raw/ # 原始 HDF 文件按 tile 和时间子目录 lst_day/ # 拆分出的 LST 子数据集 lst_night/ qc/ mosaic/ # 拼接后的 GeoTIFF clipped/ # 按研究区裁剪后的结果 log/ # 处理日志这个结构看着简单但能避免一个很常见的翻车不同时间的文件混在一起跑了批量脚本后发现 2020 年第 200 天和 2021 年第 200 天拼在了一起。我的做法是把 HDF 文件名里的日期提取出来写进输出文件名格式统一为LST_YYYYDOY.tif比如LST_2020200.tif这样后面做时间序列排序时按文件名排序就是按时间排序。3.2 抠子数据集gdal_translate 提取单波段这一步是流水线的第一环。用 gdal_translate 把 HDF 里的特定子数据集拉出来转成 GeoTIFF。这里要特别叮嘱子数据集全名很长每次手敲都会疯掉我都是先跑一遍 gdalinfo 把名字复制出来或者用脚本动态获取。gdal_translate -of GTiff \ HDF4_EOS:EOS_GRID:MOD11A1.A2020001.h23v04.006.2020001030000.hdf:MOD_Grid_Daily_1km_LST:LST_Day_1km \ ./lst_day/LST_day_2020001.tif这条命令做的事是打开 HDF 文件定位到 LST_Day_1km 子数据集写成单个 GeoTIFF。输出文件保留了原始投影信息和地理变换参数。参数说明-of GTiff指定输出格式坐标和投影信息默认继承源文件所以此时如果直接打开看到的还是正弦投影。我一般在这一步不做任何重采样只为把数据从 HDF 容器里救出来避免每次都要处理 HDF 嵌套结构。3.3 按日期拼接gdalwarp 还是 gdal_merge.py拼接有很多工具但处理 1km 分辨率的地表温度产品我推荐用 gdalwarp因为它能同时完成镶嵌和重投影一步到位。gdal_merge.py 也可以拼但它只做简单拼接不处理投影不一致的问题两片不同日期或不同轨道的影像只要投影参数有一丁点差异拼完就会出现明显的接缝。gdalwarp -t_srs EPSG:4326 -r bilinear -overwrite \ ./lst_day/LST_day_2020001_tile1.tif \ ./lst_day/LST_day_2020001_tile2.tif \ ./mosaic/LST_day_2020200_wgs84.tif这里把覆盖研究区的两个 tile 拼接并重投影到 WGS84。-t_srs EPSG:4326是目标坐标系-r bilinear是重采样方法。地表温度这种连续变量用双线性插值比最近邻更平滑但如果你在处理分类产品或者需要严格保留原始像元值的场景比如火点检测要改用-r near。重采样方法选错是新手常踩的坑后面避坑章节还会展开。再看一遍这条命令的逻辑输入是两个 tile 的白天 LST GeoTIFF输出是拼接好的 WGS84 GeoTIFF。gdalwarp 会自动读取每个文件的覆盖范围找到重叠区按默认策略做叠加。对于 MOD11A1 这种每天多轨数据有时同一天相邻轨道的重叠区数值不一致我的做法是拼完后用 QC 做掩膜后处理而不是在拼接时加羽化参数。羽化会让真正的高温异常被抹平得不偿失。3.4 按矢量边界裁剪把数据切到研究区拼接完成了如果研究区只是一个城市或者一个流域下一步是裁剪。最常见的方式是用一个 shapefile 边界做掩膜裁剪gdalwarp -cutline ./shp/study_area.shp -crop_to_cutline \ -t_srs EPSG:4326 -dstnodata -9999 \ ./mosaic/LST_day_2020200_wgs84.tif \ ./clipped/LST_day_2020200_study.tif-cutline指定边界矢量-crop_to_cutline让输出范围严格跟随边界范围。这里有一个细节边界矢量的坐标系最好和目标栅格一致如果不一致gdalwarp 会自动做坐标变换但你要确保 shapefile 里定义了正确的投影信息否则裁剪结果可能在研究区外或者完全空白。到这一步一个 MODIS 文件已经从 HDF 变成了一个带正确投影、按边界裁剪好的 GeoTIFF。但这才完成“综合处理”的一半因为还没有做质量控制。而缺少质量控制的地表温度数据就是热词里问到的“能不能直接用”的根源。4. 地表温度数据到底能不能直接用质量控制、尺度因子与无效值4.1 灰度值 ≠ 真实温度0.02 这个尺度因子是怎么用的先直接回答热词那个问题modis 下载地表温度数据可以直接用吗答案是“不可以直接画图”原因有两个。第一LST 产品里存储的是灰度值DN 值不是物理温度。MOD11A1 的白天地表温度波段灰度值乘以 0.02 才得到开尔文温度再减去 273.15 才等于摄氏度。第二影像里不全是有效观测有云遮挡、有填充值、有质量不佳的像元。如果你不处理这两件事做出来的地表温度图会同时包含大量异常低值云顶温度被当成地表温度和诡异的高值。读取 LST 灰度值的脚本我一般这样写from osgeo import gdal import numpy as np # 打开 HDF 中的 LST 子数据集 lst_ds gdal.Open( HDF4_EOS:EOS_GRID:MOD11A1.A2020001.h23v04.006.2020001030000.hdf: MOD_Grid_Daily_1km_LST:LST_Day_1km ) lst_dn lst_ds.ReadAsArray().astype(np.float32) # 尺度因子 0.02单位是开尔文 lst_k lst_dn * 0.02 lst_c lst_k - 273.15这里lst_dn是原始灰度值lst_c是摄氏温度。有一个易错点填充值 0 乘以 0.02 等于 00 再减 273.15 成了 -273.15这个值会严重拉低后续统计的最低值。所以必须先剔除填充值再算物理量。正确的顺序是先找填充值和无效值再乘尺度因子最后才做后续运算。4.2 QC 波段质量不是用眼睛看的是拿 bit 算出来的MOD11A1 质量控制波段QC_Day是 8 位无符号整数它不是一个“0 到 255 的等级分数”而是每两位 bit 记录一种质量信息的位域。最常见的做法是提取 bit0 和 bit1 判断综合质量等级0 表示质量好1 表示质量一般2 表示云遮挡3 表示云阴影。很多人拿到 QC 后直接做阈值筛选比如保留 QC 64这是完全错误的因为 QC64 的二进制是 01000000bit0-1 是 0质量好却被误杀了。正确的 QC 筛选方式是用位运算from osgeo import gdal import numpy as np qc_ds gdal.Open( HDF4_EOS:EOS_GRID:MOD11A1.A2020001.h23v04.006.2020001030000.hdf: MOD_Grid_Daily_1km_LST:QC_Day ) qc qc_ds.ReadAsArray() # bit0-1: 0良好, 1一般, 2云, 3云阴影 quality qc 0b11 valid_mask (quality 0) | (quality 1)取qc 0b11位与运算得到低两位的十进制值。这样好质量的像元即使它的高比特位有其他信息比如云检测标记、日/夜标记也不会被误过滤。这一点是地表温度处理里的分水岭理解它之后你再去读其他 MODIS 产品的 QA 文档会发现套路都一样。4.3 拼起来的完整处理温度、质量、掩膜一次搞定把上面几段组合成一个完整的最小流程from osgeo import gdal import numpy as np product_path ( HDF4_EOS:EOS_GRID:MOD11A1.A2020001.h23v04.006.2020001030000.hdf: MOD_Grid_Daily_1km_LST: ) lst_dn gdal.Open(product_path LST_Day_1km).ReadAsArray().astype(np.float32) qc gdal.Open(product_path QC_Day).ReadAsArray() # 第一步根据 QC 生成有效像元掩膜 good_pixels ((qc 0b11) 0) | ((qc 0b11) 1) # 第二步剔除填充值DN 值为 0 是无效观测 good_pixels (lst_dn ! 0) # 第三步尺度因子换算 无效值抑制 lst_k np.where(good_pixels, lst_dn * 0.02, np.nan) lst_c lst_k - 273.15 # 统计去看一眼确认数值范围合理 print(温度范围摄氏度, np.nanmin(lst_c), ~, np.nanmax(lst_c))过程的逻辑是先算质量控制得到哪些像元可信再剔掉填充值最后才做单位换算。顺序不能错如果先换算再筛选填充值 0 会污染温度统计。np.where在这里把无效像元直接置为np.nan后续做任何统计运算都不会再被它干扰。到这一步地表温度数据才算“能用了”。数据科学家经常说 Garbage in garbage out在 MODIS 这里不过滤 QC 就是 Garbage in garbage out 的说明书级案例。用直方图检查一下处理后的lst_c正常的 LST 白昼应该在 -20 到 60 摄氏度之间如果出现大量负值往往是云遮挡没滤干净如果出现大量超过 60 摄氏度的像元先检查是不是裸土或者火点再检查尺度因子是不是漏乘了。5. 综合处理中的 4 个常见坑现象、原因和排查顺序5.1 拼接后的影像出现明显的十字接缝或暗带现象用 gdalwarp 拼完两个 tile 后影像中间出现一条清晰的分界线一边亮一边暗或者在重叠区出现几何错位。原因这种情况绝大多数不是算法问题而是输入的两片影像日期不一致。MODIS 逐日产品每天由多次过境拼接而成云覆盖不同地表温度也会随时段变化。你下数据时选了同一天但下载到的是不同轨道来源的 tile它们虽然日期相同但过境时刻差了一个多小时正午地表温度变化很快拼在一起自然有色差。解决先查两个 tile 的元数据确认过境时间文件名里有精确到秒的时间字段。处理地表温度时同一天的 Terra 白天数据如果过境时间相差超过 30 分钟我一般直接放弃拼接改用 MOD11A2 的 8 天合成产品或者把研究区限定在单一过境覆盖范围内。拼接前额外检查一遍文件时间比拼完再调试快得多。5.2 QC 波段用阈值筛选后好的像元全被删了现象用qc 64之类的条件筛完影像变成大片空洞连夏天晴空万里的时候都是碎的。原因QC 是按位域存储的不是数值越高质量越差。MOD11A1 的 QC_Day 为 8 位高 bit 位存放云检测标志、邻近云等信息数值大小和质量好坏没有单调关系。我见过有人用qc 0筛数据结果影像只剩零星十几个像元因为大部分像元的 bit2-3 或 bit4-5 并不为 0但它们本身不代表质量问题。解决永远用位与运算。qc 0b11取低两位再判断是否为 0 或 1。如果要做更严格的质量控制再去看产品文档里 bit2-3 和 bit4-5 的定义不要用十进制阈值。5.3 温度图像数值在几万甚至十几万明显不是温度现象做完处理输出影像的像元值在 10000 到 30000 之间画出来像是外星地貌。原因忘了乘尺度因子直接把 DN 值当物理量使用了。MOD11A1 的 LST 灰度值通常在 7500 到 14000 之间乘 0.02 后才是 150 到 280 开尔文。这是热词“能不能直接用”里最典型的错误。解决乘尺度因子是硬性步骤写进任何脚本的前三行。我的习惯是把scale 0.02初始化在脚本顶部并加注释说明它是 MOD11A1 的固定属性。处理完输出后立刻用np.nanmin和np.nanmax打印数值范围范围不对就停下来查不要等到画图才暴露。5.4 按边界裁剪后影像全部空白或者偏移了几百公里现象用-cutline裁剪完输出的 GeoTIFF 打开是一片黑或者数据在研究区边界外像是一块被平移过的影像。原因两个坐标参考不一致。最常见的是栅格还是正弦投影而 shapefile 是 WGS84 经纬度虽然 gdalwarp 会自动做投影变换但如果 shapefile 缺少.prj文件GDAL 不知道它的坐标系就会按经纬度数值直接当坐标用结果裁剪范围跑到非洲或者大海。解决先跑gdalinfo shapefile.shp确认矢量有投影信息若没有.prj在 GIS 软件里给 shapefile 手动指定 WGS84。裁剪前先gdalinfo看一眼栅格范围和矢量范围是否大致吻合在同一个经纬度体系下两者的 extent 应该存在重叠。这一步只要养成习惯几乎可以杜绝裁剪空白的坑。5.5 时间序列里某一天的温度断崖式偏低其他天都正常现象做 2020 年逐日 LST 曲线大部分日期在 20 到 35 度之间某一天突然掉到 -5 度画出来的曲线像心电图。原因那天的研究区恰好被云覆盖而且 QC 筛选条件太宽保留了 quality1 或 2 的数据。MOD11A1 的 QC 位域里云掩膜信息在高 bit 位仅用qc 0b11会保留“一般质量”的像元而这些像元里很多是云边缘。解决做时间序列时质量控制收紧到quality 0再加一道云掩膜检查。具体做法是读 MOD11A1 附属的云掩膜波段如果产品提供或者直接用 MOD35 云产品做二次过滤。我处理长时间序列时还有个习惯对每个日期计算云覆盖比例云覆盖超过 40% 的日期直接标记为缺测不参与平均和统计。宁可少一天数据也不要一个假值混进序列里。6. 把 V1.0 变成自己的批处理骨架脚本模板与验证方法到这一步你已经有了完整的单日处理链路。把它固化成批处理骨架才算得上“软件”。我通常会把处理流程抽象成一个 Python 脚本输入是 HDF 文件列表输出是裁剪后的温标 GeoTIFF 加一张 QC 概览图。骨架结构如下from osgeo import gdal import numpy as np import glob def process_one_day(hdf_path, output_dir, cutline_shp): # 1. 列出 HDF 内所有子数据集 ds gdal.Open(hdf_path) subdatasets ds.GetSubDatasets() # 2. 按名称定位 LST 和 QC 层 def find_sds(keyword): for name, desc in subdatasets: if keyword in name: return name return None lst_name find_sds(LST_Day_1km) qc_name find_sds(QC_Day) # 3. 读取数据并做 QC 筛选与尺度换算 lst_dn gdal.Open(lst_name).ReadAsArray().astype(np.float32) qc gdal.Open(qc_name).ReadAsArray() good ((qc 0b11) 0) (lst_dn ! 0) lst_c np.where(good, lst_dn * 0.02 - 273.15, np.nan) # 4. 输出 GeoTIFF原始投影后面统一交给 gdalwarp 处理 driver gdal.GetDriverByName(GTiff) out_path f{output_dir}/LST_C_{lst_name.split(:)[-1]}.tif out_ds driver.Create(out_path, lst_c.shape[1], lst_c.shape[0], 1, gdal.GDT_Float32) out_ds.GetRasterBand(1).WriteArray(lst_c) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_ds.FlushCache() return out_path # 批量入口遍历某个日期的所有 tile for hdf in sorted(glob.glob(./raw/MOD11A1.A2020200.*.hdf)): process_one_day(hdf, ./lst_day, ./shp/study_area.shp)这段脚本把每天多个 tile 从 HDF 中拆出并换算成摄氏温度输出文件仍保留正弦投影。后续再用 gdalwarp 统一做投影转换和批量拼接。脚本里的find_sds函数解决了一个很实际的问题HDF 子数据集全名在不同版本的 MODIS 产品里可能有细微差异按关键字定位比硬编码全名更稳。骨架搭好后验证环节不能省。我的验证习惯是用气象站点的实测气温数据做参考但注意LST 是地表温度不是 1.5 米气温两者差值白天可能达到 5 到 15 摄氏度夜间差值小很多。比较时主要看两点一是 LST 和站点气温的相关系数是否显著二是高温期、低温期出现的日期是否对齐。如果你处理的是裸土区域LST 白天波动幅度会比气温大得多不要因为这个就觉得数据处理错了。还有一个我花过不少时间才养成的习惯每次批量处理完导出一张随意日期的 LST 分布图和对应 QC 掩膜图用眼睛扫一遍。MODIS 处理里的很多问题——投影错、尺度因子漏乘、云没滤净——在这种抽查下会在十分钟内暴露比事后发现数据不可用再重跑强得多。V1.0 的代码不追求一次写对追求的是出问题时能快速定位到是哪个环节。毕竟这行当的常态就是前面省下的检查时间后面都会以通宵重算的方式还回去。希望这篇笔记能帮你少踩几个我已经踩过的坑。本文还有配套的精品资源点击获取