企业官网建设流程全解析

简介：这份资源面向从事遥感与地球科学分析的科研人员和学生，提供基于MATLAB驱动MRT（MODIS Reprojection Tool）对MODIS数据进行批处理重投影与数据集提取的脚本工具，帮助解决大批量HDF-EOS格式数据在投影转换、区域裁剪和变量筛选中的效率问题。压缩包内共1个文件，为m脚本类型，体积约1KB，核心是用于自动化调用MRT命令行、指定输入输出路径与投影参数的批处理脚本，可一次性完成多产品重投影、镶嵌与感兴趣区域提取，减少重复手工操作。目前已有418人学习下载，适合具备一定遥感、GIS背景及MATLAB编程经验的用户参考，能快速搭建可复用的MODIS预处理工作流，将原始数据转换为GeoTIFF或ENVI等通用格式，便于后续在GIS软件中集成分析，也为气候研究、环境监测等场景下的数据准备提供实用起点。

1. 从一堆 HDF 到一张能用的投影栅格：MODIS 批处理到底卡在哪

如果你手头有一批 MODIS 的 HDF 文件，想把它们统一重投影、裁出研究区、再提取成 GeoTIFF 或表格，那你大概率已经听过 MRT 这个名字。MODIS Reprojection Tool 就是干这件事的老牌工具，配合一个批处理脚本，能把几十上百景影像从正弦投影转到常见的等经纬或阿尔伯斯投影，顺带做镶嵌和子集提取。听起来简单，但真正上手的人都知道，卡人的从来不是“会不会点按钮”，而是批处理参数文件写错一个字段、瓦片编号对不上、投影参数和分辨率不匹配，跑一晚上出来一堆空文件。这篇笔记就按我实际踩过的路，把 MRT 批处理重投影和数据集提取的完整链路拆开讲，从环境准备、参数文件写法、批处理脚本，到提取像元值时怎么避开坐标偏移的坑，尽量让第一次接触的人也能照着跑通，熟手也能对照检查自己的参数边界。

2. 先搞清楚 MRT 能做什么、不能做什么

2.1 MRT 的定位与替代方案对比

MRT 是专门为 MODIS 产品设计的重投影和子集工具，支持正弦投影到地理经纬度、阿尔伯斯等距、兰勃特等角、墨卡托等常见投影的转换，也能按经纬度范围或瓦片编号做空间子集，还能顺带做波段选择和格式转换。它的优势在于对 MODIS 的 HDF 结构理解得比较透，能自动识别波段、缩放因子和填充值，不用自己从头解析元数据。但它的局限也很明显：只认 MODIS 的 HDF 格式，其他遥感数据它不管；批处理靠参数文件驱动，没有图形界面那么直观；而且它依赖 Java 运行环境，不同版本对 Java 版本还有要求。

常见做法是，如果只是少量文件，用 MRT 的 GUI 跑一遍，把生成的参数文件存下来，再改成批处理模板。如果文件量大，或者需要集成到自动化流程里，就直接写参数文件加脚本循环。也有人用 GDAL 的gdalwarp替代，但 GDAL 对 MODIS 正弦投影的处理需要手动指定+proj=sinu +lon_0=0之类的参数，而且 HDF 里的子数据集要先用gdalinfo确认路径，对新手来说反而更绕。所以如果你的数据就是 MODIS 标准产品，MRT 仍然是省事的选择。

2.2 环境准备与文件目录组织

MRT 的安装包通常是一个压缩包，解压后里面有bin、data、doc等目录。Windows 下直接运行bin里的可执行文件，Linux 下需要给bin和lib里的文件加执行权限。我一般会把 MRT 放在一个不带空格和中文的路径下，比如D:\tools\MRT或/opt/MRT，因为参数文件里写路径时，空格和中文经常导致解析失败，这是血泪经验。

数据目录建议按产品类型和年份分开，比如MOD09A1_2020、MOD09A1_2021，每个目录下放原始的 HDF 文件。输出目录单独建一个，不要和输入混在一起，否则批处理循环时容易把输出文件又当成输入读进去。另外，MRT 在运行时会生成一些临时文件，确保输出目录有写入权限。

# Linux 下给 MRT 加执行权限 chmod +x /opt/MRT/bin/* chmod +x /opt/MRT/lib/* # 检查 Java 环境，MRT 一般需要 Java 8 或更高 java -version

这段命令做两件事：一是给 MRT 的可执行文件和库文件加上执行权限，否则运行时会报权限拒绝；二是确认 Java 版本，如果 Java 版本太低，MRT 启动时会直接闪退或者报UnsupportedClassVersionError。参数说明：chmod +x后面的路径要换成你实际的 MRT 安装路径，java -version输出的版本号如果低于 1.8，建议先升级 Java。

