GLASS产品下载后的数据格式是hdf格式,第一步需要将其转化为tif格式,不写代码的话一般是MRT工具批量处理比较方便,详情请看本人之前的博客,这里本人使用python脚本对其处理。

GLASS产品主要是两种产品处理得到的,分别是MODIS和AVHRR。

针对GLASS-MODIS数据集(以ET为例)批量处理代码如下,直接改路径即可,需要修改的地方我都备注了。

# -*- coding: utf-8 -*-
import arcpy
import os
import shutil
folder_path = r'G:\data_ET\GLASS_MODIS'      # 第一层文件夹(最原始文件夹)#改
os.chdir(folder_path)               # 进入第一层文件夹
current_dir = os.getcwd()           # 获取当前目录
# 获取第一层文件夹中的所有文件和文件夹的名称
files_and_folders = os.listdir(folder_path)
print(files_and_folders)
# 遍历第一层文件夹中的所有文件和文件夹
for name in files_and_folders:
    # 构建完整的文件路径
    full_path = os.path.join(folder_path, name)
    # 判断是否为文件夹
    if os.path.isdir(full_path):
        print(name, 'isdir')
    else:
        print(name, 'isfile')
    os.chdir(folder_path+"\\"+''.join(name))  # 进入第二层文件夹
    subcurrent_dir = os.getcwd()         # 获取当前目录
    print(subcurrent_dir)
    subfiles_and_folders = os.listdir(subcurrent_dir)  # 列出第二层文件夹所有内容
    print(subfiles_and_folders)
    for subname in subfiles_and_folders:     # 遍历第二层文件夹中的所有文件和文件夹
        subfull_path = os.path.join(subcurrent_dir, subname)
        if os.path.isdir(subfull_path):
            print(subfull_path, 'isdir')
        else:
            print(subfull_path, 'isfile')
        os.chdir(subfull_path)
        subsubcurrent_dir = os.getcwd()
        print(subsubcurrent_dir)
        dir_path = os.path.join(subsubcurrent_dir, 'my_hdf_to_tif\\')  # 创造存放tif的数据文件夹 #改
        os.mkdir(dir_path)  # 创建新的空的文件夹  #改
        #os.rmdir(dir_path)  # 删除空文件夹   #改
        #   shutil.rmtree(dir_path) #删除包含内容的非空文件夹
        subsubfiles_and_folders = os.listdir(subsubcurrent_dir)
        print(subsubfiles_and_folders)
        for subsubname in subsubfiles_and_folders:   # 遍历第三层文件夹中的所有文件和文件夹
            subsubfull_path = os.path.join(subsubcurrent_dir, subsubname)  #hdf文件和my_hdf_to_tif文件夹
            if os.path.isfile(subsubfull_path):
                print(subsubfull_path, 'isfile')
                subsubname2 = subsubname.rstrip(".hdf")
                subsubname2 = subsubname2 + ".tif"
                subsubfull_path2 = os.path.join(dir_path, subsubname2)  # tif文件
                print(subsubfull_path2)
                arcpy.ExtractSubDataset_management(subsubfull_path, subsubfull_path2)
            else:
                print(subsubfull_path, 'isdir')

运行效果如下,会生成my_hdf_to_tif文件夹,里面存储着生成的tif文件。

第二步,hdf转tif后,需要把研究区不同切片tif文件拼接起来。

# -*- coding: utf-8 -*-
import arcpy
import os
import glob
import shutil
import re
folder_path = r'G:\data_ET\GLASS_MODIS'      # 第一层文件夹(最原始文件夹)#改
os.chdir(folder_path)               # 进入第一层文件夹
current_dir = os.getcwd()           # 获取当前目录
# 获取第一层文件夹中的所有文件和文件夹的名称
files_and_folders = os.listdir(folder_path)
print(files_and_folders)
# 遍历第一层文件夹中的所有文件和文件夹
for name in files_and_folders:
    # 构建完整的文件路径
    full_path = os.path.join(folder_path, name)
    # 判断是否为文件夹
    if os.path.isdir(full_path):
        print(name, 'isdir')
    else:
        print(name, 'isfile')
    os.chdir(folder_path+"\\"+''.join(name))  # 进入第二层文件夹
    subcurrent_dir = os.getcwd()         # 获取当前目录
    print(subcurrent_dir)
    subfiles_and_folders = os.listdir(subcurrent_dir)  # 列出第二层文件夹所有内容
    print(subfiles_and_folders)
    for subname in subfiles_and_folders:     # 遍历第二层文件夹中的所有文件和文件夹
        subfull_path = os.path.join(subcurrent_dir, subname)
        if os.path.isdir(subfull_path):
            print(subfull_path, 'isdir')
        else:
            print(subfull_path, 'isfile')
        os.chdir(subfull_path)
        subsubcurrent_dir = os.getcwd()
        print(subsubcurrent_dir)
        output = os.path.join(subsubcurrent_dir, 'my_mosaic\\')  #改
        os.mkdir(output)  # 创建新的空的文件夹  #改
        # os.rmdir(output)  # 删除空文件夹   #改
        #   shutil.rmtree(output) #删除包含内容的非空文件夹
        dir_path = os.path.join(subsubcurrent_dir, 'my_hdf_to_tif\\') #改
        os.chdir(dir_path)
        subsubsubcurrent_dir = os.getcwd()
        print(subsubsubcurrent_dir)
        subsubsubfiles_and_folders = os.listdir(subsubsubcurrent_dir)
        print(subsubsubfiles_and_folders)
        list = []
        for subsubsubname in glob.glob(os.path.join(subsubsubcurrent_dir, '*.tif')):
            print(subsubsubname)
            list.append(subsubsubname)
        base = list[0]
        pattern = r'[.][h](.*?)[v](.*?)[.]'
        rep = '.'
        output_file = re.sub(pattern, rep, base)
        print(output_file)
        [output_file_dir, output_file_name]=os.path.split(output_file)
        print(output_file_name)
        list = ";".join(list)
        print(list)
        cellwidth = arcpy.Describe(base).meanCellWidth  # 像元宽度
        bandcount = arcpy.Describe(base).bandCount  # 获取bandCount
        Coordinate_System = arcpy.Describe(base).SpatialReference  #坐标系统
        arcpy.MosaicToNewRaster_management(list, output, output_file_name , Coordinate_System,
                                       "32_BIT_FLOAT",
                                     cellwidth, bandcount, "MAXIMUM", "FIRST")

