当前位置： 首页 > web >正文

MODIS数据下载及处理

web 2025/7/21 19:51:16

时间：2025年5月25日 19:16

地点：丹青楼

任务：从MODIS数据中获取NDVI值

本文参考博客：

MODIS数据的简介和下载（一）——MODIS数据简介-CSDN博客

MODIS数据下载及图像处理教程_modis图像-CSDN博客

第一部分：数据介绍

MODIS数据简介

MODIS是传感器而不是卫星，全名是中分辨率成像光谱仪，搭载在Terra和Aqua卫星上

第一台MODIS仪器于1999年在Terra上发射，其赤道穿越时间为上午10：30;

第二台MODIS仪器于2002年在Aqua平台上发射，其穿越赤道时间为下午1：30

Terra-和Aqua-MODIS仪器每1到2天观察一次整个地球表面，获取36个光谱波段的数据，波长从0.4 μm到14.4 μm不等。双MODIS设计旨在优化无云成像，同时最大限度地减少早晨和下午阳光产生的阴影和眩光的光学效应。

MODIS数据产品具有三种不同的分辨率（250米、500米和1公里），有助于提高我们对陆地、海洋和低层大气中发生的全球环境过程和动态的理解。

MODIS产品类型

0级产品：也称原始数据;

1级产品：指L1A数据，已经被赋予定标参数;

2级产品：经过定标定位后数据，本系统产品是国际标准的EOS-HDF格式。包含所有波段数据，是应用比较广泛的一类数据。;

3级产品：在1B数据的基础上，对由遥感器成像过程产生的边缘畸变(Bowtie效应)进行校正，产生L3级产品;

4级产品：由参数文件提供的参数，对图像进行几何纠正，辐射校正，使图像的每一点都有精确的地理编码、反射率和辐射率。L4级产品的MODIS图像进行不同时相的匹配时，误差小于1个像元。该级产品是应用级产品不可缺少的基础;

5级及以上产品：根据各种应用模型开发L5级产品。

MODIS数据36个波段的特征参数与主要应用方向


MODIS数据的典型应用：

（1）植被动态监测以及植被生理生态的多种产品遥感反演
典型的是NDVI合成产品、NPP（净初级生产力）、LAI（叶面积指数）、植被动态变化。这一部分内容主要应用生态环境监测中，尤其是生态学相关研究。

（2）地表温度反演以及相关产品遥感反演
运用MODIS的热红外波段，通过劈窗算法（类似AVHRR）可以反演地表温度。包括温度异常和林火产品。在林学上有广泛的应用。

（3）光学气溶胶厚度反演
这部分是近些年来关注的重点，由于PM2.5或者说雾霾的造成的环境问题日益严重，如何从遥感监测PM2.5是一个研究热点，现在比较普遍的使用MODIS产品来进行光学气溶胶（AOD)反演，然后反演PM2.5。

（4）其他应用
MODIS的应用还有非常多，比如像海表温度反演——事实上这方面的温度反演精度要高于地表温度反演，主要是海水从性质上说属于近似黑体；叶绿素浓度反演；离水辐射；冰雪覆盖监测（以NDSI为例）等。

命名规则

MODIS命名规则如下（摘自CSDN博客）：
MOD04 是产品名称，表示MODIS气溶胶产品。
L2 表示 产品级别，Level2。
A2005224 表示产品时间2005年第224天（以每年1月1日为第一天）。
0205 表示卫星过境时间，换算成北京时间要加8小时。
005表示产品版本，Version005，之前是v004，相比之前版本有很多改进。
2006225195920 表示的是产品处理时间。

第二部分：数据下载

下载地址：在谷歌上搜索MODIS，进入MODIS Web

可以查看到以下信息

需要注意的是：MODIS 植被指数产品只有16天合成和月合成产品

如何选择自己使用哪一种类型的数据?

（1）考虑过境时间，Terra过境时间是上午10：30，Aqua过境时间为下午1：30

（2）研究尺度 250米，500米，1000米

（3）时间分辨率：16天 或 一个月

我研究的是土壤水分监测，因此选择 MOD13Q1

之后进入网址Earthdata Search | Earthdata Search下载

下载流程不再赘述，可见参考博客（我是一个一个下载的，比较笨的方法，批量下载方法没有研究出来）

第三部分：文件格式转换及NDVI参数提取

第一步：转换为tiff文件

# -*- coding: utf-8 -*-
import os
import re
from osgeo import gdaldef mod13q1_to_tif(hdf_path, output_dir):"""转换MOD13Q1 HDF文件为NDVI TIFF"""# 自动创建输出目录os.makedirs(output_dir, exist_ok=True)# 打开HDF文件hdf_ds = gdal.Open(hdf_path, gdal.GA_ReadOnly)# 定位NDVI子数据集（根据MOD13Q1规范）ndvi_sds = Nonefor sds in hdf_ds.GetSubDatasets():sds_name, sds_desc = sds# 匹配250m NDVI数据集（MOD13Q1特定标识）if re.search(r'250m 16 days NDVI', sds_desc):ndvi_sds = sds_namebreakif not ndvi_sds:raise ValueError(f"未找到NDVI子数据集：{hdf_path}")# 构建输出文件名（保留原始命名结构）hdf_basename = os.path.basename(hdf_path)tif_name = hdf_basename.replace('.hdf', '_NDVI.tif')output_path = os.path.join(output_dir, tif_name)# 执行转换sds_ds = gdal.Open(ndvi_sds)driver = gdal.GetDriverByName('GTiff')out_ds = driver.CreateCopy(output_path, sds_ds, 0)# 强制写入元数据out_ds.SetProjection(sds_ds.GetProjection())out_ds.SetGeoTransform(sds_ds.GetGeoTransform())out_ds.FlushCache()# 关闭数据集out_ds = Nonesds_ds = Nonehdf_ds = Nonereturn output_pathif __name__ == "__main__":# 配置路径（示例路径，根据实际情况修改）input_dir = r"E:\MODIS"  # HDF文件存放目录output_dir = r"E:\output_modis"  # 输出目录# 遍历处理所有HDF文件processed_count = 0for filename in os.listdir(input_dir):if filename.lower().endswith(".hdf"):hdf_path = os.path.join(input_dir, filename)try:output_path = mod13q1_to_tif(hdf_path, output_dir)print(f"成功转换：{filename} -> {os.path.basename(output_path)}")processed_count += 1except Exception as e:print(f"转换失败 {filename}: {str(e)}")print(f"\n处理完成！共转换 {processed_count} 个文件")print(f"输出目录：{output_dir}")

