文章目录

一、功能说明
二、代码
三、工具设置
- 3.1 参数(parameters)设置
- 3.2 验证(validation)设置
四、工具界面
五、使用教程
- 5.1 准备
- 5.2 操作步骤
- 5.3 运行截图
- 5.4 输出文件预览

系列文章目录： MODIS数据处理

一、功能说明

批量加工MODIS数据，目前支持的产品包括：

MOD13：NDVI、EVI
MOD16：ET、PET

输入：
下载的hdf文件

输出：
研究区的栅格影像

二、代码

#!/usr/bin/python
# -*- coding: UTF-8 -*-
"""
@File    : easyModis.py
@Author  : salierib
@Time    : 2021/11/5 11:44
@Version: 1.0
"""import os
import time
import arcpypreset = arcpy.GetParameterAsText(0)
workspace = arcpy.GetParameterAsText(1)
hdfs = arcpy.GetParameterAsText(2)
masks = arcpy.GetParameterAsText(3)
out_coor_system = arcpy.GetParameterAsText(4)
cell_size = arcpy.GetParameterAsText(5)
sds_index = arcpy.GetParameterAsText(6)
sds_name = arcpy.GetParameterAsText(7)
pixel_type = arcpy.GetParameterAsText(8)
scale_factor = arcpy.GetParameterAsText(9)
mosaic_method = arcpy.GetParameterAsText(10)
colormap_mode = arcpy.GetParameterAsText(11)
resampling_type = arcpy.GetParameterAsText(12)
pr_prefix = arcpy.GetParameterAsText(13)
scale_prefix = arcpy.GetParameterAsText(14)
sn_prefix = arcpy.GetParameterAsText(15)
condition = arcpy.GetParameterAsText(16)hdfs = hdfs.split(";")
masks = masks.split(";")def batch_extract_sds(hdfs, out_dir, sds_index=0, suffix="NDVI"):nums = len(hdfs)num = 1for hdf in hdfs:s = time.time()base_name = os.path.splitext(os.path.basename(hdf))[0]out_tif = os.path.join(out_dir, base_name + "." + "{0}.tif".format(suffix))if not os.path.exists(out_tif):try:arcpy.ExtractSubDataset_management(hdf, out_tif, sds_index)e = time.time()arcpy.AddMessage("%d/%d | %s completed, time used %.2fs" % (num, nums, out_tif, e - s))except Exception as err:arcpy.AddMessage("%d/%d | %s errored, %s" % (num, nums, out_tif, err))else:arcpy.AddMessage("%d/%d | %s already exists" % (num, nums, out_tif))num += 1def normal_mosaic_rule(fname):return '.'.join(fname.split('.')[:2]) + '.' + '.'.join(fname.split('.')[-2:])def group_tifs(tif_names, group_func="mosaic"):if group_func == "mosaic":group_func = normal_mosaic_rulegrouped_tifs = {}for tif_name in tif_names:k = group_func(tif_name)if k in grouped_tifs:grouped_tifs[k].append(tif_name)else:grouped_tifs[k] = []grouped_tifs[k].append(tif_name)return grouped_tifsdef batch_mosaic(in_dir, out_dir, groups=None, pixel_type="16_BIT_SIGNED", mosaic_method="MAXIMUM",colormap_mode="FIRST"):tif_names = [n for n in os.listdir(in_dir) if n.endswith(".tif")]if groups is None:groups = group_tifs(tif_names, group_func="mosaic")arcpy.env.workspace = in_dirbase = tif_names[0]out_coor_system = arcpy.Describe(base).spatialReferencecell_width = arcpy.Describe(base).meanCellWidthband_count = arcpy.Describe(base).bandCountnums = len(groups)num = 1for i in groups:s = time.time()groups[i] = ';'.join(groups[i])if not os.path.exists(os.path.join(out_dir, i)):try:arcpy.MosaicToNewRaster_management(groups[i], out_dir, i, out_coor_system, pixel_type, cell_width,band_count, mosaic_method, colormap_mode)e = time.time()arcpy.AddMessage("%d/%d | %s completed, time used %.2fs" % (num, nums, i, e - s))except Exception as err:arcpy.AddMessage("%d/%d | %s errored, %s" % (num, nums, i, err))else:arcpy.AddMessage("%d/%d | %s already exists" % (num, nums, i))num = num + 1def batch_project_raster(rasters, out_dir, prefix="pr_", out_coor_system="WGS_1984.prj",resampling_type="NEAREST", cell_size="#"):if arcpy.CheckExtension("Spatial") != "Available":arcpy.AddMessage("Error!!! Spatial Analyst is unavailable")nums = len(rasters)num = 1for raster in rasters:s = time.time()raster_name = os.path.split(raster)[1]out_raster = os.path.join(out_dir, prefix + raster_name)if not os.path.exists(out_raster):try:arcpy.ProjectRaster_management(raster, out_raster, out_coor_system, resampling_type, cell_size, "#","#", "#")e = time.time()arcpy.AddMessage("%d/%d | %s completed, time used %.2fs" % (num, nums, out_raster, e - s))except Exception as err:arcpy.AddMessage("%d/%d | %s errored, %s" % (num, nums, out_raster, err))else:arcpy.AddMessage("%d/%d | %s already exists" % (num, nums, raster))num = num + 1def batch_clip_raster(rasters, out_dir, masks):nums = len(rasters) * len(masks)num = 1for i, mask in enumerate(masks):mask_name = os.path.splitext(os.path.basename(mask))[0]for raster in rasters:s = time.time()old_raster_name = os.path.splitext(os.path.basename(raster))[0]new_raster_name = "{0}_{1}.tif".format(mask_name, old_raster_name.split("_")[-1])out_raster = os.path.join(out_dir, new_raster_name)if not os.path.exists(out_raster):try:arcpy.Clip_management(raster, "#", out_raster, mask, "#", "ClippingGeometry")e = time.time()arcpy.AddMessage("%d/%d | %s completed, time used %.2fs" % (num, nums, out_raster, e - s))except Exception as err:arcpy.AddMessage("%d/%d | %s errored, %s" % (num, nums, out_raster, err))else:arcpy.AddMessage("%d/%d | %s already exists" % (num, nums, out_raster))num += 1def batch_multiply(rasters, out_dir, scale_factor=0.0001, prefix="scaled_"):scale_factor = str(scale_factor)arcpy.CheckOutExtension("Spatial")if arcpy.CheckExtension("Spatial") != "Available":arcpy.AddMessage("Error!!! Spatial Analyst is unavailable")nums = len(rasters)num = 1for raster in rasters:s = time.time()raster_name = os.path.split(raster)[1]out_raster = os.path.join(out_dir, prefix + raster_name)if not os.path.exists(out_raster):try:arcpy.gp.Times_sa(raster, scale_factor, out_raster)e = time.time()arcpy.AddMessage("%d/%d | %s completed, time used %.2fs" % (num, nums, out_raster, e - s))except Exception as err:arcpy.AddMessage("%d/%d | %s errored, %s" % (num, nums, out_raster, err))else:arcpy.AddMessage("%d/%d | %s already exists" % (num, nums, out_raster))num = num + 1def batch_setnull(rasters, out_dir, condition="VALUE>65528", prefix="sn_"):arcpy.CheckOutExtension("Spatial")nums = len(rasters)num = 1for raster in rasters:s = time.time()raster_name = os.path.split(raster)[1]out_raster = os.path.join(out_dir, prefix + raster_name)if not os.path.exists(out_raster):try:arcpy.gp.SetNull_sa(raster, raster, out_raster, condition)e = time.time()arcpy.AddMessage("%d/%d | %s completed, time used %.2fs" % (num, nums, out_raster, e - s))except Exception as err:arcpy.AddMessage("%d/%d | %s errored, %s" % (num, nums, out_raster, err))else:arcpy.AddMessage("%d/%d | %s already exists" % (num, nums, out_raster))num = num + 1def find_tifs(in_dir):return [os.path.join(in_dir, fname) for fname in os.listdir(in_dir) if fname.endswith(".tif")]def localtime():return time.strftime("%Y-%m-%d %H:%M:%S", time.localtime())def mod13preprocess(workspace, hdfs, masks, out_coor_system, cell_size="#",sds_index=0, sds_name="NDVI",pixel_type="16_BIT_SIGNED", mosaic_method="LAST", colormap_mode="FIRST",pr_prefix="pr_", resampling_type="NEAREST",scale_prefix="", scale_factor=0.0001):if not os.path.exists(workspace):os.mkdir(workspace)dir_names = ["1_extract", "2_mosaic", "3_reproject", "4_clip", "5_scale"]dirs = [os.path.join(workspace, name) for name in dir_names]for dir in dirs:if not os.path.exists(dir):os.mkdir(dir)# step1s = time.time()arcpy.AddMessage("Starting step 1/5: extract subdataset into {0}... {1}".format(dirs[0], localtime()))batch_extract_sds(hdfs, dirs[0], sds_index=sds_index, suffix=sds_name)e = time.time()arcpy.AddMessage("Time for step1 = {0} seconds. {1}\n".format(e - s, localtime()))# step2s = time.time()arcpy.AddMessage("Starting step 2/5: mosaic raster into {0}... {1}".format(dirs[1], localtime()))batch_mosaic(dirs[0], dirs[1], pixel_type=pixel_type, mosaic_method=mosaic_method, colormap_mode=colormap_mode)e = time.time()arcpy.AddMessage("Time for step2 = {0} seconds. {1}\n".format(e - s, localtime()))# step3s = time.time()arcpy.AddMessage("Starting step 3/5: reproject raster into {0}... {1}".format(dirs[2], localtime()))rasters = find_tifs(dirs[1])batch_project_raster(rasters, dirs[2], prefix=pr_prefix, out_coor_system=out_coor_system, resampling_type=resampling_type, cell_size=cell_size)e = time.time()arcpy.AddMessage("Time for step3 = {0} seconds. {1}\n".format(e - s, localtime()))# step4s = time.time()arcpy.AddMessage("Starting step 4/5: clip raster into {0}... {1}".format(dirs[3], localtime()))rasters = find_tifs(dirs[2])batch_clip_raster(rasters, dirs[3], masks=masks)e = time.time()arcpy.AddMessage("Time for step4 = {0} seconds. {1}\n".format(e - s, localtime()))# step5s = time.time()arcpy.AddMessage("Starting step 5/5:raster times scale factor into {0}... {1}".format(dirs[4], localtime()))tifs = find_tifs(dirs[3])batch_multiply(tifs, out_dir=dirs[4], prefix=scale_prefix, scale_factor=scale_factor)e = time.time()arcpy.AddMessage("Time for step5 = {0} seconds. {1}\n".format(e - s, localtime()))def mod16preprocess(workspace, hdfs, masks, out_coor_system, cell_size="#",sds_index=0, sds_name="ET",pixel_type="16_BIT_UNSIGNED", mosaic_method="LAST", colormap_mode="FIRST",pr_prefix="pr_", resampling_type="NEAREST",sn_prefix="sn_", condition="VALUE > 65528",scale_prefix="", scale_factor=0.1):if not os.path.exists(workspace):os.mkdir(workspace)dir_names = ["1_extract", "2_mosaic", "3_reproject", "4_clip", "5_setn", "6_scale"]dirs = [os.path.join(workspace, name) for name in dir_names]for dir in dirs:if not os.path.exists(dir):os.mkdir(dir)# step1s = time.time()arcpy.AddMessage("Starting step 1/6: extract subdataset into {0}... {1}".format(dirs[0], localtime()))batch_extract_sds(hdfs, dirs[0], sds_index=sds_index, suffix=sds_name)e = time.time()arcpy.AddMessage("Time for step1 = {0} seconds. {1}\n".format(e - s, localtime()))# step2s = time.time()arcpy.AddMessage("Starting step 2/6: mosaic raster into {0}... {1}".format(dirs[1], localtime()))batch_mosaic(dirs[0], dirs[1], pixel_type=pixel_type, mosaic_method=mosaic_method, colormap_mode=colormap_mode)e = time.time()arcpy.AddMessage("Time for step2 = {0} seconds. {1}\n".format(e - s, localtime()))# step3s = time.time()arcpy.AddMessage("Starting step 3/6: reproject raster into {0}... {1}".format(dirs[2], localtime()))rasters = find_tifs(dirs[1])batch_project_raster(rasters, dirs[2], prefix=pr_prefix, out_coor_system=out_coor_system, resampling_type=resampling_type, cell_size=cell_size)e = time.time()arcpy.AddMessage("Time for step3 = {0} seconds. {1}\n".format(e - s, localtime()))# step4s = time.time()arcpy.AddMessage("Starting step 4/6: clip raster into {0}... {1}".format(dirs[3], localtime()))rasters = find_tifs(dirs[2])batch_clip_raster(rasters, dirs[3], masks=masks)e = time.time()arcpy.AddMessage("Time for step4 = {0} seconds. {1}\n".format(e - s, localtime()))# step5s = time.time()arcpy.AddMessage("Starting step 5/6: exclude invalid value, into {0}... {1}".format(dirs[4], localtime()))rasters = find_tifs(dirs[3])batch_setnull(rasters, dirs[4], condition=condition, prefix=sn_prefix)e = time.time()arcpy.AddMessage("Time for step5 = {0} seconds. {1}\n".format(e - s, localtime()))# step6s = time.time()arcpy.AddMessage("Starting step 6/6: raster times scale factor into {0}... {1}".format(dirs[5], localtime()))tifs = find_tifs(dirs[4])batch_multiply(tifs, out_dir=dirs[5], prefix=scale_prefix, scale_factor=scale_factor)e = time.time()arcpy.AddMessage("Time for step5 = {0} seconds. {1}\n".format(e - s, localtime()))if preset in ["MOD13_NDVI", "MOD13_EVI"]:mod13preprocess(workspace=workspace,hdfs=hdfs,masks=masks,out_coor_system=out_coor_system,cell_size=cell_size,sds_index=sds_index,sds_name=sds_name,pixel_type=pixel_type,mosaic_method=mosaic_method,colormap_mode=colormap_mode,resampling_type=resampling_type,scale_factor=scale_factor,scale_prefix=scale_prefix,pr_prefix=pr_prefix)
elif preset in ["MOD16_ET", "MOD16_PET"]:mod16preprocess(workspace=workspace,hdfs=hdfs,masks=masks,out_coor_system=out_coor_system,cell_size=cell_size,sds_index=sds_index,sds_name=sds_name,pixel_type=pixel_type,mosaic_method=mosaic_method,colormap_mode=colormap_mode,resampling_type=resampling_type,scale_factor=scale_factor,scale_prefix=scale_prefix,pr_prefix=pr_prefix,sn_prefix=sn_prefix,condition=condition)

