批量下载并处理GLASS(modis/AVHRR)产品数据(二)
·
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文件夹:


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

所有评论(0)