第二步：拼接tiff文件

# -*- coding: utf-8 -*-
import os
import re
from osgeo import gdal
from collections import defaultdictdef mosaic_mod13q1(input_dir, output_path):"""拼接MOD13Q1 NDVI数据"""# 按轨道位置分组文件tile_groups = defaultdict(list)pattern = re.compile(r"h(\d{2})v(\d{2})")  # 匹配轨道号hXXvXX# 遍历所有TIFF文件for fname in os.listdir(input_dir):if fname.endswith("_NDVI.tif"):match = pattern.search(fname)if match:tile_groups[(match.group(1), match.group(2))].append(os.path.join(input_dir, fname))# 创建虚拟栅格vrt_path = os.path.splitext(output_path)[0] + ".vrt"vrt_options = gdal.BuildVRTOptions(resampleAlg='near',separate=False,addAlpha=False)# 按轨道位置排序文件（西南到东北）sorted_files = []for h in sorted(set(h for (h, v) in tile_groups.keys())):for v in sorted(set(v for (h, v) in tile_groups.keys())):if (h, v) in tile_groups:sorted_files.extend(tile_groups[(h, v)])# 生成VRTvrt = gdal.BuildVRT(vrt_path, sorted_files, options=vrt_options)vrt = None  # 关闭VRT# 执行拼接warp_options = gdal.WarpOptions(format='GTiff',creationOptions=['TILED=YES','COMPRESS=DEFLATE','BIGTIFF=IF_SAFER'],srcNodata=-3000,  # MOD13Q1 NDVI无效值dstNodata=-3000,multithread=True,warpMemoryLimit=2048)mosaic = gdal.Warp(output_path, vrt_path, options=warp_options)mosaic.FlushCache()mosaic = None# 清理临时文件os.remove(vrt_path)return output_pathif __name__ == "__main__":# 配置路径input_dir = r"E:\output_modis"  # 转换后的TIFF目录output_path = r"E:\output_modis1\MOD13Q1_Global_NDVI.tif"  # 输出路径# 创建输出目录os.makedirs(os.path.dirname(output_path), exist_ok=True)# 执行拼接print("开始拼接...")result = mosaic_mod13q1(input_dir, output_path)# 验证结果ds = gdal.Open(result)print("\n拼接结果信息：")print(f"文件大小：{ds.RasterXSize}x{ds.RasterYSize}")print(f"分辨率：{ds.GetGeoTransform()[1]} 米")print(f"投影：{ds.GetProjection()[:50]}...")ds = Noneprint(f"\n拼接完成！输出文件：{output_path}")

代码是参阅了别的大佬的博客，自己用deepseek修改之后得到的，自己的能力还是比较有限
最终的结果：（数据没有下完，所以会有一定的缺失）

http://www.xdnf.cn/news/9664.html

相关文章：

  • 电商平台 API、数据抓取与爬虫技术的区别及优势分析
  • linux目录
  • CTFSHOW-WEB-36D杯
  • Unity数字人开发笔记——人物模型
  • 【Redis】热点key问题，的原因和处理，一致性哈希，删除大key的方法
  • 【C语言】深入理解C语言中的自定义数据类型：struct、union与enum
  • 大话软工笔记—基本概念
  • 三视图重建 笔记
  • python入门day02
  • 制导与导航总述、分类介绍、MATLABdemo
  • PROFIBUS转PROFINET网关：饲料行业的通信桥梁
  • LeetCode 543 二叉树的直径
  • 使用Miniconda管理Python环境
  • MS3494模拟矩阵开关
  • transformer-PositionalEncoding (对数空间计算实现)
  • 行业案例 | OPPO借助Azure AI Speech国际服务实现音频文件智能转录
  • 基于MATLAB的二维圆形随机骨料生成程序
  • APL Photonics封面成果：KAUST用五边形激光腔刷新物理随机数生成极限——800Gb/s！
  • Selenium 测试框架 - JavaScript
  • Xamarin入门笔记（Xamarin已经被MAUI取代）
  • 利益相关者意见分歧，如何决策
  • 在线临床指标分类信息表转甜甜圈矩阵图
  • 将git最后一次提交把涉及到的文件按原来目录结构提取出来
  • LLM中的Loss与Logits详解
  • 【leetcode】206. 反转链表
  • Linux Shellcode开发（Stager Reverse Shell）
  • 简述MySQL优化锁方面你有什么建议？
  • 彰显国产力量|暴雨亮相2025 C3安全峰会
  • Guava限频器RateLimiter的使用示例
  • STM32学习第一课--工程建立（云端备份与自我复盘）