遥感数据处理这件事最怕的不是算法难而是重复劳动。我最早接触MODIS数据的时候一个研究区、一年12个月、4种产品光是HDF文件的拼接、投影、裁剪、格式转换就能耗掉整整两天中间还得盯着MRT那个老界面反复点选参数稍不留神就报错重来。后来被逼得没办法硬是把整套流程用Python脚本串了起来配合ArcGIS做后处理和出图现在同样的工作量压缩到半小时以内而且可以批量跑、无人值守。这篇内容就是把这套流程完整拆开讲清楚从HDF原始文件到NDVI、ET成品图每一步为什么这么做、参数怎么定、坑在哪里都会讲到。适合做生态遥感、农业监测、水文分析方向的同学也适合任何需要批量处理栅格数据的GIS从业者参考。1. 为什么放弃MRT转向脚本化处理1.1 MRT的便利与它的天花板MRTMODIS Reprojection Tool在早期确实是处理MODIS数据的标配工具图形界面友好支持批量导入HDF、选择波段、设定投影和输出格式。对于偶尔处理几个文件的人来说它够用。但问题在于一旦数据量上去MRT的短板就暴露得非常明显。首先是批处理能力弱。MRT虽然提供了命令行版本但参数文件.prm的编写并不直观尤其是涉及多个产品、多个波段组合的时候每个组合都要单独写一个prm文件管理起来很混乱。其次是投影和重采样选项有限MRT内置的投影类型虽然覆盖了常见的几种但如果你想用自定义的Albers等积投影或者特定的地理坐标系配置起来就很别扭。再就是运行稳定性MRT在处理大区域、长时间序列数据时偶尔会卡死或者输出异常值而且错误提示非常模糊排查起来费时费力。我印象最深的一次是处理青藏高原区域连续8年的MOD13A2 NDVI数据MRT跑到第三年的时候突然报了一个Error in processing tile的提示没有任何具体信息。后来逐个文件排查才发现是某一个HDF文件的第13波段存在数据异常MRT直接崩溃。如果换成脚本处理这种问题可以在代码里加异常捕获跳过问题文件继续跑不会影响整体进度。1.2 脚本化处理的三个核心优势转向Python脚本之后最直接的感受是可控性完全不一样了。具体来说脚本化处理带来三个层面的提升第一是流程可复用。一套脚本写好之后换一个研究区、换一个时间段只需要改几个参数就能直接跑不需要重新配置工具。这对于需要反复做实验、调整参数的研究场景来说效率提升是数量级的。第二是异常可追溯。代码里可以加日志记录、异常捕获、中间结果检查每一步都有据可查。哪个文件处理失败了、失败原因是什么、跳过了哪些数据全都清清楚楚。这比MRT那种黑箱式报错要友好太多。第三是与其他工具链无缝衔接。Python处理完的中间结果可以直接喂给ArcGIS做空间分析或者用GDAL、Rasterio做进一步处理整个链条是打通的。而MRT的输出往往还需要额外的格式转换步骤才能进入下一步分析。提示如果你的数据量在10个文件以内MRT确实够用。但只要超过这个量级或者需要反复处理不同区域的数据脚本化就是必然选择。1.3 整体技术路线概览这套流程的整体思路是这样的HDF原始文件 → 波段提取与拼接 → 投影转换与重采样 → 研究区裁剪 → NDVI/ET计算 → 成品图输出。其中前四步用Python脚本完成后两步根据具体需求可以在Python里做也可以导入ArcGIS做后处理和制图。具体用到的工具和库包括环节工具/库作用HDF读取PyHDF / GDAL读取MODIS HDF4格式数据波段提取GDAL提取指定波段并转为GeoTIFF拼接GDAL Warp / Rasterio多景影像拼接投影转换GDAL Warp重投影到目标坐标系裁剪Rasterio / ArcPy按研究区边界裁剪NDVI计算NumPy / Rasterio波段运算ET计算NumPy / Rasterio基于SEBAL或其他模型制图输出ArcGIS / Matplotlib成品图渲染与导出这条路线的好处是每一环都可以独立调试出了问题容易定位。而且Python生态里的这些库都是成熟稳定的社区支持也好遇到问题基本都能找到解决方案。2. 环境搭建与依赖库的取舍2.1 Python环境的选择处理MODIS数据对Python环境的要求其实不算高但有几个库的安装比较讲究。我建议用Anaconda来管理环境因为GDAL、Rasterio这些库在Windows下用pip安装经常出问题而conda源里有预编译好的版本省去很多麻烦。创建一个独立的环境conda create -n modis python3.9 conda activate modis为什么选Python 3.9而不是更新的版本因为GDAL和PyHDF在某些新版本Python上的兼容性还不稳定3.9是目前最稳妥的选择。我试过3.11PyHDF编译报错折腾了半天还是退回3.9。2.2 核心依赖库安装核心依赖库包括以下几个conda install -c conda-forge gdal rasterio pyhdf numpy matplotlib这里重点说一下PyHDF和GDAL的分工。PyHDF专门用来读取HDF4格式的文件MODIS的很多产品比如MOD13、MOD16都是HDF4格式用PyHDF读取最直接。GDAL虽然也能读HDF但对HDF4的支持需要额外编译有时候会出问题。所以我的做法是用PyHDF读数据和元信息用GDAL做投影转换和格式输出。Rasterio是基于GDAL的Python封装API更友好适合做裁剪和波段运算。NumPy不用说了做数组运算的基础。注意如果你在Windows下用conda安装GDAL失败可以试试先安装conda install -c conda-forge gdal3.4指定版本或者用mamba替代conda解决依赖冲突的能力更强。2.3 ArcGIS的角色定位ArcGIS在这套流程里主要承担两个角色后处理和制图输出。后处理包括一些Python脚本不太方便做的操作比如栅格计算器的复杂表达式、邻域分析、重分类等。制图输出则是ArcGIS的强项做出来的图专业、美观适合直接放到论文或报告里。如果你用的是ArcGIS Pro它自带的Python环境是ArcGIS Pro的conda环境可以直接调用ArcPy。但要注意ArcGIS Pro的Python环境和你自己创建的conda环境是分开的如果要在ArcPy里调用GDAL需要在ArcGIS Pro的环境里也安装GDAL。我的建议是分开处理Python脚本用自己的环境跑输出GeoTIFF后再导入ArcGIS做后续处理避免环境冲突。3. HDF文件的批量读取与波段提取3.1 MODIS HDF文件的内部结构MODIS的HDF文件本质上是一个容器里面包含了多个数据集SDSScientific Data Sets。以MOD13A2NDVI产品为例一个HDF文件里包含NDVI归一化植被指数EVI增强植被指数VI_Quality质量波段pixel_reliability像元可靠性sur_refl_b01到sur_refl_b07七个地表反射率波段角度信息太阳天顶角、观测天顶角等每个数据集都有自己的维度和数据类型。NDVI通常是int16类型需要乘以缩放因子scale factor才能得到真实值。MOD13A2的NDVI缩放因子是0.0001也就是说原始值10000对应真实NDVI值1.0。用PyHDF查看文件结构from pyhdf.SD import SD, SDC hdf SD(MOD13A2.A2020001.h26v05.006.2020015000000.hdf, SDC.READ) datasets hdf.datasets() for name, info in datasets.items(): print(f数据集: {name}, 维度: {info[1]}, 类型: {info[3]})这段代码会列出文件里所有数据集的名称、维度和数据类型。第一次处理某个产品的时候建议先跑一下这个确认波段名称和结构因为不同产品的波段命名规则不一样。3.2 批量读取与波段提取脚本批量处理的核心逻辑是遍历文件夹下所有HDF文件 → 逐个读取指定波段 → 输出为GeoTIFF。这里有个关键点MODIS的HDF文件本身带有地理定位信息但PyHDF读取出来的只是纯数组没有地理坐标。所以需要额外读取经纬度信息或者用GDAL来读取带地理信息的波段。我的做法是用GDAL直接读取HDF中的子数据集这样输出的GeoTIFF自带地理坐标import os import glob from osgeo import gdal def extract_band(hdf_path, band_name, output_path): 从HDF文件中提取指定波段并输出为GeoTIFF # 构建子数据集路径 subdataset fHDF4_EOS:EOS_GRID:{hdf_path}:MODIS_Grid_16DAY_1km_VI:{band_name} # 打开子数据集 ds gdal.Open(subdataset) if ds is None: print(f无法打开: {subdataset}) return False # 输出为GeoTIFF gdal.Translate(output_path, ds, formatGTiff) ds None return True # 批量处理 input_dir rD:\MODIS\MOD13A2 output_dir rD:\MODIS\NDVI_tif os.makedirs(output_dir, exist_okTrue) hdf_files glob.glob(os.path.join(input_dir, *.hdf)) for hdf_file in hdf_files: filename os.path.basename(hdf_file).replace(.hdf, .tif) output_path os.path.join(output_dir, filename) extract_band(hdf_file, NDVI, output_path) print(f已处理: {filename})这里的关键是子数据集路径的构建。HDF4_EOS:EOS_GRID是GDAL读取HDF4的固定前缀后面的MODIS_Grid_16DAY_1km_VI是网格名称不同产品的网格名称不一样。MOD13A2是MODIS_Grid_16DAY_1km_VIMOD16A2ET产品是MODIS_Grid_8Day_1km_LSTE或者MODIS_Grid_8Day_1km_ET具体要看产品文档。提示如果不确定网格名称可以用gdalinfo命令查看HDF文件的信息里面会列出所有子数据集的完整路径。3.3 缩放因子与无效值处理提取出来的NDVI数据是int16类型值域通常在-2000到10000之间。要得到真实的NDVI值需要乘以缩放因子0.0001。同时MODIS数据有专门的无效值填充fill value通常是-3000需要处理掉。import numpy as np import rasterio def apply_scale_factor(input_tif, output_tif, scale0.0001, fill_value-3000): 应用缩放因子并处理无效值 with rasterio.open(input_tif) as src: data src.read(1).astype(np.float32) profile src.profile.copy() # 处理无效值 data[data fill_value] np.nan # 应用缩放因子 data data * scale # 更新数据类型 profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_tif, w, **profile) as dst: dst.write(data, 1)这一步看起来简单但不做缩放的话后续所有分析都是错的。我见过不少初学者直接拿原始int16值去做时序分析结果NDVI值域变成-2000到10000完全没法用。所以这个步骤一定要养成习惯。4. 拼接、投影转换与研究区裁剪4.1 多景影像的拼接逻辑MODIS数据是按瓦片tile组织的全球被划分为36×18个瓦片每个瓦片覆盖约1100km×1100km。如果你的研究区跨越多个瓦片就需要先拼接。比如青藏高原区域通常涉及h25v05、h26v05、h25v06、h26v06四个瓦片。拼接用GDAL的Warp工具最方便from osgeo import gdal def mosaic_tiles(input_files, output_path): 拼接多个瓦片 options gdal.WarpOptions( formatGTiff, resampleAlggdal.GRA_NearestNeighbour, srcNodata-3000, dstNodata-3000 ) gdal.Warp(output_path, input_files, **options)这里重采样方法的选择很关键。对于NDVI这种连续值数据理论上用双线性插值GRA_Bilinear更平滑但会引入不属于原始数据的新值。对于分类数据或者质量波段必须用最近邻GRA_NearestNeighbour。我的习惯是NDVI用双线性质量波段用最近邻这样既保证平滑又不破坏分类信息。4.2 投影转换的参数设定MODIS原始数据用的是正弦投影Sinusoidal这种投影在低纬度地区变形较小但在高纬度地区面积变形严重。做区域分析的时候通常需要转换到更适合的投影比如Albers等积投影或者UTM投影。以中国区域为例常用的Albers投影参数是albers_proj projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs def reproject(input_tif, output_tif, dst_proj): 投影转换 options gdal.WarpOptions( formatGTiff, dstSRSdst_proj, resampleAlggdal.GRA_Bilinear, xRes1000, yRes1000, srcNodata-3000, dstNodata-3000 ) gdal.Warp(output_tif, input_tif, **options)分辨率的选择要根据原始数据来定。MOD13A2的空间分辨率是1km所以重采样到1000m是合理的。如果你重采样到500m虽然文件变大了但信息量并没有增加反而引入了插值误差。这一点很多人容易忽略觉得分辨率越高越好其实不是。4.3 按研究区边界精确裁剪裁剪需要两部分数据待裁剪的栅格和研究区边界矢量。用Rasterio的mask功能可以很方便地实现import rasterio from rasterio.mask import mask import geopandas as gpd def clip_raster(input_tif, boundary_shp, output_tif): 按矢量边界裁剪栅格 # 读取边界 gdf gpd.read_file(boundary_shp) geometries gdf.geometry.values with rasterio.open(input_tif) as src: # 裁剪 out_image, out_transform mask(src, geometries, cropTrue) out_meta src.meta.copy() # 更新元数据 out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(output_tif, w, **out_meta) as dst: dst.write(out_image)这里有个常见的坑如果研究区边界是经纬度坐标WGS84而栅格是投影坐标Albers直接裁剪会报错或者结果为空。解决办法是先统一坐标系把边界矢量转换到和栅格一样的投影下gdf gdf.to_crs(src.crs)这一步看起来简单但实际工作中经常忘记导致裁剪结果为空排查半天才发现是坐标系不匹配。注意裁剪的时候cropTrue会把栅格裁剪到边界的最小外接矩形边界外的像元会被设为nodata。如果你需要精确到边界形状还需要额外做掩膜处理。5. NDVI与ET成品的计算与后处理5.1 NDVI的时序合成与最大值合成法NDVI数据通常是16天合成的一年有23期。做年度分析的时候常用的方法是最大值合成法MVCMaximum Value Composite即取一年中每个像元的最大NDVI值这样可以最大程度消除云污染和大气影响。import numpy as np import rasterio import glob def mvc_composite(ndvi_files, output_path): 最大值合成 # 读取所有期数据 arrays [] for f in sorted(ndvi_files): with rasterio.open(f) as src: arrays.append(src.read(1)) profile src.profile.copy() # 堆叠并取最大值 stack np.stack(arrays, axis0) max_ndvi np.nanmax(stack, axis0) # 输出 profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_path, w, **profile) as dst: dst.write(max_ndvi.astype(np.float32), 1)MVC的逻辑很直观一年中NDVI最高的那期大概率是植被生长最旺盛、云污染最少的时候。这个方法虽然简单但在植被动态监测中非常有效也是文献里最常用的方法之一。5.2 ET产品的处理要点MODIS的ET产品MOD16A2和NDVI产品在处理上有几个不同点需要注意第一是数据单位。MOD16A2的ET数据单位是kg/m²/8day也就是每8天的蒸散量。要换算成日均ET需要除以8。如果要换算成mm由于水的密度是1kg/L1kg/m²等于1mm所以数值上是一样的。第二是质量波段。MOD16A2有专门的质量控制波段ET_QC需要根据QC值筛选有效数据。QC值是一个8位整数不同的位代表不同的质量信息。通常的做法是只保留QC值小于某个阈值的像元。第三是无效值处理。MOD16A2的填充值是32767需要单独处理。def process_et(et_tif, qc_tif, output_tif): 处理ET数据筛选有效像元 with rasterio.open(et_tif) as src: et src.read(1).astype(np.float32) profile src.profile.copy() with rasterio.open(qc_tif) as src: qc src.read(1) # 处理填充值 et[et 32767] np.nan # 根据QC筛选保留QC小于64的像元即质量较好的 et[qc 64] np.nan # 转换为日均ET除以8 et et / 8.0 profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_tif, w, **profile) as dst: dst.write(et, 1)5.3 在ArcGIS中做后处理与制图Python处理完的GeoTIFF导入ArcGIS后可以做进一步的后处理和制图。几个常用的操作栅格计算器做归一化、重分类、阈值提取等。比如把NDVI小于0.1的像元设为0非植被大于0.1的保留。符号化NDVI通常用绿-黄-红的渐变色带ET用蓝-绿-黄的色带。ArcGIS的色带库里有现成的也可以自定义。布局出图添加指北针、比例尺、图例、经纬网格导出为高分辨率PNG或PDF。这一步是ArcGIS的强项做出来的图比Matplotlib要专业得多。提示如果要在ArcGIS里批量出图可以用ArcPy写脚本结合地图文档.mxd模板自动替换数据源并导出。这在做多期对比图的时候特别有用。6. 批量处理中的踩坑记录与性能优化6.1 内存溢出与分块处理处理大区域、长时间序列数据的时候最容易遇到的问题就是内存溢出。尤其是做MVC合成的时候如果把23期1km分辨率的全国数据全部读进内存光是一个波段就要占几个GB加上中间变量很容易爆内存。解决办法是分块处理。Rasterio支持按窗口读取数据def mvc_composite_block(ndvi_files, output_path, block_size1024): 分块最大值合成 with rasterio.open(ndvi_files[0]) as src: profile src.profile.copy() height, width src.height, src.width profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_path, w, **profile) as dst: for i in range(0, height, block_size): for j in range(0, width, block_size): # 计算当前块的范围 win_height min(block_size, height - i) win_width min(block_size, width - j) window rasterio.windows.Window(j, i, win_width, win_height) # 读取所有期的当前块 blocks [] for f in ndvi_files: with rasterio.open(f) as src: blocks.append(src.read(1, windowwindow)) # 取最大值 stack np.stack(blocks, axis0) max_block np.nanmax(stack, axis0) # 写入 dst.write(max_block.astype(np.float32), 1, windowwindow)分块处理虽然代码复杂一点但内存占用可以控制在几百MB以内处理全国数据也没问题。块大小的选择要根据你的内存来定一般1024×1024或者2048×2048比较合适。6.2 坐标系不匹配导致的裁剪失败这个坑我在前面提过但值得再强调一次。裁剪失败最常见的原因就是坐标系不匹配。症状是裁剪结果为空或者只有一小块。排查方法是用gdalinfo查看栅格的坐标系用ogrinfo查看矢量的坐标系确认两者是否一致如果不一致用gdf.to_crs()转换矢量或者用gdal.Warp转换栅格。我个人的习惯是统一用投影坐标系因为投影坐标系下的距离和面积计算更准确。6.3 并行处理加速批量任务如果数据量很大单线程处理太慢可以用Python的multiprocessing库做并行。比如批量提取波段的时候每个文件独立处理天然适合并行from multiprocessing import Pool import os def process_single_file(hdf_file): 处理单个文件 filename os.path.basename(hdf_file).replace(.hdf, .tif) output_path os.path.join(output_dir, filename) extract_band(hdf_file, NDVI, output_path) return filename if __name__ __main__: hdf_files glob.glob(os.path.join(input_dir, *.hdf)) with Pool(processes4) as pool: results pool.map(process_single_file, hdf_files) print(f完成 {len(results)} 个文件)进程数根据你的CPU核心数来定一般设为核心数的一半到全部。注意GDAL在多进程下有时候会有问题如果遇到崩溃可以改用ThreadPool或者降低进程数。6.4 常见错误与快速排查表错误现象可能原因解决办法裁剪结果为空坐标系不匹配统一矢量与栅格坐标系NDVI值域异常未应用缩放因子乘以0.0001内存溢出一次性读取全部数据分块处理HDF读取失败网格名称错误用gdalinfo查看子数据集路径投影转换后变形重采样方法不当连续值用双线性分类值用最近邻并行处理崩溃GDAL多进程冲突降低进程数或改用线程这张表是我在实际工作中总结出来的基本上覆盖了80%以上的常见问题。遇到报错的时候先对照这张表排查能省不少时间。7. 从脚本到工作流我的实际使用体会整套流程跑通之后我把它封装成了一个配置文件驱动的工作流。所有的输入路径、输出路径、投影参数、分辨率、缩放因子都写在一个YAML文件里主脚本读取配置后自动执行。这样换研究区的时候只需要改配置文件不用动代码。# config.yaml input_dir: D:/MODIS/MOD13A2 output_dir: D:/MODIS/output product: MOD13A2 band: NDVI scale_factor: 0.0001 fill_value: -3000 target_proj: projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm resolution: 1000 boundary: D:/data/study_area.shp这种配置驱动的方式最大的好处是可追溯。半年后回头看某个结果是怎么生成的翻出配置文件就一目了然不用去回忆当时用了什么参数。另外一个小技巧是中间结果保留。很多人为了省磁盘空间处理完就把中间文件删了。我的建议是至少保留波段提取后的GeoTIFF因为后续如果要调整投影或者裁剪范围可以直接从这一步重新跑不用从头再来。磁盘空间现在很便宜但重新处理的时间成本很高。最后说一个关于ArcGIS和Python配合的心得。我现在的习惯是Python做批量和计算ArcGIS做交互和出图。Python脚本负责把原始数据变成干净的、带地理坐标的GeoTIFF然后导入ArcGIS做符号化、布局、导出。两者各司其职效率最高。如果硬要用ArcPy做所有事情反而会因为ArcGIS的Python环境限制而束手束脚。这套流程我用了三年多从最初的几十个文件到现在每年处理上万景数据稳定性一直很好。核心逻辑没有变过只是根据具体需求做了些微调。如果你也在做类似的批量栅格处理希望这些经验能帮你少走些弯路。
Python批量处理MODIS遥感数据：从HDF到NDVI/ET成品图全流程
遥感数据处理这件事最怕的不是算法难而是重复劳动。我最早接触MODIS数据的时候一个研究区、一年12个月、4种产品光是HDF文件的拼接、投影、裁剪、格式转换就能耗掉整整两天中间还得盯着MRT那个老界面反复点选参数稍不留神就报错重来。后来被逼得没办法硬是把整套流程用Python脚本串了起来配合ArcGIS做后处理和出图现在同样的工作量压缩到半小时以内而且可以批量跑、无人值守。这篇内容就是把这套流程完整拆开讲清楚从HDF原始文件到NDVI、ET成品图每一步为什么这么做、参数怎么定、坑在哪里都会讲到。适合做生态遥感、农业监测、水文分析方向的同学也适合任何需要批量处理栅格数据的GIS从业者参考。1. 为什么放弃MRT转向脚本化处理1.1 MRT的便利与它的天花板MRTMODIS Reprojection Tool在早期确实是处理MODIS数据的标配工具图形界面友好支持批量导入HDF、选择波段、设定投影和输出格式。对于偶尔处理几个文件的人来说它够用。但问题在于一旦数据量上去MRT的短板就暴露得非常明显。首先是批处理能力弱。MRT虽然提供了命令行版本但参数文件.prm的编写并不直观尤其是涉及多个产品、多个波段组合的时候每个组合都要单独写一个prm文件管理起来很混乱。其次是投影和重采样选项有限MRT内置的投影类型虽然覆盖了常见的几种但如果你想用自定义的Albers等积投影或者特定的地理坐标系配置起来就很别扭。再就是运行稳定性MRT在处理大区域、长时间序列数据时偶尔会卡死或者输出异常值而且错误提示非常模糊排查起来费时费力。我印象最深的一次是处理青藏高原区域连续8年的MOD13A2 NDVI数据MRT跑到第三年的时候突然报了一个Error in processing tile的提示没有任何具体信息。后来逐个文件排查才发现是某一个HDF文件的第13波段存在数据异常MRT直接崩溃。如果换成脚本处理这种问题可以在代码里加异常捕获跳过问题文件继续跑不会影响整体进度。1.2 脚本化处理的三个核心优势转向Python脚本之后最直接的感受是可控性完全不一样了。具体来说脚本化处理带来三个层面的提升第一是流程可复用。一套脚本写好之后换一个研究区、换一个时间段只需要改几个参数就能直接跑不需要重新配置工具。这对于需要反复做实验、调整参数的研究场景来说效率提升是数量级的。第二是异常可追溯。代码里可以加日志记录、异常捕获、中间结果检查每一步都有据可查。哪个文件处理失败了、失败原因是什么、跳过了哪些数据全都清清楚楚。这比MRT那种黑箱式报错要友好太多。第三是与其他工具链无缝衔接。Python处理完的中间结果可以直接喂给ArcGIS做空间分析或者用GDAL、Rasterio做进一步处理整个链条是打通的。而MRT的输出往往还需要额外的格式转换步骤才能进入下一步分析。提示如果你的数据量在10个文件以内MRT确实够用。但只要超过这个量级或者需要反复处理不同区域的数据脚本化就是必然选择。1.3 整体技术路线概览这套流程的整体思路是这样的HDF原始文件 → 波段提取与拼接 → 投影转换与重采样 → 研究区裁剪 → NDVI/ET计算 → 成品图输出。其中前四步用Python脚本完成后两步根据具体需求可以在Python里做也可以导入ArcGIS做后处理和制图。具体用到的工具和库包括环节工具/库作用HDF读取PyHDF / GDAL读取MODIS HDF4格式数据波段提取GDAL提取指定波段并转为GeoTIFF拼接GDAL Warp / Rasterio多景影像拼接投影转换GDAL Warp重投影到目标坐标系裁剪Rasterio / ArcPy按研究区边界裁剪NDVI计算NumPy / Rasterio波段运算ET计算NumPy / Rasterio基于SEBAL或其他模型制图输出ArcGIS / Matplotlib成品图渲染与导出这条路线的好处是每一环都可以独立调试出了问题容易定位。而且Python生态里的这些库都是成熟稳定的社区支持也好遇到问题基本都能找到解决方案。2. 环境搭建与依赖库的取舍2.1 Python环境的选择处理MODIS数据对Python环境的要求其实不算高但有几个库的安装比较讲究。我建议用Anaconda来管理环境因为GDAL、Rasterio这些库在Windows下用pip安装经常出问题而conda源里有预编译好的版本省去很多麻烦。创建一个独立的环境conda create -n modis python3.9 conda activate modis为什么选Python 3.9而不是更新的版本因为GDAL和PyHDF在某些新版本Python上的兼容性还不稳定3.9是目前最稳妥的选择。我试过3.11PyHDF编译报错折腾了半天还是退回3.9。2.2 核心依赖库安装核心依赖库包括以下几个conda install -c conda-forge gdal rasterio pyhdf numpy matplotlib这里重点说一下PyHDF和GDAL的分工。PyHDF专门用来读取HDF4格式的文件MODIS的很多产品比如MOD13、MOD16都是HDF4格式用PyHDF读取最直接。GDAL虽然也能读HDF但对HDF4的支持需要额外编译有时候会出问题。所以我的做法是用PyHDF读数据和元信息用GDAL做投影转换和格式输出。Rasterio是基于GDAL的Python封装API更友好适合做裁剪和波段运算。NumPy不用说了做数组运算的基础。注意如果你在Windows下用conda安装GDAL失败可以试试先安装conda install -c conda-forge gdal3.4指定版本或者用mamba替代conda解决依赖冲突的能力更强。2.3 ArcGIS的角色定位ArcGIS在这套流程里主要承担两个角色后处理和制图输出。后处理包括一些Python脚本不太方便做的操作比如栅格计算器的复杂表达式、邻域分析、重分类等。制图输出则是ArcGIS的强项做出来的图专业、美观适合直接放到论文或报告里。如果你用的是ArcGIS Pro它自带的Python环境是ArcGIS Pro的conda环境可以直接调用ArcPy。但要注意ArcGIS Pro的Python环境和你自己创建的conda环境是分开的如果要在ArcPy里调用GDAL需要在ArcGIS Pro的环境里也安装GDAL。我的建议是分开处理Python脚本用自己的环境跑输出GeoTIFF后再导入ArcGIS做后续处理避免环境冲突。3. HDF文件的批量读取与波段提取3.1 MODIS HDF文件的内部结构MODIS的HDF文件本质上是一个容器里面包含了多个数据集SDSScientific Data Sets。以MOD13A2NDVI产品为例一个HDF文件里包含NDVI归一化植被指数EVI增强植被指数VI_Quality质量波段pixel_reliability像元可靠性sur_refl_b01到sur_refl_b07七个地表反射率波段角度信息太阳天顶角、观测天顶角等每个数据集都有自己的维度和数据类型。NDVI通常是int16类型需要乘以缩放因子scale factor才能得到真实值。MOD13A2的NDVI缩放因子是0.0001也就是说原始值10000对应真实NDVI值1.0。用PyHDF查看文件结构from pyhdf.SD import SD, SDC hdf SD(MOD13A2.A2020001.h26v05.006.2020015000000.hdf, SDC.READ) datasets hdf.datasets() for name, info in datasets.items(): print(f数据集: {name}, 维度: {info[1]}, 类型: {info[3]})这段代码会列出文件里所有数据集的名称、维度和数据类型。第一次处理某个产品的时候建议先跑一下这个确认波段名称和结构因为不同产品的波段命名规则不一样。3.2 批量读取与波段提取脚本批量处理的核心逻辑是遍历文件夹下所有HDF文件 → 逐个读取指定波段 → 输出为GeoTIFF。这里有个关键点MODIS的HDF文件本身带有地理定位信息但PyHDF读取出来的只是纯数组没有地理坐标。所以需要额外读取经纬度信息或者用GDAL来读取带地理信息的波段。我的做法是用GDAL直接读取HDF中的子数据集这样输出的GeoTIFF自带地理坐标import os import glob from osgeo import gdal def extract_band(hdf_path, band_name, output_path): 从HDF文件中提取指定波段并输出为GeoTIFF # 构建子数据集路径 subdataset fHDF4_EOS:EOS_GRID:{hdf_path}:MODIS_Grid_16DAY_1km_VI:{band_name} # 打开子数据集 ds gdal.Open(subdataset) if ds is None: print(f无法打开: {subdataset}) return False # 输出为GeoTIFF gdal.Translate(output_path, ds, formatGTiff) ds None return True # 批量处理 input_dir rD:\MODIS\MOD13A2 output_dir rD:\MODIS\NDVI_tif os.makedirs(output_dir, exist_okTrue) hdf_files glob.glob(os.path.join(input_dir, *.hdf)) for hdf_file in hdf_files: filename os.path.basename(hdf_file).replace(.hdf, .tif) output_path os.path.join(output_dir, filename) extract_band(hdf_file, NDVI, output_path) print(f已处理: {filename})这里的关键是子数据集路径的构建。HDF4_EOS:EOS_GRID是GDAL读取HDF4的固定前缀后面的MODIS_Grid_16DAY_1km_VI是网格名称不同产品的网格名称不一样。MOD13A2是MODIS_Grid_16DAY_1km_VIMOD16A2ET产品是MODIS_Grid_8Day_1km_LSTE或者MODIS_Grid_8Day_1km_ET具体要看产品文档。提示如果不确定网格名称可以用gdalinfo命令查看HDF文件的信息里面会列出所有子数据集的完整路径。3.3 缩放因子与无效值处理提取出来的NDVI数据是int16类型值域通常在-2000到10000之间。要得到真实的NDVI值需要乘以缩放因子0.0001。同时MODIS数据有专门的无效值填充fill value通常是-3000需要处理掉。import numpy as np import rasterio def apply_scale_factor(input_tif, output_tif, scale0.0001, fill_value-3000): 应用缩放因子并处理无效值 with rasterio.open(input_tif) as src: data src.read(1).astype(np.float32) profile src.profile.copy() # 处理无效值 data[data fill_value] np.nan # 应用缩放因子 data data * scale # 更新数据类型 profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_tif, w, **profile) as dst: dst.write(data, 1)这一步看起来简单但不做缩放的话后续所有分析都是错的。我见过不少初学者直接拿原始int16值去做时序分析结果NDVI值域变成-2000到10000完全没法用。所以这个步骤一定要养成习惯。4. 拼接、投影转换与研究区裁剪4.1 多景影像的拼接逻辑MODIS数据是按瓦片tile组织的全球被划分为36×18个瓦片每个瓦片覆盖约1100km×1100km。如果你的研究区跨越多个瓦片就需要先拼接。比如青藏高原区域通常涉及h25v05、h26v05、h25v06、h26v06四个瓦片。拼接用GDAL的Warp工具最方便from osgeo import gdal def mosaic_tiles(input_files, output_path): 拼接多个瓦片 options gdal.WarpOptions( formatGTiff, resampleAlggdal.GRA_NearestNeighbour, srcNodata-3000, dstNodata-3000 ) gdal.Warp(output_path, input_files, **options)这里重采样方法的选择很关键。对于NDVI这种连续值数据理论上用双线性插值GRA_Bilinear更平滑但会引入不属于原始数据的新值。对于分类数据或者质量波段必须用最近邻GRA_NearestNeighbour。我的习惯是NDVI用双线性质量波段用最近邻这样既保证平滑又不破坏分类信息。4.2 投影转换的参数设定MODIS原始数据用的是正弦投影Sinusoidal这种投影在低纬度地区变形较小但在高纬度地区面积变形严重。做区域分析的时候通常需要转换到更适合的投影比如Albers等积投影或者UTM投影。以中国区域为例常用的Albers投影参数是albers_proj projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs def reproject(input_tif, output_tif, dst_proj): 投影转换 options gdal.WarpOptions( formatGTiff, dstSRSdst_proj, resampleAlggdal.GRA_Bilinear, xRes1000, yRes1000, srcNodata-3000, dstNodata-3000 ) gdal.Warp(output_tif, input_tif, **options)分辨率的选择要根据原始数据来定。MOD13A2的空间分辨率是1km所以重采样到1000m是合理的。如果你重采样到500m虽然文件变大了但信息量并没有增加反而引入了插值误差。这一点很多人容易忽略觉得分辨率越高越好其实不是。4.3 按研究区边界精确裁剪裁剪需要两部分数据待裁剪的栅格和研究区边界矢量。用Rasterio的mask功能可以很方便地实现import rasterio from rasterio.mask import mask import geopandas as gpd def clip_raster(input_tif, boundary_shp, output_tif): 按矢量边界裁剪栅格 # 读取边界 gdf gpd.read_file(boundary_shp) geometries gdf.geometry.values with rasterio.open(input_tif) as src: # 裁剪 out_image, out_transform mask(src, geometries, cropTrue) out_meta src.meta.copy() # 更新元数据 out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(output_tif, w, **out_meta) as dst: dst.write(out_image)这里有个常见的坑如果研究区边界是经纬度坐标WGS84而栅格是投影坐标Albers直接裁剪会报错或者结果为空。解决办法是先统一坐标系把边界矢量转换到和栅格一样的投影下gdf gdf.to_crs(src.crs)这一步看起来简单但实际工作中经常忘记导致裁剪结果为空排查半天才发现是坐标系不匹配。注意裁剪的时候cropTrue会把栅格裁剪到边界的最小外接矩形边界外的像元会被设为nodata。如果你需要精确到边界形状还需要额外做掩膜处理。5. NDVI与ET成品的计算与后处理5.1 NDVI的时序合成与最大值合成法NDVI数据通常是16天合成的一年有23期。做年度分析的时候常用的方法是最大值合成法MVCMaximum Value Composite即取一年中每个像元的最大NDVI值这样可以最大程度消除云污染和大气影响。import numpy as np import rasterio import glob def mvc_composite(ndvi_files, output_path): 最大值合成 # 读取所有期数据 arrays [] for f in sorted(ndvi_files): with rasterio.open(f) as src: arrays.append(src.read(1)) profile src.profile.copy() # 堆叠并取最大值 stack np.stack(arrays, axis0) max_ndvi np.nanmax(stack, axis0) # 输出 profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_path, w, **profile) as dst: dst.write(max_ndvi.astype(np.float32), 1)MVC的逻辑很直观一年中NDVI最高的那期大概率是植被生长最旺盛、云污染最少的时候。这个方法虽然简单但在植被动态监测中非常有效也是文献里最常用的方法之一。5.2 ET产品的处理要点MODIS的ET产品MOD16A2和NDVI产品在处理上有几个不同点需要注意第一是数据单位。MOD16A2的ET数据单位是kg/m²/8day也就是每8天的蒸散量。要换算成日均ET需要除以8。如果要换算成mm由于水的密度是1kg/L1kg/m²等于1mm所以数值上是一样的。第二是质量波段。MOD16A2有专门的质量控制波段ET_QC需要根据QC值筛选有效数据。QC值是一个8位整数不同的位代表不同的质量信息。通常的做法是只保留QC值小于某个阈值的像元。第三是无效值处理。MOD16A2的填充值是32767需要单独处理。def process_et(et_tif, qc_tif, output_tif): 处理ET数据筛选有效像元 with rasterio.open(et_tif) as src: et src.read(1).astype(np.float32) profile src.profile.copy() with rasterio.open(qc_tif) as src: qc src.read(1) # 处理填充值 et[et 32767] np.nan # 根据QC筛选保留QC小于64的像元即质量较好的 et[qc 64] np.nan # 转换为日均ET除以8 et et / 8.0 profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_tif, w, **profile) as dst: dst.write(et, 1)5.3 在ArcGIS中做后处理与制图Python处理完的GeoTIFF导入ArcGIS后可以做进一步的后处理和制图。几个常用的操作栅格计算器做归一化、重分类、阈值提取等。比如把NDVI小于0.1的像元设为0非植被大于0.1的保留。符号化NDVI通常用绿-黄-红的渐变色带ET用蓝-绿-黄的色带。ArcGIS的色带库里有现成的也可以自定义。布局出图添加指北针、比例尺、图例、经纬网格导出为高分辨率PNG或PDF。这一步是ArcGIS的强项做出来的图比Matplotlib要专业得多。提示如果要在ArcGIS里批量出图可以用ArcPy写脚本结合地图文档.mxd模板自动替换数据源并导出。这在做多期对比图的时候特别有用。6. 批量处理中的踩坑记录与性能优化6.1 内存溢出与分块处理处理大区域、长时间序列数据的时候最容易遇到的问题就是内存溢出。尤其是做MVC合成的时候如果把23期1km分辨率的全国数据全部读进内存光是一个波段就要占几个GB加上中间变量很容易爆内存。解决办法是分块处理。Rasterio支持按窗口读取数据def mvc_composite_block(ndvi_files, output_path, block_size1024): 分块最大值合成 with rasterio.open(ndvi_files[0]) as src: profile src.profile.copy() height, width src.height, src.width profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_path, w, **profile) as dst: for i in range(0, height, block_size): for j in range(0, width, block_size): # 计算当前块的范围 win_height min(block_size, height - i) win_width min(block_size, width - j) window rasterio.windows.Window(j, i, win_width, win_height) # 读取所有期的当前块 blocks [] for f in ndvi_files: with rasterio.open(f) as src: blocks.append(src.read(1, windowwindow)) # 取最大值 stack np.stack(blocks, axis0) max_block np.nanmax(stack, axis0) # 写入 dst.write(max_block.astype(np.float32), 1, windowwindow)分块处理虽然代码复杂一点但内存占用可以控制在几百MB以内处理全国数据也没问题。块大小的选择要根据你的内存来定一般1024×1024或者2048×2048比较合适。6.2 坐标系不匹配导致的裁剪失败这个坑我在前面提过但值得再强调一次。裁剪失败最常见的原因就是坐标系不匹配。症状是裁剪结果为空或者只有一小块。排查方法是用gdalinfo查看栅格的坐标系用ogrinfo查看矢量的坐标系确认两者是否一致如果不一致用gdf.to_crs()转换矢量或者用gdal.Warp转换栅格。我个人的习惯是统一用投影坐标系因为投影坐标系下的距离和面积计算更准确。6.3 并行处理加速批量任务如果数据量很大单线程处理太慢可以用Python的multiprocessing库做并行。比如批量提取波段的时候每个文件独立处理天然适合并行from multiprocessing import Pool import os def process_single_file(hdf_file): 处理单个文件 filename os.path.basename(hdf_file).replace(.hdf, .tif) output_path os.path.join(output_dir, filename) extract_band(hdf_file, NDVI, output_path) return filename if __name__ __main__: hdf_files glob.glob(os.path.join(input_dir, *.hdf)) with Pool(processes4) as pool: results pool.map(process_single_file, hdf_files) print(f完成 {len(results)} 个文件)进程数根据你的CPU核心数来定一般设为核心数的一半到全部。注意GDAL在多进程下有时候会有问题如果遇到崩溃可以改用ThreadPool或者降低进程数。6.4 常见错误与快速排查表错误现象可能原因解决办法裁剪结果为空坐标系不匹配统一矢量与栅格坐标系NDVI值域异常未应用缩放因子乘以0.0001内存溢出一次性读取全部数据分块处理HDF读取失败网格名称错误用gdalinfo查看子数据集路径投影转换后变形重采样方法不当连续值用双线性分类值用最近邻并行处理崩溃GDAL多进程冲突降低进程数或改用线程这张表是我在实际工作中总结出来的基本上覆盖了80%以上的常见问题。遇到报错的时候先对照这张表排查能省不少时间。7. 从脚本到工作流我的实际使用体会整套流程跑通之后我把它封装成了一个配置文件驱动的工作流。所有的输入路径、输出路径、投影参数、分辨率、缩放因子都写在一个YAML文件里主脚本读取配置后自动执行。这样换研究区的时候只需要改配置文件不用动代码。# config.yaml input_dir: D:/MODIS/MOD13A2 output_dir: D:/MODIS/output product: MOD13A2 band: NDVI scale_factor: 0.0001 fill_value: -3000 target_proj: projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm resolution: 1000 boundary: D:/data/study_area.shp这种配置驱动的方式最大的好处是可追溯。半年后回头看某个结果是怎么生成的翻出配置文件就一目了然不用去回忆当时用了什么参数。另外一个小技巧是中间结果保留。很多人为了省磁盘空间处理完就把中间文件删了。我的建议是至少保留波段提取后的GeoTIFF因为后续如果要调整投影或者裁剪范围可以直接从这一步重新跑不用从头再来。磁盘空间现在很便宜但重新处理的时间成本很高。最后说一个关于ArcGIS和Python配合的心得。我现在的习惯是Python做批量和计算ArcGIS做交互和出图。Python脚本负责把原始数据变成干净的、带地理坐标的GeoTIFF然后导入ArcGIS做符号化、布局、导出。两者各司其职效率最高。如果硬要用ArcPy做所有事情反而会因为ArcGIS的Python环境限制而束手束脚。这套流程我用了三年多从最初的几十个文件到现在每年处理上万景数据稳定性一直很好。核心逻辑没有变过只是根据具体需求做了些微调。如果你也在做类似的批量栅格处理希望这些经验能帮你少走些弯路。
代码包深度解析与实践指南)
