modis数据下载测试处理

怪我平庸至极—xyw

已于 2025-01-17 15:36:58 修改

阅读量166

点赞数 1

文章标签：前端 javascript linux

于 2025-01-13 10:37:51 首次发布

本文链接：https://blog.youkuaiyun.com/weixin_65080186/article/details/145109103

版权

% Earthdata 登录信息
% Earthdata 登录信息（使用 Token 代替密码认证） 
token = 'eyJ0eXAiOiJKV1QiLCJvcmlnaW4iOiJFYXJ0aGRhdGEgTG9naW4iLCJzaWciOiJlZGxqd3RwdWJrZXlfb3BzIiwiYWxnIjoiUlMyNTYifQ.eyJ0eXBlIjoiVXNlciIsInVpZCI6IjEzNDM1NTM1MyIsImV4cCI6MTczODA1ODg5NSwiaWF0IjoxNzMyODc0ODk1LCJpc3MiOiJodHRwczovL3Vycy5lYXJ0aGRhdGEubmFzYS5nb3YiLCJpZGVudGl0eV9wcm92aWRlciI6ImVkbF9vcHMifbGV2ZWwiPxffEvybVNz_w83pQfg3520JrXIyJti86xbZHhb-bd2lizffWRbFp3wZCYk6WFvVl4won5f9DROFvik9pKwngFNQARZ0GyoOss1v_2ji6IpHMyER5nq_mH0xFKNhJxxzBg1uE0jBObffFzTrZOlKTblxzi8No4tutoqZjQXhP7aK8lHKUsQNgRG1w'; %
% 打开文件
filename1 = 'F:\modis_data\mod_1\mod2024.txt';  % 替换为你的文件名
fid = fopen(filename1, 'r');
% 检查文件是否成功打开
if fid == -1
    error('文件无法打开');
end

lines = {};  % 创建一个空的 cell 数组
line_num = 1;
while ~feof(fid)
    line = fgetl(fid);
    lines{line_num} = line;  % 将每一行内容保存到 cell 数组中
    line_num = line_num + 1;
end
temp=[];
count=0;
for i=1:size(lines,2)
    i
    filename = lines{1,i};
    % 执行下载
    try
        options = weboptions('HeaderFields', {'Authorization', ['Bearer ' token]});
        outfilename = filename(64:108);
        websave(outfilename, filename, options);
%         disp(['Downloaded: ', outfilename]);
    catch
        count=count+1;
        temp(count,1)=i;
        disp(['Failed to download: ', filename]);
    end
end

https://n5eil01u.ecs.nsidc.org/DP4/MOSA/MYD10A1.061/2018.01.31/MYD10A1.A2018031.h08v09.061.2021314061615.hdf

import requests

# Earthdata 登录信息（Token 代替密码认证）
token = 'eyJ0eXAiOiJKV1QiLCJvcmlnaW4iOiJFYXJ0aGRhdGEgTG9naW4iLCJzaWciOiJlZGxqd3RwdWJrZXlfb3BzIiwiYWxnIjoiUlMyNTYifQ.eyJ0eXBlIjoiVXNlciIsInVpZCI6IjEzNDM1NTM1MyIsImV4cCI6MTczODA1ODg5NSwiaWF0IjoxNzMyODc0ODk1LCUYPxffEvybVNz_w83pQfg3520JrXIyJti86xbZHhb-bd2lizffWRbFp3wZCYk6WFvVl4won5f9DROFvik9pKwngFwNQARZ0GyoOss1v_2ji6IpHMyER5nq_mH0xFKNhJxxzBg1uE0jBObffFzTrZOlKTblxzi8No4tutoqZjQXhP7aK8lHKUsQNgRG1w'

# 文件路径
filename1 = 'F:/modis_data/mod_1/mod2024.txt'  # 替换为你的文件名

# 打开文件
try:
    with open(filename1, 'r') as file:
        lines = file.readlines()
except FileNotFoundError:
    raise Exception("文件无法打开")

# 处理每一行
temp = []
count = 0

for i, line in enumerate(lines):
    line = line.strip()  # 移除行尾的换行符
    filename = line
    try:
        # 执行下载
        headers = {'Authorization': f'Bearer {token}'}
        outfilename = filename[63:109]  # 从 64 到 109 索引位置（MATLAB 中是 64 到 108）
        
        response = requests.get(filename, headers=headers, stream=True)
        response.raise_for_status()  # 如果请求失败，会抛出异常
        
        # 保存文件
        with open(outfilename, 'wb') as out_file:
            for chunk in response.iter_content(chunk_size=8192):
                out_file.write(chunk)
        
        print(f"Downloaded: {outfilename}")
    except requests.exceptions.RequestException:
        count += 1
        temp.append(i + 1)  # 保存失败的行数（1-based index）
        print(f"Failed to download: {filename}")

