AnsweredAssumed Answered

Python routine to calculate MODIS NDVI anomaly

Question asked by vianna on Jun 14, 2018
Latest reply on Jun 20, 2018 by xander_bakker

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