2.3 参数文件的结构与关键字段

MRT 的批处理参数文件是一个纯文本文件，通常以.prm结尾。它由若干行组成，每行一个键值对，格式是KEY = VALUE。核心字段包括：

这些字段里最容易翻车的是SPECTRAL_SUBSET的括号和数字个数。MODIS 不同产品的波段数不一样，比如 MOD09A1 有 7 个波段，SPECTRAL_SUBSET就要写 7 个数字，选中的波段写 1，不选的写 0。如果数字个数不对，MRT 会报错或者输出空文件。另一个坑是OUTPUT_PROJECTION_PARAMETERS，不同投影需要的参数个数不同，地理经纬度只需要一个参数（通常是 0），阿尔伯斯需要多个，写错了投影结果会扭曲。

3. 写一个能复用的批处理参数模板

3.1 从单文件参数到批处理模板

最稳妥的做法是先用 MRT 的 GUI 跑一个文件，把生成的.prm文件保存下来，然后用文本编辑器打开，把里面和具体文件相关的字段改成变量占位符。比如INPUT_FILENAME和OUTPUT_FILENAME改成{input}和{output}，其他字段保持不变。这样你就得到了一个模板，后续用脚本替换占位符生成每个文件的参数文件。

我一般会保留一份template.prm，里面除了输入输出路径，其他都是固定值。如果研究区范围固定，SPATIAL_SUBSET_UL_CORNER和SPATIAL_SUBSET_LR_CORNER也写死；如果每个文件的范围不一样，就把这两个字段也做成变量。重采样方法根据数据类型选：分类数据用NEAREST_NEIGHBOR，连续数据用BILINEAR或CUBIC_CONVOLUTION。输出投影如果选GEOGRAPHIC，OUTPUT_PROJECTION_PARAMETERS通常写( 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 )，这是 MRT 的固定格式，不用深究每个 0 的含义，照抄就行。

3.2 用 Python 生成参数文件并调用 MRT

下面是一个 Python 脚本，它会遍历输入目录下的所有 HDF 文件，为每个文件生成一个.prm参数文件，然后调用 MRT 的resample命令执行重投影。

import os import subprocess from string import Template # 配置路径 mrt_bin = "/opt/MRT/bin/resample" # MRT 可执行文件路径 input_dir = "/data/MOD09A1_2020" # 输入 HDF 目录 output_dir = "/data/output" # 输出目录 template_file = "template.prm" # 参数模板文件 # 读取模板 with open(template_file, "r") as f: template_content = f.read() # 研究区范围（经纬度） ul_lon, ul_lat = 110.0, 40.0 lr_lon, lr_lat = 120.0, 30.0 # 遍历 HDF 文件 for filename in os.listdir(input_dir): if not filename.endswith(".hdf"): continue input_path = os.path.join(input_dir, filename) # 输出文件名去掉 .hdf 后缀，加 _reproj base_name = os.path.splitext(filename)[0] output_path = os.path.join(output_dir, base_name + "_reproj") # 替换模板中的占位符 prm_content = Template(template_content).substitute( input=input_path, output=output_path, ul_lon=ul_lon, ul_lat=ul_lat, lr_lon=lr_lon, lr_lat=lr_lat ) # 写入临时参数文件 prm_file = os.path.join(output_dir, base_name + ".prm") with open(prm_file, "w") as f: f.write(prm_content) # 调用 MRT cmd = [mrt_bin, "-p", prm_file] result = subprocess.run(cmd, capture_output=True, text=True) if result.returncode != 0: print(f"处理 {filename} 失败：{result.stderr}") else: print(f"处理 {filename} 完成")

这段脚本的逻辑很直接：读模板、替换占位符、写参数文件、调 MRT。关键点在于Template的占位符要和模板文件里的变量名一致，比如模板里写INPUT_FILENAME = $input，脚本里就传input=input_path。参数说明：mrt_bin要换成你实际的resample路径，Windows 下是resample.exe；ul_lon等四个变量是研究区范围，根据你的需求改；subprocess.run的capture_output=True会捕获 MRT 的输出，方便排查错误。如果 MRT 报错，先看result.stderr里的信息，常见的是路径不存在、波段数不对、投影参数格式错误。

3.3 模板文件示例与字段解释

下面是一个针对 MOD09A1 的模板文件示例，输出地理经纬度投影，选第 1 到第 3 波段，重采样用最近邻。