# 打印下载失败的文件行号
print(f"Total failures: {count}")
if count > 0:
    print("Failed to download from lines:", temp)

import os
import shutil

# 源文件夹路径
source_folder = '/path/to/source/folder'

# 目标文件夹路径
mod2017_folder = '/path/to/destination/mod2017'
myd2017_folder = '/path/to/destination/myd2017'
mod2018_folder = '/path/to/destination/mod2018'
myd2018_folder = '/path/to/destination/myd2018'

# 确保目标文件夹存在
os.makedirs(mod2017_folder, exist_ok=True)
os.makedirs(myd2017_folder, exist_ok=True)
os.makedirs(mod2018_folder, exist_ok=True)
os.makedirs(myd2018_folder, exist_ok=True)

# 遍历源文件夹中的文件
for filename in os.listdir(source_folder):
    # 构造文件的完整路径
    file_path = os.path.join(source_folder, filename)
    
    # 跳过目录，确保只处理文件
    if os.path.isdir(file_path):
        continue

    # 根据文件名将文件分类到不同的文件夹
    if filename.startswith('MOD10A1.A2017'):
        shutil.move(file_path, os.path.join(mod2017_folder, filename))
    elif filename.startswith('MYD10A1.A2017'):
        shutil.move(file_path, os.path.join(myd2017_folder, filename))
    elif filename.startswith('MOD10A1.A2018'):
        shutil.move(file_path, os.path.join(mod2018_folder, filename))
    elif filename.startswith('MYD10A1.A2018'):
        shutil.move(file_path, os.path.join(myd2018_folder, filename))

print("文件已分类并移动到相应的文件夹。")

import os
from pymodis import convertmodis_gdal
from datetime import datetime, timedelta

# 设置输入和输出目录
input_dir = 'D:/MODIS_DATA'  # 替换为你的输入数据目录
output_dir = 'D:/MODIS_OUTPUT'  # 替换为你的输出目录

# 设置日期范围（使用开始和结束日期）
start_date = datetime(2022, 1, 1)  # 起始日期
end_date = datetime(2022, 1, 16)  # 截止日期

# 设置初始变量
current_day = start_date
day1 = 1  # 初始天数
epsg = 4326  # 使用 WGS84 坐标系
resampling = 'NEAREST_NEIGHBOR'  # 设置重采样方法

# 循环处理每一天的数据
while current_day <= end_date:
    day_str = current_day.strftime('%Y%j')  # 获取当前日期的 Julian 日期格式
    print(f"正在处理 {day_str} 的数据...")

    # 获取当前日期下的输入文件
    input_files = [os.path.join(input_dir, f) for f in os.listdir(input_dir) if f.startswith(f"MYD10A1.A{day_str}") and f.endswith(".hdf")]

    if not input_files:
        print(f"没有找到 {day_str} 的数据文件，跳过")
        current_day += timedelta(days=1)
        continue

    # 设置输出文件前缀
    prefix = os.path.join(output_dir, f"MOSAIC_{day_str}")  # 输出文件的前缀

    # 创建 convertModisGDAL 对象
    converter = convertmodis_gdal.convertModisGDAL(
        hdfname=input_files,
        prefix=prefix,
        subset='NDSI_Snow_Cover',  # 提取 NDSI_Snow_Cover 波段
        res=None,  # 保持原始分辨率
        outformat='GTiff',  # 输出格式为 GeoTIFF
        epsg=epsg,  # 使用 WGS84 坐标系
    )

    # 执行数据转换（拼接）
    converter.run()
    print(f"{day_str} 拼接完成！")

    # 日期增加 1 天
    current_day += timedelta(days=1)



input_files = [os.path.join(input_dir, f) for f in os.listdir(input_dir) if f.startswith(f"MYD10A1.A{day_str}") and f.endswith(".hdf")]

from pymodis import convertmodis_gdal
import os

# 设置数据目录和输出目录
data_dir = r'D:\ModisTools\MRT\data'  # 原始MODIS数据目录
output_dir = r'D:\ModisTools\MRT\result'  # 输出拼接后的文件目录

# 获取所有符合条件的 HDF 文件
files_by_date = defaultdict(list)
for filename in os.listdir(data_dir):
    if filename.endswith('.hdf') and 'MYD10A1' in filename:
        # 提取文件名中的日期部分（假设日期部分是 'A2018031'，即文件名中的第 7 到第 14 个字符）
        date = filename[6:14]  # 例如：从 'MYD10A1.A2018031' 中提取 'A2018031'
        files_by_date[date].append(os.path.join(data_dir, filename))

