AnsweredAssumed Answered

Python routine to calculate MODIS NDVI anomaly

Question asked by vianna on Jun 14, 2018
Latest reply on Jun 21, 2018 by vianna

Xander Bakker

Python script to calculate NDVI anomalies

https://geonet.esri.com/people/xander_bakker

 

I´d like to automate all my workflow, from downloaded original raster datasets to anomaly calculator:

 

3 years data and folder structure can be downloaded here.

 

Question 1: What is the best option: to run separated scripts or create a unique code?

Question 2: How to translate the model builder to python script (Python script 2)?

Question 3: If the best option is to run separated, how to make it in a logical sequence (script 1 -> script 2 -> script 3)

 

Summary workflow:

  1. Select and download NDVI filtered raster datasets;
    1. Store downloaded original raster datasets at X:\Modis250\Originais\
    2. Organizing files and folders: (Python script 1)
    3. Rename original raster datasets (MODIS.ndvi.YYYYjjj.yL6000.BOKU.tif) to filename standard “ndvi_YYYY_dec_dd.tif”, where “YYYY” is the year, “jjj” is the Julian day and “dd” is the decade (ten days period).
    4. Transfer renamed files from X:\Modis250\Originais to X:\Modis250\NDVI10dias
  2. Pre-processing raster datasets: (Python script 2)
    1. Project renamed raster datasets ("ndvi_YYYY_dec_dd.tif") from Lat/Lon WGS 1984 to Lat/Lon SIRGAS2000 using "SAD_1969_To_WGS_1984_14 + SAD_1969_To_SIRGAS_2000_1" geographic transformation.
    2. Overwrite spurious data values: NDVI = Con("%NDVI%"<1,0,Con("%NDVI%">9998,9999,"%NDVI%"))
    3. Standardize raster datasets: NDVI = (NDVI – Mean(NDVI))/Std(NDVI)
    4. Scale raster datasets to 0 – 100: NDVI = ((NDVI-Min(NDVI))/(Max(NDVI)-Min(NDVI)) * 100
    5. Mask standardized raster datasets by the area of interest (X:\Modis250\SRTM\mde_sc_90.tif) and output to X:\Modis250\NDVI_10dias_norm
  3. Processing data (Python script 3)
    1. Calculate NDVI mean for each decadal: mean_dec_dd.tif = MEAN (ndvi_YYY1_dec_dd.tif, ndvi_YYY2_dec_dd.tif, ndvi_YYY3_dec_dd.tif,…)
    2. Calculate NDVI anomaly for each year by decade: anom_yyyy_dec_dd.tif = ndvi_YYYY_dec_dd.tif - mean_dec_dd.tif

 

Python scripts 1 and 3 are working. The scrip 2 (Pre-processing) is a model in a toolbox.

 

Python script 1: written by myself, maybe need to be improved!

#Rename original files
import arcpy, datetime, glob, os
arcpy.env.overwriteOutput = True
arcpy.env.workspace=r'O:\Modis250\Originais'
infolder='O:\Modis250\Originals'
outfolder='O:\Modis250\NDVI_10dias'
list1=arcpy.ListRasters("*.tif")
list1.sort()
   
for raster in list1:
     source_path = os.path.join(infolder, raster)
     oldFilename=raster
     dd = int(oldFilename[15:18])
     if dd < 11: dd = 1
     elif dd < 21: dd = 2
     elif dd < 31: dd = 3
     elif dd < 41: dd = 4
     elif dd < 51: dd = 5
     elif dd < 61: dd = 6
     elif dd < 71: dd = 7
     elif dd < 81: dd = 8
     elif dd < 91: dd = 9
     elif dd < 101: dd = 10
     elif dd < 111: dd = 11
     elif dd < 121: dd = 12
     elif dd < 131: dd = 13
     elif dd < 141: dd = 14
     elif dd < 151: dd = 15
     elif dd < 161: dd = 16
     elif dd < 171: dd = 17
     elif dd < 181: dd = 18
     elif dd < 191: dd = 19
     elif dd < 201: dd = 20
     elif dd < 211: dd = 21
     elif dd < 221: dd = 22
     elif dd < 231: dd = 23
     elif dd < 241: dd = 24
     elif dd < 251: dd = 25
     elif dd < 261: dd = 26
     elif dd < 271: dd = 27
     elif dd < 281: dd = 28
     elif dd < 291: dd = 29
     elif dd < 301: dd = 30
     elif dd < 311: dd = 31
     elif dd < 321: dd = 32
     elif dd < 331: dd = 33
     elif dd < 341: dd = 34
     elif dd < 351: dd = 35
     else: dd = 36
     dd = str(dd)
     newFilename="ndvi_" + oldFilename[11:15] + "_dec_" + dd + ".tif"
     destination_path=os.path.join(outfolder, newFilename)
     os.rename(source_path, destination_path)

 

Model Builder -> Python script 2

 

Python script 3: Python script to calculate NDVI anomalies

def main():
    import arcpy
    import os
    arcpy.env.overwriteOutput = True

    # Check out SA license
    arcpy.CheckOutExtension("spatial")

    # settings
    ws_norm = r'D:\Modis250\NDVI_10dias_norm'
    ws_mean = r'D:\Modis250\NDVI_10dias_mean'
    ws_anom = r'D:\Modis250\NDVI_10dias_anom'

    # create a list of all raster in norm folder
    arcpy.env.workspace = ws_norm
    lst_ndvi_ras = arcpy.ListRasters()
     lst_ndvi_ras.sort()
    print ("lst_ndvi_ras: {}".format(lst_ndvi_ras))

    # get list of years
    decades = GetListOfDecades(lst_ndvi_ras)
     decades.sort()
    print ("decades: {}".format(decades))

    # loop through each year
    for decade in decades:
        # Get rasters for given year
        lst_ndvi_decade = GetListOfRasterForGivenDecade(lst_ndvi_ras, decade)
        print("lst_ndvi_decade: {}".format(lst_ndvi_decade))

        # calculate mean raster for decade
        lst_ndvi_decade_paths = [os.path.join(ws_norm, r) for r in lst_ndvi_decade]
        print("lst_ndvi_decade_paths: {}".format(lst_ndvi_decade_paths))
        mean_ras = arcpy.sa.CellStatistics(lst_ndvi_decade_paths, "MEAN", "DATA")

        # store mean ras for decade
        # mean_dec_d_.tif
        out_name_mean = "mean_dec_{}.tif".format(decade)  ###
        out_name_path_mean =  os.path.join(ws_mean, out_name_mean)
        print ("out_name_path_mean: {}".format(out_name_path_mean))
        mean_ras.save(out_name_path_mean)

        # loop through each raster in decade
        for ndvi_ras in lst_ndvi_decade:
            # calculate anom for each raster in decade
            # anom_2002_dec_1.tif = norm_ndvi_2002_dec_1.tif - mean_dec_1_.tif;
            print ("ndvi_ras: {}".format(ndvi_ras))
            # minus_ras = arcpy.sa.Minus(ndvi_ras, out_name_path_mean)
            out_name_path_norm = os.path.join(ws_norm, ndvi_ras)
            minus_ras = arcpy.Raster(out_name_path_mean) - arcpy.Raster(out_name_path_norm)

            # store anom raster
            out_name_anom = ndvi_ras.replace("norm_ndvi", "anom")
            out_name_path_anom = os.path.join(ws_anom, out_name_anom)
            minus_ras.save(out_name_path_anom)
            print ("out_name_path_anom: {}".format(out_name_path_anom))


    # return SA license
    arcpy.CheckInExtension("spatial")

def GetListOfDecades(lst):
    # norm_ndvi_yyyy_dec_d.tif
    lst1 = list(set([r.split("_")[4] for r in lst]))
    return [r.split('.')[0] for r in lst1]

def GetListOfRasterForGivenDecade(lst, decade):
    lst1 = []
    search = "_{}.tif".format(decade)
    for r in lst:
        if search in r.lower():
            lst1.append(r)
    return lst1


if __name__ == '__main__':
    main()

Outcomes