INPUT_FILENAME = $input SPECTRAL_SUBSET = ( 1 1 1 0 0 0 0 ) SPATIAL_SUBSET_TYPE = INPUT_LAT_LONG SPATIAL_SUBSET_UL_CORNER = ( $ul_lat $ul_lon ) SPATIAL_SUBSET_LR_CORNER = ( $lr_lat $lr_lon ) OUTPUT_FILENAME = $output RESAMPLING_TYPE = NEAREST_NEIGHBOR OUTPUT_PROJECTION_TYPE = GEOGRAPHIC OUTPUT_PROJECTION_PARAMETERS = ( 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 ) DATUM = WGS84 OUTPUT_PIXEL_SIZE = 0.0025

注意SPATIAL_SUBSET_UL_CORNER的顺序是纬度在前、经度在后，这是 MRT 的约定，写反了会裁到地球另一边。OUTPUT_PIXEL_SIZE的单位是度，0.0025 度大约是 250 米，对应 MOD09A1 的原始分辨率。如果输出投影换成阿尔伯斯，OUTPUT_PROJECTION_PARAMETERS需要填更多参数，具体格式可以在 MRT 的文档里找到，或者用 GUI 生成一次后照抄。

4. 重投影后的数据提取与坐标对齐

4.1 用 GDAL 读取重投影结果并提取像元值

MRT 输出的是 GeoTIFF，每个波段一个文件，文件名通常带波段后缀。提取像元值时，可以用 GDAL 的 Python 绑定，也可以用rasterio。我一般用rasterio，因为它的 API 更简洁。下面是一个提取指定经纬度处像元值的例子。

import rasterio import numpy as np # 打开重投影后的 TIFF tif_path = "/data/output/MOD09A1_2020_001_reproj.tif" with rasterio.open(tif_path) as src: # 目标经纬度 lon, lat = 115.0, 35.0 # 将经纬度转换为行列号 row, col = src.index(lon, lat) # 读取该位置的像元值 value = src.read(1, window=((row, row+1), (col, col+1))) print(f"像元值：{value[0, 0]}") # 检查坐标参考系 print(f"CRS: {src.crs}") print(f"分辨率: {src.res}")

这段代码的核心是src.index(lon, lat)，它把地理坐标转成栅格的行列号。注意src.index返回的是(row, col)，不是(col, row)，顺序别搞反。src.read(1, window=...)读取第 1 波段在指定窗口内的值，窗口大小是 1x1，所以返回一个二维数组，取[0, 0]就是单个像元值。如果src.crs不是EPSG:4326，说明重投影时投影类型选错了，需要检查参数文件里的OUTPUT_PROJECTION_TYPE。

4.2 批量提取与表格输出

如果要在多个点位上批量提取，可以把点位存成 CSV，循环读取每个点的经纬度，然后写入结果表格。

import csv import rasterio points = [(115.0, 35.0), (116.0, 36.0), (117.0, 34.0)] tif_path = "/data/output/MOD09A1_2020_001_reproj.tif" with rasterio.open(tif_path) as src: with open("extracted.csv", "w", newline="") as f: writer = csv.writer(f) writer.writerow(["lon", "lat", "value"]) for lon, lat in points: row, col = src.index(lon, lat) value = src.read(1, window=((row, row+1), (col, col+1)))[0, 0] writer.writerow([lon, lat, value])

这个脚本把每个点的经纬度和提取值写进 CSV。参数说明：points列表里是经纬度元组，tif_path是重投影后的文件路径。如果点位很多，建议先检查每个点是否在栅格范围内，src.index对超出范围的点会返回超出行列号的值，读取时会报错。可以在循环里加一个判断：if 0 <= row < src.height and 0 <= col < src.width。

4.3 坐标偏移的排查方法

重投影后提取的值对不上，最常见的原因是坐标参考系不一致。比如 MRT 输出的是地理经纬度，但rasterio打开时如果没正确识别 CRS，src.index就会按默认的像素坐标算，导致偏移。排查方法是先打印src.crs和src.bounds，确认 CRS 是EPSG:4326，边界范围和研究区一致。如果 CRS 不对，可以用rasterio.warp重新投影，或者检查 MRT 参数文件里的DATUM和OUTPUT_PROJECTION_TYPE是否匹配。

另一个常见问题是OUTPUT_PIXEL_SIZE设得太大或太小，导致重采样后像元值被平滑或出现锯齿。对于分类数据，OUTPUT_PIXEL_SIZE最好和原始分辨率一致，重采样方法用NEAREST_NEIGHBOR；对于连续数据，可以适当放大像元，但不要超过原始分辨率的 2 倍，否则信息损失明显。

5. 避坑与常见问题排查

5.1 批处理中途报错但不知道哪个文件出问题

现象：脚本跑了一晚上，早上发现输出目录里只有一部分文件，日志里一堆错误但分不清是哪个文件。原因：MRT 的错误信息有时只输出到标准错误，而脚本没有捕获或者没有打印文件名。解决：在循环里给每个文件加一个标识，比如print(f"正在处理 {filename}")，并且在subprocess.run之后检查returncode，把stderr和文件名一起打印。更稳妥的做法是每个文件处理完后检查输出文件是否存在且大小大于 0，否则记录到失败列表。