# 使用 convertmodis_gdal 拼接相同日期的文件
for date, files in files_by_date.items():
    if len(files) > 1:  # 如果有多个文件，则进行拼接
        # 输出文件路径
        output_file = os.path.join(output_dir, f'mosaic_output_{date}.hdf')

        # 使用 convertmodis_gdal 拼接数据
        try:
            convertmodis_gdal.convert(input_files=files, output_file=output_file, output_format='HDF4Image')

            print(f'拼接完成：{output_file}')
        except Exception as e:
            print(f'处理 {date} 时出现错误: {e}')


[root@himac appfolder_2]# python Stitching.py
Traceback (most recent call last):
  File "Stitching.py", line 2, in <module>
    from pymodis import modis
ImportError: cannot import name 'modis' from 'pymodis' (/root/anaconda3/envs/blockfl/lib/python3.8/site-packages/pymodis/__init__.py)

[root@himac appfolder_2]# python hp.py
['__builtins__', '__cached__', '__doc__', '__file__', '__loader__', '__name__', '__package__', '__path__', '__spec__', '__version__', '__warningregistry__', 'convertmodis', 'convertmodis_gdal', 'downmodis', 'optparse_required', 'parsemodis', 'print_function', 'productmodis', 'qualitymodis', 'warnings']

import os
import glob
from collections import defaultdict
from datetime import datetime

# 定义文件路径和输出路径
input_dir = "/home/data8/platformtest/backend/appfolder/appfolder_2/datadownload/mod2017"
output_path = "/home/data8/platformtest/backend/appfolder/appfolder_2/datadownload/stitchmod2017"

# 定义重投影参数 (假设已经设置好)
reproject_options = (
    (0, 0, 1, 1),  # 示例的仿射变换矩阵 (gt_out)
    (2400, 2400),  # 图像大小 (xSize, ySize)
    "EPSG:4326"     # 目标坐标系
)

# 每一天的文件命名模式
satellite_bands = [
    "h22v04.061", "h23v04.061", "h24v04.061", "h25v04.061",
    "h23v05.061", "h24v05.061", "h25v05.061", "h26v05.061", "h27v05.061",
    "h23v06.061", "h24v06.061", "h25v06.061", "h26v06.061", "h27v06.061"
]

# 获取所有文件列表
file_groups = defaultdict(list)

# 获取所有文件，并按日期分组
for band in satellite_bands:
    # 根据卫星带和日期生成文件模式
    file_pattern = os.path.join(input_dir, f"MOD10A1.A*.{band}*.hdf")
    all_files = glob.glob(file_pattern)
    
    for file in all_files:
        # 提取文件名中的日期 (A2022001 格式)
        basename = os.path.basename(file)
        date_str = basename.split('.')[1][1:]  # 提取如 "2022001" 这部分
        date = datetime.strptime(date_str, "%Y%j").strftime("%Y-%m-%d")  # 转换为 YYYY-MM-DD 格式
        
        # 将文件按日期分组
        file_groups[date].append(file)

# 为每个日期调用 stitch_and_reproject 函数
for date, files in file_groups.items():
    product = "MODIS"
    
    # 调用 stitch_and_reproject 函数
    try:
        output_file = stitch_and_reproject(
            file_list=files,
            output_path=output_path,
            product=product,
            date=date,
            reproject_options=reproject_options
        )
        print(f"Output file for {date} saved at: {output_file}")
    except Exception as e:
        print(f"Error processing {date}: {e}")

def calc_longtime_snow_mask(dep_list: list, dem_mask = GDAL_MASKS.dem_mask) -> np.ndarray:
    '''
    用于计算长期积雪掩模（snow mask），根据给定的多个图像数据，结合 DEM（数字高程模型）掩码和设定的阈值规则，生成一个最终的掩模图像。这个掩模图像会标记出积雪、土地和云层等区域。
    '''
    if dep_list is None or len(dep_list) <= 0:
        return None
    days = len(dep_list)
    res = np.zeros(dem_mask.shape) + CLOUD_CODE
    snow_cnt = np.zeros(dem_mask.shape)
    land_cnt = np.zeros(dem_mask.shape)
    cloud_cnt = np.zeros(dem_mask.shape)
    for img in dep_list:
        if img is None:
            continue
        snow_cnt += img == SNOW_CODE
        land_cnt += img == LAND_CODE
        cloud_cnt += img == CLOUD_CODE
    res[dem_mask >= LONGTIME_ELEVATION_DECISIONS[0]] = SNOW_CODE
    res[(snow_cnt + cloud_cnt >= days * LONGTIME_SNOW_CONTROL_RATIO) & (dem_mask > LONGTIME_ELEVATION_DECISIONS[1]) & (dem_mask < LONGTIME_ELEVATION_DECISIONS[0])] = SNOW_CODE
    res[(cloud_cnt <= days * LONGTIME_LAND_CONTROL_RATIO) & (land_cnt + cloud_cnt == days) & (dem_mask <= LONGTIME_ELEVATION_DECISIONS[1])] = LAND_CODE

    return res