系列文章目录： 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 arcpy

preset = 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 = 1
    for 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 += 1


def 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_rule
    grouped_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_tifs


def 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_dir
    base = tif_names[0]
    out_coor_system = arcpy.Describe(base).spatialReference
    cell_width = arcpy.Describe(base).meanCellWidth
    band_count = arcpy.Describe(base).bandCount
    nums = len(groups)
    num = 1
    for 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 + 1


def 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 = 1
    for 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 + 1


def batch_clip_raster(rasters, out_dir, masks):
    nums = len(rasters) * len(masks)
    num = 1
    for 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 += 1


def 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 = 1
    for 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 + 1


def batch_setnull(rasters, out_dir, condition="VALUE>65528", prefix="sn_"):
    arcpy.CheckOutExtension("Spatial")
    nums = len(rasters)
    num = 1
    for 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 + 1


def 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)

    # step1
    s = 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()))

    # step2
    s = 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()))

    # step3
    s = 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()))

    # step4
    s = 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()))

    # step5
    s = 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)

    # step1
    s = 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()))

    # step2
    s = 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()))

    # step3
    s = 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()))

    # step4
    s = 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()))

    # step5
    s = 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()))

    # step6
    s = 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 controlling
  the 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 is
    called when the tool is opened."""
    for i in range(6, 17):
        self.params[i].category = "Advanced options"
    return

  def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""
    if self.params[0].altered:
        if self.params[0].value in ["MOD16_ET", "MOD16_PET"]:
            self.params[16].enabled = 1
        else:
            self.params[16].enabled = 0
    if self.params[0].value == "MOD13_NDVI":
        #self.params[5].value = "250 250"  # cell_size
        self.params[6].value = 0  # sds_index
        self.params[7].value = "NDVI"  # sds_name
        self.params[8].value = "16_BIT_SIGNED"
        self.params[9].value = 0.0001
    elif self.params[0].value == "MOD13_EVI":
        #self.params[5].value = "250 250"  # cell_size
        self.params[6].value = 1  # sds_index
        self.params[7].value = "EVI"  # sds_name
        self.params[8].value = "16_BIT_SIGNED"
        self.params[9].value = 0.0001
    elif self.params[0].value == "MOD16_ET":
        #self.params[5].value = "500 500"  # cell_size
        self.params[6].value = 0  # sds_index
        self.params[7].value = "ET"  # sds_name
        self.params[8].value = "16_BIT_UNSIGNED"
        self.params[9].value = 0.1
        self.params[16].value = "VALUE > 65528"
    elif self.params[0].value == "MOD16_PET":
        #self.params[5].value = "500 500"  # cell_size
        self.params[6].value = 2  # sds_index
        self.params[7].value = "PET"  # sds_name
        self.params[8].value = "16_BIT_UNSIGNED"
        self.params[9].value = 0.1
        self.params[16].value = "VALUE > 65528"
    return

  def updateMessages(self):
    """Modify the messages created by internal validation for each tool
    parameter.  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 输出文件预览

---