三、工具设置

关于如何创建并使用脚本，可以参考这本书的第13章
[1]赞德伯根. 面向ArcGIS的Python脚本编程[M]. 人民邮电出版社, 2014.

3.1 参数(parameters)设置

3.2 验证(validation)设置

import arcpy
class ToolValidator(object):"""Class for validating a tool's parameter values and controllingthe behavior of the tool's dialog."""def __init__(self):"""Setup arcpy and the list of tool parameters."""self.params = arcpy.GetParameterInfo()def initializeParameters(self):"""Refine the properties of a tool's parameters.  This method iscalled when the tool is opened."""for i in range(6, 17):self.params[i].category = "Advanced options"returndef updateParameters(self):"""Modify the values and properties of parameters before internalvalidation is performed.  This method is called whenever a parameterhas been changed."""if self.params[0].altered:if self.params[0].value in ["MOD16_ET", "MOD16_PET"]:self.params[16].enabled = 1else:self.params[16].enabled = 0if self.params[0].value == "MOD13_NDVI":#self.params[5].value = "250 250"  # cell_sizeself.params[6].value = 0  # sds_indexself.params[7].value = "NDVI"  # sds_nameself.params[8].value = "16_BIT_SIGNED"self.params[9].value = 0.0001elif self.params[0].value == "MOD13_EVI":#self.params[5].value = "250 250"  # cell_sizeself.params[6].value = 1  # sds_indexself.params[7].value = "EVI"  # sds_nameself.params[8].value = "16_BIT_SIGNED"self.params[9].value = 0.0001elif self.params[0].value == "MOD16_ET":#self.params[5].value = "500 500"  # cell_sizeself.params[6].value = 0  # sds_indexself.params[7].value = "ET"  # sds_nameself.params[8].value = "16_BIT_UNSIGNED"self.params[9].value = 0.1self.params[16].value = "VALUE > 65528"elif self.params[0].value == "MOD16_PET":#self.params[5].value = "500 500"  # cell_sizeself.params[6].value = 2  # sds_indexself.params[7].value = "PET"  # sds_nameself.params[8].value = "16_BIT_UNSIGNED"self.params[9].value = 0.1self.params[16].value = "VALUE > 65528"returndef updateMessages(self):"""Modify the messages created by internal validation for each toolparameter.  This method is called after internal validation."""return