效果如下,将生成my_mosaic文件夹:

第三步,重投影,例如投影到WGS1984

# -*- coding: utf-8 -*-
import arcpy
import os
import glob
import shutil
import re
folder_path = r'G:\data_ET\GLASS_MODIS'      # 第一层文件夹(最原始文件夹)
os.chdir(folder_path)               # 进入第一层文件夹
current_dir = os.getcwd()           # 获取当前目录
# 获取第一层文件夹中的所有文件和文件夹的名称
files_and_folders = os.listdir(folder_path)
print(files_and_folders)
# 遍历第一层文件夹中的所有文件和文件夹
for name in files_and_folders:
    # 构建完整的文件路径
    full_path = os.path.join(folder_path, name)
    # 判断是否为文件夹
    if os.path.isdir(full_path):
        print(name, 'isdir')
    else:
        print(name, 'isfile')
    os.chdir(folder_path+"\\"+''.join(name))  # 进入第二层文件夹
    subcurrent_dir = os.getcwd()         # 获取当前目录
    print(subcurrent_dir)
    subfiles_and_folders = os.listdir(subcurrent_dir)  # 列出第二层文件夹所有内容
    print(subfiles_and_folders)
    for subname in subfiles_and_folders:     # 遍历第二层文件夹中的所有文件和文件夹
        subfull_path = os.path.join(subcurrent_dir, subname)
        if os.path.isdir(subfull_path):
            print(subfull_path, 'isdir')
        else:
            print(subfull_path, 'isfile')
        os.chdir(subfull_path)
        subsubcurrent_dir = os.getcwd()
        print(subsubcurrent_dir)
        output = os.path.join(subsubcurrent_dir, 'my_project_raster\\')
        os.mkdir(output)  # 创建新的空的文件夹  #改
        # os.rmdir(output)  # 删除空文件夹   #改
        #   shutil.rmtree(output) #删除包含内容的非空文件夹
        dir_path = os.path.join(subsubcurrent_dir, 'my_mosaic\\')
        os.chdir(dir_path)
        subsubsubcurrent_dir = os.getcwd()
        print(subsubsubcurrent_dir)
        subsubsubfiles_and_folders = os.listdir(subsubsubcurrent_dir)
        print(subsubsubfiles_and_folders)
        list = []
        for subsubsubname in glob.glob(os.path.join(subsubsubcurrent_dir, '*.tif')):
            print(subsubsubname)
            list.append(subsubsubname)
        print(list)
        input_file = ' '.join(list)
        print(input_file)
        [file_dir, file_name] = os.path.split(input_file)
        output_file = output + file_name
        print(output_file)
        Coordinate_System = arcpy.SpatialReference(4326)  # WGS-1984
        # Coordinate_System = 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],\
        #            PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]];\
        #            -400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119522E-09;\
        #           0.001;0.001;IsHighPrecision'
        arcpy.ProjectRaster_management(input_file, output_file, Coordinate_System, \
                                       "BILINEAR", "", "", "#", "#")

效果如下,将生成my_project_ratser文件夹:

Logo

DAMO开发者矩阵,由阿里巴巴达摩院和中国互联网协会联合发起,致力于探讨最前沿的技术趋势与应用成果,搭建高质量的交流与分享平台,推动技术创新与产业应用链接,围绕“人工智能与新型计算”构建开放共享的开发者生态。

更多推荐