5.2 输出文件为空或只有几 KB

现象：MRT 运行没有报错，但输出的 TIFF 文件只有几 KB，打开后全是 NoData。原因：SPATIAL_SUBSET_UL_CORNER和SPATIAL_SUBSET_LR_CORNER的经纬度写反了，或者研究区范围和影像实际覆盖范围没有交集。解决：先用gdalinfo查看原始 HDF 的经纬度范围，确认研究区在影像范围内。另外检查SPATIAL_SUBSET_TYPE是否设成了INPUT_LAT_LONG，如果设成INPUT_PIXEL_LINE，但填的是经纬度，就会裁到错误位置。

5.3 波段选择错误导致输出波段数不对

现象：明明选了 3 个波段，输出却只有 1 个或者 7 个。原因：SPECTRAL_SUBSET的括号里数字个数和产品波段数不匹配。比如 MOD09A1 有 7 个波段，你只写了 3 个数字，MRT 会按默认行为处理，可能只输出第一个波段。解决：查清楚产品的波段数，SPECTRAL_SUBSET里写够数字，选中的写 1，不选的写 0。不确定的话，先用 GUI 选一次，看它生成的参数文件里SPECTRAL_SUBSET是怎么写的。

5.4 投影参数格式错误导致重投影失败

现象：MRT 报错Invalid projection parameters或者输出结果扭曲。原因：OUTPUT_PROJECTION_PARAMETERS的括号、数字个数或顺序不对。不同投影需要的参数不同，地理经纬度通常 15 个 0，阿尔伯斯需要更多。解决：用 GUI 生成一次对应投影的参数文件，直接复制OUTPUT_PROJECTION_PARAMETERS那一行。不要手动改数字，除非你清楚每个参数的含义。

5.5 路径中有空格或中文导致解析失败

现象：MRT 报错File not found，但文件明明存在。原因：参数文件里的路径包含空格或中文，MRT 解析时把空格当成了分隔符。解决：把所有路径改成不带空格和中文的英文路径，或者用引号把路径括起来。我一般会在项目开始前就把数据目录改成纯英文无空格的路径，省得后面折腾。

6. 进阶：把 MRT 批处理嵌入自动化流程的几个技巧

如果你已经能跑通单次批处理，下一步可以考虑把它嵌入更自动化的流程。我自己的习惯是用一个主脚本做调度，MRT 只负责重投影和子集，后续的提取、统计、入库用 Python 单独处理。这样分工明确，出问题也容易定位。

一个具体的技巧是：不要每次都用 MRT 重新生成参数文件，而是把参数文件模板和文件清单分开管理。文件清单可以是一个 CSV，里面包含输入路径、输出路径、研究区范围等字段，主脚本读 CSV 生成参数文件。这样修改研究区范围时只需要改 CSV，不用动脚本代码。

另一个技巧是并行处理。MRT 本身是单线程的，但你可以用 Python 的multiprocessing或者concurrent.futures同时跑多个 MRT 进程。注意不要开太多，一般设为 CPU 核心数的一半，因为 MRT 运行时也会占内存。下面是一个简单的并行示例：

from concurrent.futures import ProcessPoolExecutor import subprocess def run_mrt(prm_file): cmd = ["/opt/MRT/bin/resample", "-p", prm_file] result = subprocess.run(cmd, capture_output=True, text=True) return prm_file, result.returncode, result.stderr # 假设 prm_files 是已经生成好的参数文件列表 with ProcessPoolExecutor(max_workers=4) as executor: for prm_file, code, err in executor.map(run_mrt, prm_files): if code != 0: print(f"{prm_file} 失败：{err}")

这段代码用ProcessPoolExecutor并行调用 MRT，max_workers=4表示同时跑 4 个进程。参数说明：prm_files是参数文件路径列表，run_mrt函数返回参数文件名、返回码和错误信息。并行处理能显著缩短大批量文件的处理时间，但要注意磁盘 I/O 和内存占用，如果输出目录在同一块磁盘上，并行写入可能会成为瓶颈。

最后说一个验证方法：重投影完成后，随机抽几个文件，用gdalinfo检查 CRS、分辨率和范围，再用rasterio提取几个已知点的值，和原始 HDF 里的值对比。如果偏差在合理范围内（比如最近邻重采样应该完全一致），说明流程没问题。如果偏差很大，回头检查投影参数和重采样方法。这个验证步骤我每次都会做，虽然麻烦，但能避免后面用错数据。希望帮到你。

本文还有配套的精品资源，点击获取