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

[DanPatterson_Retired]

I just assembled your sample to print out the variables.

There is an error as you can see in the names

workspace = "D:\\Modis250\\NDVI_10dias_norm" #Workspace with decendial (ten days period) NDVI
n = 1
d='{}'.format(n)
y = str(2004)
list_ndvi_ras = "norm_ndvi_*_"+d+".tif", "tif"
list_ndvi_year = "norm_ndvi_"+y+"_dec_"+d+".tif", "tif"
ndvi_year = workspace + str(list_ndvi_year)
NDVIMeanOut = ("D:\\Modis250\\NDVI_10dias_anom\\mean_dec_{}.tif".format(d))



workspace
'D:\\Modis250\\NDVI_10dias_norm'

list_ndvi_ras
('norm_ndvi_*_1.tif', 'tif')

list_ndvi_year
('norm_ndvi_2004_dec_1.tif', 'tif')

ndvi_year
"D:\\Modis250\\NDVI_10dias_norm('norm_ndvi_2004_dec_1.tif', 'tif')"‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[LuizVianna]

The problem is occurring in the "Minus" operation, the same as occurred in my initial code.

[LuizVianna]

My initial code...

# Import arcpy module
import arcpy
import os
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(n)
d='{}'.format(n)
y = str(2004)
list_ndvi_ras=arcpy.ListRasters("norm_ndvi_*_"+d+".tif", "tif")#List of decendial NDVI for all years to calculate the dcendial mean
list_ndvi_year=arcpy.ListRasters("norm_ndvi_"+y+"_dec_"+d+".tif", "tif")#List of decendial NDVI rasters for the last year (y)
ndvi_year = arcpy.env.workspace+str(list_ndvi_year)
#NDVIMeanOut = "D:\\Modis250\\NDVI_10dias_mean\\"+"mean_dec_"+d+".tif" #Define the decendial NDVI mean output (could be temporary)
NDVIMeanOut = ("D:\\Modis250\\NDVI_10dias_anom\\mean_dec_{}.tif".format(d))
NDVIMean = arcpy.gp.CellStatistics_sa(list_ndvi_ras, 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
NDVIAnomOut = ("D:\\Modis250\\NDVI_10dias_anom\\anom_{}_dec_{}.tif".format(y, d))
NDVIAnom=arcpy.gp.Minus_3d(ndvi_year, NDVIMean, NDVIAnomOut) #Calculate the last year decendial NDVI anomaly and record the output tifs
#NDVIAnom = ndvi_year - NDVIMean
#NDVIAnom.save(NDVIAnomOut)
#print ("Decade "+d+" of "+y+" ok")
print("Decade {} of {} ok".format(d, y))

The error:

Traceback (most recent call last):
File "<string>", line 27, in <module>
File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing\_base.py", line 506, in <lambda>
return lambda *args: val(*gp_fixargs(args, True))
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000865: Input raster or constant value 1: D:\Modis250\NDVI_10dias_norm['norm_ndvi_2004_dec_1.tif'] does not exist.
Failed to execute (Minus).

[LuizVianna]

The code below runs perfectly at ArcMap 10.5, however it is not working at ArcGis Pro.

At ArcGisPro,

It prints the "lst_ndvi_ras";

prints the decades (unsorted);

prints the "lst_ndvi_decade";

prints the "out_name_path_mean"

prints the "ndvi_ras"

stores the first tif file of mean,

and present the following error:

Traceback (most recent call last):
File "<string>", line 71, in <module>
File "<string>", line 44, in main
File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Functions.py", line 4271, in Minus
in_raster_or_constant2)
File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Utils.py", line 53, in swapper
result = wrapper(*args, **kwargs)
File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Functions.py", line 4268, in Wrapper
["Minus", in_raster_or_constant1, in_raster_or_constant2])

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()
    print ("lst_ndvi_ras:", lst_ndvi_ras)

    # get list of years
    decades = GetListOfDecades(lst_ndvi_ras)
    print ("decades:", 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:", lst_ndvi_decade)

        # calculate mean raster for decade
        mean_ras = arcpy.sa.CellStatistics(lst_ndvi_decade, "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:", 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:", ndvi_ras)
            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)
            print ("out_name_path_anom:", 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()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[XanderBakker]

So to be clear... did you try the code I provided or did you use parts of it?

[LuizVianna]

The code you provided. Exactly as above. I only inserted the "()" at lines with "print".

[XanderBakker · Accepted solution]

I tried it too and had some different results in Pro (including ArcGIS Pro closing very efficiently). I also obtained an error with the Cell statistics in some cases and decided to remove the minus and include the paths in the list of rasters that are used in the Cell Statistics. This seems to be working (you will need to change the paths back again).

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'
    ws_norm = r'C:\GeoNet\Modis250\NDVI_10dias_norm'
    ws_mean = r'C:\GeoNet\Modis250\NDVI_10dias_mean'
    ws_anom = r'C:\GeoNet\Modis250\NDVI_10dias_anom'


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

    # get list of years
    decades = GetListOfDecades(lst_ndvi_ras)
    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()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[JoeBorgione]

...had some different results in Pro (including ArcGIS Pro closing very efficiently)....

Ah haa.... the efficient closing feature....

[LuizVianna]

Now it worked! Thank you!

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

Summary workflow:

1. Select and download NDVI filtered raster datasets;
   1. Organizing files and folders: (Python script 1)
   2. Store downloaded original raster datasets at X:\Modis250\Originals
   3. Rename original raster datasets to filename standard “ndvi_YYYY_dec_dd.tif”, where “YYYY” is the year and “dd” is the decade (ten days period).
   4. Transfer renamed files from X:\Modis250\Originals to X:\Modis250\NDVI10dias
2. Pre-processing raster datasets: (Python script 2)
   1. Project original raster datasets from Lat/Lon WGS 1984 to Lat/Lon SIRGAS2000
   2. Overwrite spurious data: NDVI = (if NDVI <1, 0 and if NDVI >=9999, 9999)
   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))
   5. Mask standardized raster datasets by the area of interest 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 = MINUS(ndvi_YYYY_dec_dd.tif, mean_dec_dd.tif)

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

Can we still discussing this on this topic or I need to start another one?

Best regards!

[XanderBakker]

Perhaps if it is a new question, it would be good to start a new thread, but it is always good to include a link to this thread (to understand what else is used in the process) and to tag me if you want me to help you with the other question (although I prefer python scripts over model builder).