四、工具界面

五、使用教程

工作目录建议使用英文，ArcGIS建议使用英文版，否则可能报错ascii error。

5.1 准备

以MYD16A3加工新疆地区的ET500m栅格为例，需要准备的输入文件包括：

MYD16A3数据(.hdf)
新疆地区边界文件(.shp)

5.2 操作步骤

具体的处理步骤如下：

1、将预设参数设置为MOD16_ET
2、新建一个文件夹作为工作目录，建议使用固态。S:\output
3、选择要处理的hdf文件
4、选择裁剪边界文件
5、设置目标坐标系，强烈建议将目标坐标系与边界文件所在的坐标系保持一致
6、设置栅格分辨率为500m
7、如果是新手，不建议修改Advanced options中的参数，其参数会根据所选的预设自动修改

5.3 运行截图

运行消息主要显示分步骤的用时等信息

5.4 输出文件预览

【MODIS数据处理#14】拼接、投影、裁剪一键完成，比MRT更方便的ArcGIS脚本工具（含代码）相关推荐

windows下9款一键快速搭建PHP本地运行环境的好工具(含php7.0环境)
推荐几款一键快速搭建PHP本地运行环境的好工具(含php7.0及apache,nigix,mysql) 首推phpstudy2016和wampServer3.0.6 理由支持php7.0 目前 ...
【MODIS数据处理#13】使用Arcpy一键加工长时间序列MODIS数据
文章目录一.前言 1.1 回顾 1.2 在Pycharm使用Arcpy的方法二.功能简介三.代码 3.1 文件搜索file_picker.py 3.2 批处理工具batch_tools.py 2 ...
【MODIS数据处理#15】分享一个自制的MODIS数据处理工具箱
文章目录一.下载地址二.工具箱内容三.配置教程四.使用教程后记整理了本人自制的MODIS数据批处理脚本工具,以ArcGIS共享工具箱(.tbx)的方式免费分享给大家.所有工具都有详细的说明 ...
IDL实现MODIS Grid（正弦投影）产品的重投影及拼接处理
IDL实现MODIS Grid(正弦投影)产品的重投影及拼接处理前言 map_proj_image函数使用关键单个文件的重投影示例多个文件的重投影+拼接后记前言关于MODIS正弦投影产品的 ...
[MODIS数据处理#9]例四：基于MCD12Q2数据集初步分析中国植被物候空间分布特征
文章目录一.MCD12Q2介绍二.问题描述三.预处理 1.提取子数据集 2.波段提取 3.拼接 4.栅格投影 5.裁剪四.获取生长季长度五.获取生长季始.末所在日期系列文章目录: MODI ...
[MODIS数据处理#8]批量将ET栅格的时间分辨率从8-day转换为monthly的一种思路
文章目录一.问题描述二.数据预处理 1.下载数据 2.提取子数据集.镶嵌.重投影 3.裁剪 4.排除特殊值区域三.按月镶嵌至新栅格 1.批量按月镶嵌脚本 1.1.功能介绍 1.2.脚本代码 1. ...
[MODIS数据处理#5]例二：将ET栅格的时间分辨率从8-day转换为monthly
一.下载目前EarthData提供的MOD16产品只有MOD16A2.MOD16A3两种 MOD16A2: 8-day,500m MOD16A3: yearly, 500m GMAO美国蒙大拿大学森 ...
【MODIS数据处理#12】例七：基于MOD09Q1数据集合成NDVI
文章目录一.MOD09Q1数据介绍二.处理步骤 2.1 MRT工具预处理 2.2 去除无效值 2.3 计算NDVI 系列文章目录: MODIS数据处理一.MOD09Q1数据介绍官网介绍链接:h ...
ArcGIS：如何进行栅格数据的拼接和裁剪、坡度坡向的提取、地形透视图的建立、等高线的提取、剖面图的创建？
目录 01 说明 02 实验目的及要求 03 实验设备及软件平台 04 实验内容与步骤 4.1 DEM 数据拼接和裁剪 4.2 地形属性的提取 4.3 透视图的建立(均在ArcScence中操作) 4 ...

【MODIS数据处理#14】拼接、投影、裁剪一键完成，比MRT更方便的ArcGIS脚本工具（含代码）