Python实现NDVI的MVC算法

因为所下数据为NDVI旬数据，合并为月数据时采用国际通用最大值合成法（MVC）将旬数据合并成月数据，因数据较多，遂考虑用代码解决。主要分为遍历文件，读取制定名称的栅格文件，以及ArcGIS工具箱中的Raster Calculator的con函数调用等几大部分。因为用的for循环结构，所以原有文件的命名要有一定的规律可循。如果源文件的名称较杂，可以用bat结合excel先进行批量重命名（具体方法网上有，本人亦可能追加补充）。

话不多说，上代码：

import os
import sys
import arcpy
import numpy
from arcpy import Raster
os.chdir("D:/###/ndvi20171128/raster")  
os.getcwd()

data1_list = os.listdir("./ndvi0")
data2_list = os.listdir("./ndvi5")

for year in range(1982,2013):
    import numpy
    # Caution: the "number" in below sentence will be start from 1 to 12 which less than 13
    for number in range(1,13):
        # Assign data from filename list
        data1 = [data for data in data1_list
                 if data[0:] == (str(year) + "_" + str(number) + "_ndvi" + ".tif")]
        # Verifying if the data in a sure date is empty 
        if len(data1) == 0:
            continue
        else:
            data1 = data1[0]

        data2 = [data for data in data2_list
                 if data[0:] == (str(year) + "_" + str(number) + ".5_ndvi" + ".tif")]
        
        if len(data2) == 0:
            continue
        else:
            data2 = data2[0]
        
        # Local variables:
        raster1 = Raster("ndvi0/" + data1)
        raster2 = Raster("ndvi5/" + data2)
        rastercalc6 = os.path.join("D:/###/ndvi20171128/raster/mndvi","mndvi" + str(year) + str(number) + ".img")
        
        # Process: Raster Calculator
        arcpy.gp.RasterCalculator_sa("Con((raster1 >=  - 1) &(raster2 >=  - 1) ,Con(raster1 > raster2,raster1,raster2))", rastercalc6)

注：

1. 本文数据组织格式为：

1982_1_ndvi.tif; 1982_1.5_ndvi.tif; 1982_2_ndvi.tif; ……

2. 最好在ArcGIS或Catalog中的python窗口中运行，以免出错。

3. 因为调用了Raster Calculator工具的原因，在ArcGIS中运行会将结果直接加载到数据窗口中，导致大数据处理的后续运行会变得异常缓慢，建议在Catalog中运行相应代码。