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")
Solved! Go to Solution.
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')"
The problem is occurring in the "Minus" operation, the same as occurred in my initial code.
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).
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()
So to be clear... did you try the code I provided or did you use parts of it?
The code you provided. Exactly as above. I only inserted the "()" at lines with "print".
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()
...had some different results in Pro (including ArcGIS Pro closing very efficiently)....
Ah haa.... the efficient closing feature....
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:
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!
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).