整体思路：HDF文件转tif文件（利用MRT） -- 用shp文件进行掩膜裁剪（两种方法） -- 设置空值 -- 8天的数据转为一个月 -- 单位转换

工具：MRT , GIS, Python

nasa对mod16数据的介绍：https://www.ntsg.umt.edu/project/modis/mod16.php

一：将HDF文件转为tif文件并进行拼接，得到一副tif图片

使用MRT进行

第一步：生成.prm文件

第二步：打开命令行，输下面的两个命令：

1

java -jar D:\XYX\MRT\bin\MRTBatch.jar -d D:\Data\mod09gq -p D:\Data\mod09gq\mod09gq.prm -o D:\Data\mod09gq

红的：hdf文件的位置

黄的：.prm文件的位置

蓝的：文件输出的位置

2

再输入MRTBatch.bat，然后等待处理

二：用GIS进行裁剪

1.导入流域shp文件

2.Toolboxes-spatial analyst tools-extraction-extract by mask-鼠标右键batch（批处理）

Input raster：输入tif文件   input raster or feature mask data：掩膜（流域shp） output raster：输出tif文件（需要自己明明为XXX.tif，该方法较麻烦）

三：用python进行裁剪

1.将源代码生成py文件

2.在catalog中添加mytoolbox（已完成）

3.右键mytoolbox-add-script-填写内容

py代码： 

#!/usr/bin/python
# -*- coding: UTF-8 -*-import os
import arcpy
from arcpy import envraster_path = arcpy.GetParameterAsText(0)
mask = arcpy.GetParameterAsText(1)
d_dir_name = arcpy.GetParameterAsText(2)def clip_batch(new_dir_name, old_prefix="Pr_", new_prefix="Clip_"):# type: () -> objectenv.workspace = raster_pathraster_List = arcpy.ListRasters("*", "tif")for raster in raster_List:if raster[:3] != old_prefix:out = "\\" + new_dir_name + "\\" + new_prefix + raster[:]elif raster[:3] == old_prefix:out = "\\" + new_dir_name + "\\" + new_prefix + raster[3:]arcpy.Clip_management(raster, "#", out, mask, "0", "ClippingGeometry")arcpy.AddMessage("Step2:Clip_"+raster[:]+"has done.")arcpy.SetProgressorPosition()arcpy.AddMessage("Step2:Completed")arcpy.ResetProgressor()
dir_name = d_dir_name
out_path = raster_path + "\\" + dir_name
ex = [2]
while True:if os.path.exists(out_path):dir_name = dir_name + "(" + str(ex[-1]) + ")"out_path = raster_path + "\\" + dir_nameex.append(ex[-1] + 1)continueelse:break
os.makedirs(out_path)
arcpy.AddMessage("Step1:Creating new folder named " + str(dir_name))
arcpy.AddMessage("Step1:Completed")
clip_batch(dir_name)

四：用python设置空置

方法同上

VALUE>10000 (value大于10000时置空)

#!/usr/bin/python
# -*- coding: UTF-8 -*-
import os
import arcpy
from arcpy import env
from arcpy.sa import *raster_path = arcpy.GetParameterAsText(0)
condition = arcpy.GetParameterAsText(1)
new_dir_name = arcpy.GetParameterAsText(2)def setN_batch(new_dir_name, new_prefix = "setN_"):env.workspace = raster_pathraster_list = arcpy.ListRasters()for inraster in raster_list:outraster = "\\" + new_dir_name + "\\" + new_prefix + inrasterarcpy.gp.SetNull_sa(inraster, inraster, outraster, condition)arcpy.AddMessage("Step2:"+ str(new_prefix) + inraster[:] + "has done.")arcpy.SetProgressorPosition()arcpy.AddMessage("Step2:Completed")arcpy.ResetProgressor()dir_name = new_dir_name
out_path = raster_path + "\\" + dir_name
os.makedirs(out_path)
arcpy.AddMessage("Step1:Creating new folder named " + str(dir_name))
arcpy.AddMessage("Step1:Completed")
if arcpy.CheckExtension("Spatial") == "Available":setN_batch(dir_name)
else:arcpy.AddMessage("Error!!! Spatial Analyst is unavailable")

五：用python将8天的数据转为月数据

