【ArcGIS自定义脚本工具】MOD13Q1数据集处理：批量将16-day分辨率的NDVI栅格转换为monthly分辨率

系列文章目录： ArcGIS自定义脚本编程

---

一、功能介绍

利用Arcmap中的镶嵌至新栅格工具，镶嵌运算符选择MAXIMUM，将由MOD13Q1数据集预处理得到的16-day、250m分辨率的NDVI栅格批量转为monthly逐月分辨率。

进行镶嵌至新栅格时，平闰年1 ~ 12月份NDVI的输入栅格可参考下图进行选择：

---

二、脚本代码

#!/usr/bin/python
# -*- coding: UTF-8 -*-

import os
import time
import arcpy


def 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, 17],
               2: [17, 33, 49],
               3: [49, 65, 81],
               4: [81, 97, 113],
               5: [113, 129, 145],
               6: [145, 161, 177],
               7: [177, 193, 209],
               8: [209, 225, 241],
               9: [241, 257, 273],
               10: [273, 289],
               11: [305, 321],
               12: [321, 337, 353]}
months_run = {1: [1, 17],
              2: [17, 33, 49],
              3: [49, 65, 81],
              4: [81, 97, 113],
              5: [113, 129, 145],
              6: [145, 161, 177],
              7: [177, 193, 209],
              8: [209, 225, 241],
              9: [241, 257, 273],
              10: [273, 289, 305],
              11: [305, 321],
              12: [321, 337, 353]}

in_path = arcpy.GetParameterAsText(0)  # r"H:\taobao\suisdio\loess"
out_path = arcpy.GetParameterAsText(1)
pixel_type = arcpy.GetParameterAsText(2)
mosaic_method = arcpy.GetParameterAsText(3)
colormap_mode = arcpy.GetParameterAsText(4)

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] = fname
    if 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).bandCount

groups = {}
for year in avail_years:
    months_temp = months_ping
    if year % 4 == 0:
        months_temp = months_run
    for 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()
    if len(groups[i]) == 0:
        arcpy.AddMessage("{0}/{1} | {2} is None".format(num, nums, i))
        num = num + 1
        continue
    groups[i] = ';'.join(groups[i])
    if not os.path.exists(os.path.join(out_path, i)):
        arcpy.MosaicToNewRaster_management(groups[i], out_path, i, out_coor_system, pixel_type, cell_width, band_count,
                                           mosaic_method, 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

---

三、工具参数

---

四、工具界面

---