北师大GLASS LAI数据的处理流程及使用北师大GLASS LAI数据进行月度数据合成（最大值）

由于驱动数据的要求，之前下载了北师大的GLASS LAI数据 1KM数据进行处理，具体是经过了

1.下载，地址http://glass-product.bnu.edu.cn/index.html

2.数据的格式转换及拼接，原始数据并非TIF格式，而为hdf 可以使用 arcpy.ExtractSubDataset_management工具进行批量转化，之后再进行拼接，这样就可以得到全球的TIF格式LAI数据

3数据的裁剪，根据自己研究区域 对数据进行裁剪（把这一步放在前面可以减少后续处理所需的处理量）

import arcpy
import os
Coordinate_System = "PROJCS['WGS_1984_Web_Mercator_Auxiliary_Sphere',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator_Auxiliary_Sphere'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Standard_Parallel_1',0.0],PARAMETER['Auxiliary_Sphere_Type',0.0],UNIT['Meter',1.0]]"#这里使用的是Web84墨卡托投影坐标
inputdir=r"——输入文件夹——"
output_path = r"——输出文件夹——"
mask = "——掩膜文件——"
for time in range(2000,2019):
    seconddir = inputdir + os.sep + str(time)
    thirddirs = os.listdir(seconddir)
    for thirddir in thirddirs:
        wholedir = seconddir + os.sep + thirddir
        print wholedir
        arcpy.env.workspace = wholedir
        rasterList = arcpy.ListRasters("*", "TIF")
        print rasterList
        for raster in rasterList:
            outname = output_path + os.sep + str(time)+raster[0:3] + "Web84" + ".tif"
            if os.path.exists(outname):
                print (outname + "have exist")
            else:
                arcpy.ProjectRaster_management(raster, outname, Coordinate_System, "NEAREST", "#", "#", "#", "#")
                print (outname + "have been done")
            out = output_path + os.sep + str(time)+raster[0:3] + ".tif"
            arcpy.Clip_management(outname, "#", out, mask, "-999", "ClippingGeometry","NO_MAINTAIN_EXTENT")
            print out + " clip has done!"

4数据的坐标系处理（可选）：如果有需要。

5这里前置经过拼接、裁剪和坐标系处理，并且每张TIF的文件名已经处理成了年份+日的格式，如2000年的第一天为“2000001.tif”并且2000-2017年所有的数据都放在输入文件夹下，那么下面的脚本即可快速进行MVC最大值合成。

import arcpy
#输出坐标系信息
Coordinate_System = "PROJCS['WGS_1984_Web_Mercator_Auxiliary_Sphere',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator_Auxiliary_Sphere'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Standard_Parallel_1',0.0],PARAMETER['Auxiliary_Sphere_Type',0.0],UNIT['Meter',1.0]]"

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 = r"H:\GLASS\LAI1km_yangtze" #所有原始TIF文件输入路径
out_path = r"H:\GLASS\LAI1km"  #输出路径
arcpy.env.workspace = in_path #将arcpy工作空间转为输入路径

for year in range(2001, 2018):
    #若日期为润年
    if year % 4 == 0:
        print (str(year) + "是润年")
        for month in range(1, 13):
            tif = []
            for day in months_run[month]:
                day_num = str(day).zfill(3)  # 数字补0，补成想要的位数这里为3位
                tif.append(str(year) + str(day_num) + ".tif")
            print ("对"+str(year)+"年"+str(month)+"月执行最大值合成")
            arcpy.MosaicToNewRaster_management(tif, out_path, str(year)+'_'+str(month).zfill(2)+".tif", Coordinate_System, "8_BIT_UNSIGNED", "#", 1,
                                               "MAXIMUM", "FIRST")
    else:
        print (str(year) + "不是润年")
        for month in range(1, 13):
            tif = []
            print month
            for day in months_ping[month]:
                day_num = str(day).zfill(3)  # 数字补0，补成想要的位数这里为3位
                tif.append(str(year) + str(day_num) + ".tif")
            print ("对"+str(year)+"年"+str(month)+"月执行最大值合成")
            arcpy.MosaicToNewRaster_management(tif, out_path, str(year)+'_'+str(month).zfill(2)+".tif", Coordinate_System, "8_BIT_UNSIGNED", "#", 1,
                                               "MAXIMUM", "FIRST")

结果部分如下：

参考文章

[MODIS数据处理#8]批量将ET栅格的时间分辨率从8-day转换为monthly的一种思路