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]

Luiz  ... could you format your code Code Formatting the basics please

And the error message? 

The last line?

#print ("Decade "+d+" of "+y+" ok")

# ---- maybe try this instead, it handles type conversion to string

print("Decade {} of {} ok".format(d, y))‍‍‍‍‍

[LuizVianna]

Thank you Dan!

The problem is in line "NDVIAnom=arcpy.gp.Minus_3d(NDVINormY, NDVIMean, NDVIAnomOut) #Calculate the last year decendial NDVI anomaly and record the output tifs"

I´ve also tried "NDVIAnom = NDVINormY - NDVIMean
NDVIAnom.save(NDVIAnomOut)"

The errors are:

For the arcpy.gp.Minus_3d approach:

Traceback (most recent call last):
File "<string>", line 22, 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))
RuntimeError: Object: Error in executing tool

For the "NDVIAnom = NDVINormY - NDVIMean

NDVIAnom.save(NDVIAnomOut)" approach

Traceback (most recent call last):
File "<string>", line 23, in <module>
TypeError: unsupported operand type(s) for -: 'list' and 'geoprocessing server result object'

[JoeBorgione]

Luiz- I've reformatted your script as Dan suggested to make it a little easier to read; in this case you say your problem is on line 33.  The error indicates your arguments are in error.  Looking upstream, my guess is Line 30, and your definition of NDVIAnomOut.  Try picking out the individual components that make of the string and set them to variables rather than stringing the whole thing inside quotes: that's what I  do and then using Dan's approach of '{}'.format(variable_name) works quite nicely...

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

Thank you Joe!

I´m not expert in Python and I´m just copying some codes from internet and trying to make them work.

I changed lines 24 and line 30 to:

24   NDVIMeanOut = ("D:\\Modis250\\NDVI_10dias_anom\\mean_dec_{}.tif".format(d))

30   NDVIAnomOut = ("D:\\Modis250\\NDVI_10dias_anom\\anom_{}_dec_{}.tif".format(y, d))

and when I run the code, after line 33

33   NDVIAnom=arcpy.gp.Minus_3d(NDVINormY, NDVIMean, NDVIAnomOut)

the error is:

Traceback (most recent call last):
File "<string>", line 24, 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))
RuntimeError: Object: Error in executing tool

If I try:

NDVIAnom = NDVINormY - NDVIMean
NDVIAnom.save(NDVIAnomOut)

the error is

Traceback (most recent call last):
File "<string>", line 25, in <module>
TypeError: unsupported operand type(s) for -: 'list' and 'geoprocessing server result object'

New code is: (sorry, I dont know how to format it)

# 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(n)
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)
NDVIMeanOut = ("D:\\Modis250\\NDVI_10dias_anom\\mean_dec_{}.tif".format(d))
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
NDVIAnomOut = ("D:\\Modis250\\NDVI_10dias_anom\\anom_{}_dec_{}.tif".format(y, d))
#NDVIAnom=arcpy.gp.Minus_3d(NDVINormY, NDVIMean, NDVIAnomOut) #Calculate the last year decendial NDVI anomaly and record the output tifs
NDVIAnom = NDVINormY - NDVIMean
NDVIAnom.save(NDVIAnomOut)
#print ("Decade "+d+" of "+y+" ok")
print("Decade {} of {} ok".format(d, y))

[XanderBakker]

Hi Joe Borgione ,

What does not  make sense is:

• There is for loop on line 15 that would require indentation afterwards.
• During the loop a variable "n" is set that is not used, causing a loop of 37 iteration doing the same thing
• the Minus_3d has the output set as parameter "Minus_3d (in_raster_or_constant1, in_raster_or_constant2, out_raster)" and this should not be used in combination as "another_output = Minus_3d (in_raster_or_constant1, in_raster_or_constant2, out_raster)"
• The first parameter of the Minus_3d is a  list of rasters and should be a single raster

[JoeBorgione]

Good catch Xander!  I've indented the entire loop. 

Luiz- as Xander mentions, your loop has a malfunction in it as you don't mention the the incrementer (n) to move through the loop.

Xander Bakker‌

Luiz Vianna‌

for n in range(1,37):
  print '{}'.format(n)‍‍

[XanderBakker]

I can imagine that the idea is to calculate for each raster in the list with "decendial NDVI rasters for the last year (y)" the difference needs to be calculated with the mean. The mean should be calculated outside the look and the loop should go through each item in the list of rasters.

Look at the following structure:

year_range = range(2004, 2011) # to get the values from 2004 until 2010
for year in year_range:
    print year

    # for each year:
    # - create list of all NDVI for current year: lst_ndvi_ras
    # - calculate the mean of all NDVI for current year: mean_ndvi_year

    for ndvi_ras in lst_ndvi_ras:
        # calculate anomaly and save result with unique name‍‍‍‍‍‍‍‍‍‍

I assume you will have a range of years. Each year will have a number of NDVI rasters (is that true?). How are the NDVI rasters for each year stored? Each year has its own folder or the year is part of the name? 

[LuizVianna]

Hi Xander!

Let me try to explain what I need: There are 423 NDVI tifs in folder 'D: \ Modis250 \ NDVI_10dias_norm'. The tifs are named by year (from 2002 to 2018) and "decade" (36 ten days periods) as: norm_ndvi_yyyy_dec_d.tif, where "yyyy" is the year and the "d" is the "decade".

The loop calculates the mean for each decade, considering all years:

Mean (norm_ndvi_2002_dec_1.tif, norm_ndvi_2003_dec_1.tif, ..., norm_ndvi_2018_dec_1.tif ).

The 36 means are storaged in 'D:\Modis250\NDVI_10dias_mean' like 'mean_dec_d_.tif'.

The loop also calculates the anomaly for a choosing year ("y" in the code), subtracting each decade NDVI from those year by the NDVI mean for each decade:

anom_2002_dec_1.tif = norm_ndvi_2002_dec_1.tif - mean_dec_1_.tif;

anom_2002_dec_2.tif = norm_ndvi_2002_dec_2.tif - mean_dec_2_.tif;

anom_yyyy_dec_d.tif = norm_ndvi_yyyy_dec_d.tif - mean_dec_d_.tif;

The anomaly results shoud be storages in 'D:\Modis250\NDVI_10dias_anom' like anom_yyyy_dec_d.tif

The code seems to be ok until the mean calculation, however the 'Minus' operation is not working.

How can I set the input rasters for the "Minus" operation?

Thank you!

[LuizVianna]

The 36 "decades" NDVI means could be temporary... i dont need to store them, just use to calculate the anomaly...