import os
import time
import arcpydef show_files(path, out_files, suffix=".tif", out_type="path"):file_list = os.listdir(path)for file in file_list:cur_path = os.path.join(path, file)if os.path.isdir(cur_path):show_files(cur_path, out_files, out_type=out_type)else:if file.endswith(suffix):if out_type == "path":out_files.append(cur_path)elif out_type == "name":out_files.append(file)else:raise Exception("please select correct out_type value：path ；name")months_ping = {1: [1, 9, 17, 25],2: [25, 33, 41, 49, 57],3: [57, 65, 73, 81, 89],4: [89, 97, 105, 113],5: [121, 129, 137, 145],6: [145, 153, 161, 169, 177],7: [177, 185, 193, 201, 209],8: [209, 217, 225, 233, 241],9: [241, 249, 257, 265, 273],10: [273, 281, 289, 297],11: [305, 313, 321, 329],12: [329, 337, 345, 353, 361]}
months_run = {1: [1, 9, 17, 25],2: [25, 33, 41, 49, 57],3: [57, 65, 73, 81, 89],4: [89, 97, 105, 113, 121],5: [121, 129, 137, 145],6: [153, 161, 169, 177],7: [177, 185, 193, 201, 209],8: [209, 217, 225, 233, 241],9: [241, 249, 257, 265, 273],10: [273, 281, 289, 297, 305],11: [305, 313, 321, 329],12: [329, 337, 345, 353, 361]}in_path = arcpy.GetParameterAsText(0)  # r"H:\taobao\suisdio\loess"
out_path = arcpy.GetParameterAsText(1)
pixel_type = arcpy.GetParameterAsText(2)
colormap_mode = arcpy.GetParameterAsText(3)arcpy.env.workspace = in_path
all_tifs = []
show_files(in_path, all_tifs, out_type="name")
name_map = {}
avail_years = []
for fname in all_tifs:y = int(fname.split(".")[1][1:5])d = int(fname.split(".")[1][-3:])keyname = "{0}.{1}".format(y, d)name_map[keyname] = fnameif y not in avail_years:avail_years.append(y)base = all_tifs[0]
out_coor_system = arcpy.Describe(base).spatialReference
cell_width = arcpy.Describe(base).meanCellWidth
band_count = arcpy.Describe(base).bandCountgroups = {}
for year in avail_years:months_temp = months_pingif year % 4 == 0:months_temp = months_runfor month in months_temp.keys():new_name = "A{0}M{1}.tif".format(year, month)target_days = months_temp[month]groups[new_name] = []for i in target_days:fkey = "{0}.{1}".format(year, i)if fkey in name_map:groups[new_name].append(name_map[fkey])nums = len(groups)
num = 1
for i in groups:s = time.time()month = int(i.split(".")[0][6:])  # A2000M2.tifif len(groups[i]) == 0:arcpy.AddMessage("{0}/{1} | {2} is None".format(num, nums, i))num = num + 1continuegroups[i] = ';'.join(groups[i])if not os.path.exists(os.path.join(out_path, i)):if month != 12:arcpy.MosaicToNewRaster_management(groups[i], out_path, i, out_coor_system, pixel_type, cell_width, band_count,"MEAN", colormap_mode)else:arcpy.MosaicToNewRaster_management(groups[i], out_path, i, out_coor_system, pixel_type, cell_width, band_count,"SUM", colormap_mode)e = time.time()arcpy.AddMessage("{0}/{1} | {2} Completed, time used {3}s".format(num, nums, i, e - s))else:e = time.time()arcpy.AddMessage("{0}/{1} | {2} existed, , time used {3}s".format(num, nums, i, e - s))num = num + 1

六：用python统一乘0.1 ，转换单位为mm/day

#!/usr/bin/python
# -*- coding: UTF-8 -*-import arcpy
from arcpy import env
from arcpy.sa import *
import osraster_path = arcpy.GetParameterAsText(0)
new_dir_name = arcpy.GetParameterAsText(1)def multiply_batch(new_dir_name, new_prefix="united_"):env.workspace = raster_pathraster_list = arcpy.ListRasters()for inraster in raster_list:year = int(inraster.split(".")[0][1:5]) # 年份在文件名中的位置month = int(inraster.split(".")[0][6:]) # 日序在文件名中的位置if month == 12:if year % 4 == 0:scale_factor = 0.1 / 38else:scale_factor = 0.1 / 37else:scale_factor = 0.1 / 8outraster = "\\" + new_dir_name + "\\" + new_prefix + inrasterarcpy.gp.Times_sa(inraster, scale_factor, outraster)arcpy.AddMessage("Step2:" + str(new_prefix) + inraster[:] + "  has done.")arcpy.SetProgressorPosition()arcpy.AddMessage("Step2:Completed")arcpy.ResetProgressor()dir_name = new_dir_name
out_path = raster_path + "\\" + dir_name
os.makedirs(out_path)
arcpy.AddMessage("Step1:Creating new folder named " + str(dir_name))
arcpy.AddMessage("Step1:Completed")
if arcpy.CheckExtension("Spatial") == "Available":multiply_batch(dir_name)
else:arcpy.AddMessage("Error!!! Spatial Analyst is unavailable")

参考链接：

https://blog.csdn.net/qq_37948866/article/details/111184481

https://blog.csdn.net/qq_37948866/article/details/99889247

https://blog.csdn.net/qq_37948866/article/details/100098259

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！