Python script to calculate NDVI anomalies.

[LuizVianna]

Why the last line of this code is not working? It is the line of "Minus_3d" calculation?

This is a python script to calculate NDVI anomalies.

Thanks.

# Import arcpy module
import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

#Define workspaces
arcpy.env.workspace="D:\\Modis250\\NDVI_10dias_norm" #Workspace with decendial (ten days period) NDVI
#Filenames standard is "norm_ndvi_yyyy_dec_d.tif" where "yyyy" is the year and "d" is the ten days period: d=1 = january 1-10; d=2 = january 11-20; ...

for n in range (1,37):
d=str(1)
y = str(2004)
NDVINormAll=arcpy.ListRasters("norm_ndvi_*_"+d+".tif", "tif")#List of decendial NDVI for all years to calculate the dcendial mean
NDVINormY=arcpy.ListRasters("norm_ndvi_"+y+"_dec_"+d+".tif", "tif")#List of decendial NDVI rasters for the last year (y)
NDVIMeanOut = "D:\\Modis250\\NDVI_10dias_mean\\"+"mean_dec_"+d+".tif" #Define the decendial NDVI mean output (could be temporary)
NDVIMean = arcpy.gp.CellStatistics_sa(NDVINormAll, NDVIMeanOut,"MEAN", "NODATA")#Calculate the decendial NDVI mean (could be temporary)
NDVIAnomOut = "D:\\Modis250\\NDVI_10dias_anom\\"+"anom_"+y+"_dec_"+d+".tif" #Define the last year decendial NDVI anomaly output
NDVIAnom=arcpy.gp.Minus_3d(NDVINormY, NDVIMean, NDVIAnomOut) #Calculate the last year decendial NDVI anomaly and record the output tifs
print ("Decade "+d+" of "+y+" ok")

[LuizVianna]

https://community.esri.com/message/778938-python-routine-to-calculate-modis-ndvi-anomaly 

[LuizVianna]

Thank you, Xander!

Unfortunatly I did´n understood this code, so I´ve tried to run the script step by step. However, I dont know how. I copied and paste fom line 1 to 16, executed the code and then I typed "ws_mean", the error was:

Traceback (most recent call last):
File "<string>", line 1, in <module>
NameError: name 'ws_mean' is not defined

I expected for 'D:\Modis250\NDVI_10dias_mean'

If I copy all the script and run, I had this error:

Traceback (most recent call last):
File "<string>", line 65, in <module>
File "<string>", line 27, in main
NameError: name 'lst_ndvi_year' is not defined

I´m working on ArcGis Pro 2.1.2, an I think my Python is 2.7 (I also have ArcGis 10.5).

[XanderBakker]

You shouldn't run this by parts, since the code is within the "def main()" and "if __name__ == '__main__':" the variable will not be set.

You could try to run this by parts, without the outer parts

[XanderBakker]

Something like this:

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

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()

# get list of years
decades = GetListOfDecades(lst_ndvi_ras)

# loop through each year
for decade in decades:
    # Get rasters for given year
    lst_ndvi_decade = GetListOfRasterForGivenDecade(lst_ndvi_ras, decade)

    # calculate mean raster for decade
    mean_ras = arcpy.sa.CellStatistics(lst_ndvi_year, "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)
    mean_ras.save(out_name_path_mean)

    # loop through each raster in decade
    for ndvi_ras in lst_ndvi_ras:
        # calculate anom for each raster in decade
        # anom_2002_dec_1.tif = norm_ndvi_2002_dec_1.tif - mean_dec_1_.tif;
        minus_ras = arcpy.sa.Minus(ndvi_ras, out_name_path_mean)

        # 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)


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

... and/or include some print statements to see what is going on.