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:
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()
Solved! Go to Solution.
Perhaps "shutil.move()" could do the trick. See python - Force Overwrite in Os.Rename - Stack Overflow
Script 1 working now! Thank you!
Now, can you please help me with the script 2?
This is the script from Model Builder:
# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# Script2_Model_Builder.py
# Created on: 2018-06-14 17:09:24.00000
# (generated by ArcGIS/ModelBuilder)
# Usage: Script2_Model_Builder <InFolder> <Wildcard> <Name> <Mask> <OutFolder1> <OutFolder2> <Raster_Format>
# Description:
# ---------------------------------------------------------------------------
# Import arcpy module
import arcpy
# Load required toolboxes
arcpy.ImportToolbox("Model Functions")
# Script arguments
InFolder = arcpy.GetParameterAsText(0)
if InFolder == '#' or not InFolder:
InFolder = "O:\\Modis250\\NDVI_10dias" # provide a default value if unspecified
Wildcard = arcpy.GetParameterAsText(1)
Name = arcpy.GetParameterAsText(2)
if Name == '#' or not Name:
Name = "ndvi_2001_dec_10" # provide a default value if unspecified
Mask = arcpy.GetParameterAsText(3)
if Mask == '#' or not Mask:
Mask = "O:\\Modis250\\SRTM\\mde_sc_90.tif" # provide a default value if unspecified
OutFolder1 = arcpy.GetParameterAsText(4)
if OutFolder1 == '#' or not OutFolder1:
OutFolder1 = "O:\\Modis250\\Temp\\mask_%Name%" # provide a default value if unspecified
OutFolder2 = arcpy.GetParameterAsText(5)
if OutFolder2 == '#' or not OutFolder2:
OutFolder2 = "O:\\Modis250\\NDVI_10dias_norm\\norm_%Name%" # provide a default value if unspecified
Raster_Format = arcpy.GetParameterAsText(6)
if Raster_Format == '#' or not Raster_Format:
Raster_Format = "TIF" # provide a default value if unspecified
# Local variables:
Recursive = "false"
Raster = "O:\\Modis250\\NDVI_10dias\\ndvi_2001_dec_10.tif"
NDVI_Project = "O:\\Modis250\\Temp\\ndvi_2001_dec_10_ProjectRast1.tif"
mean = "7770.0395015223"
std = "924.13000791157"
RDS_STD = "O:\\Modis250\\Temp\\%Name%_std.tif"
min = "-8.4079504013062"
max = "2.4119555950165"
# Set Geoprocessing environments
arcpy.env.scratchWorkspace = "O:\\Modis250\\Temp"
arcpy.env.workspace = "O:\\Modis250\\MODIS_13Q1_V6_filtered"
# Process: Iterate Rasters
arcpy.IterateRasters_mb(InFolder, Wildcard, Raster_Format, Recursive)
# Process: Project Raster
arcpy.ProjectRaster_management(Raster, NDVI_Project, "GEOGCS['GCS_SIRGAS_2000',DATUM['D_SIRGAS_2000',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]", "NEAREST", "1.8000795216263E-03 1.80029836029412E-03", "'SAD_1969_To_WGS_1984_14 + SAD_1969_To_SIRGAS_2000_1'", "", "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
# Process: Raster Calculator
tempEnvironment0 = arcpy.env.mask
arcpy.env.mask = Mask
arcpy.gp.RasterCalculator_sa("Con(\"%NDVI_Project%\"<1,0,Con(\"%NDVI_Project%\">9998,9999,\"%NDVI_Project%\"))", OutFolder1)
arcpy.env.mask = tempEnvironment0
# Process: Get Raster Properties
arcpy.GetRasterProperties_management(OutFolder1, "MEAN", "")
# Process: Get Raster Properties (2)
arcpy.GetRasterProperties_management(OutFolder1, "STD", "")
# Process: Raster Calculator (2)
arcpy.gp.RasterCalculator_sa("(\"%OutFolder1%\" - float(%mean%)) / float(%std%)", RDS_STD)
# Process: Get Raster Properties (4)
arcpy.GetRasterProperties_management(RDS_STD, "MINIMUM", "")
# Process: Get Raster Properties (3)
arcpy.GetRasterProperties_management(RDS_STD, "MAXIMUM", "")
# Process: Raster Calculator (3)
arcpy.gp.RasterCalculator_sa("(\"%RDS_STD%\" - float(%min%)) / (float(%max%) - float(%min%)) * 100", OutFolder2)
I notice in your description that you apply the mask after the min, max, mean and std have been defined. These values will be different before and after applying the mask. Is using the values before applying the mask correct?
You may have noticed that the exported script contains some "anomalies". The use of an iteration in model builder normally does not get exported correctly into the output script. As you can see there are local variables (mean, std, min, max, even the last two overwrite exiting functions) that are defined before the iteration and would be used inside the iteration (if the code would run). Don't get me wrong, in Model Builder it will work correctly, but as script it will not work.
So this requires a different approach. As mentioned in other discussions, if normally takes more time to understand a model builder exported script and adjust it then it will take to start from scratch.
I notice in your description that you apply the mask after the min, max, mean and std have been defined. These values will be different before and after applying the mask. Is using the values before applying the mask correct?
You are right! I´m trying to write a new code to apply the mask before the statistics and after projecting, but some issues are happening.
Te error:
Traceback (most recent call last):
File "<string>", line 48, in <module>
File "<string>", line 26, in main
NameError: name 'ExtractByMask' is not defined
The line 4 is not running also, is it important for the ExtractByMask?
The new code:
def main():
import arcpy
from arcpy import env
#from arcpy.sa import *
import os
arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True
# Set Geoprocessing environments
ws_scr = r'D:\Modis250\Temp' #Scratch workspace
ws_in = r'D:\Modis250\NDVI_10dias' #Input workspace
ws_out = r'D:\Modis250\NDVI_10dias_norm' #Output workspace
mask=r'D:\Modis250\SRTM\mde_sc_90.tif' #Raster used as mask
arcpy.env.scratchWorkspace = ws_scr
arcpy.env.workspace = ws_in
lst_ras = arcpy.ListRasters("*.tif")
# Reproject original rasters
for raster1 in lst_ras:
print ("Processing: " + raster1)
outRaster1 = ws_scr+"\\rep_"+raster1 #Temporary file
arcpy.ProjectRaster_management(raster1, outRaster1,"GEOGCS['GCS_SIRGAS_2000',DATUM['D_SIRGAS_2000',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]", "NEAREST", "1.8000795216263E-03 1.80029836029412E-03", "'SAD_1969_To_WGS_1984_14 + SAD_1969_To_SIRGAS_2000_1'", "", "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
print (outRaster1 + " projected")
#outRaster2 = ws_scr+"\\mask_"+raster1
outRaster2 = ExtractByMask(outRaster1, mask) #Temporary file
outRaster2.save(ws_scr+"\\mask_"+raster1)
print (outRaster2 + " masked")
#outRaster3 = ws_scr+"\\aju_"+raster1
outRaster3 = Con(outRaster2<1,0,Con(outRaster2>9998,9999,outRaster2))#Temporary file
outRaster3.save(ws_scr+"\\aju_"+raster1)
print (outRaster3 + " ajusted")
mean = arcpy.GetRasterProperties_management(outRaster3, "MEAN")
std = arcpy.GetRasterProperties_management(outRaster3, "STD")
outRaster4 = ((outRaster3-mean)/std)#Temporary file
outRaster4.save(ws_scr+"\\std_"+raster1)
print (outRaster4 + " standardized")
max = arcpy.GetRasterProperties_management(outRaster4, "MAXIMUM")
min = arcpy.GetRasterProperties_management(outRaster4, "MINIMUM")
outRaster = ((outRaster4-min)/(max-min))#Final result
outRaster.save(ws_out+"\\"+"norm_"+raster1)
print (outRaster + " processed")
if __name__ == '__main__':
main()
Some small adjustments. Can you see where this produces an error and if anything is printed?
def main():
import arcpy
import os
arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True
# Set Geoprocessing environments
ws_scr = r'D:\Modis250\Temp' #Scratch workspace
ws_in = r'D:\Modis250\NDVI_10dias' #Input workspace
ws_out = r'D:\Modis250\NDVI_10dias_norm' #Output workspace
mask = r'D:\Modis250\SRTM\mde_sc_90.tif' #Raster used as mask
arcpy.env.scratchWorkspace = ws_scr
arcpy.env.workspace = ws_in
lst_ras = arcpy.ListRasters("*.tif")
# Reproject original rasters
for raster1 in lst_ras:
print ("Processing: " + raster1)
outRaster1 = os.path.join(ws_scr, "rep_{}".join(raster1)) # Temporary file
arcpy.ProjectRaster_management(raster1, outRaster1,"GEOGCS['GCS_SIRGAS_2000',DATUM['D_SIRGAS_2000',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]", "NEAREST", "1.8000795216263E-03 1.80029836029412E-03", "'SAD_1969_To_WGS_1984_14 + SAD_1969_To_SIRGAS_2000_1'", "", "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
print (outRaster1 + " projected")
#outRaster2 = ws_scr+"\\mask_"+raster1
outRaster2 = arcpy.sa.ExtractByMask(outRaster1, mask) # Temporary file
outRaster2.save(os.path.join(ws_scr, "mask_{}".format(raster1)))
print (outRaster2 + " masked")
#outRaster3 = ws_scr+"\\aju_"+raster1
outRaster3 = arcpy.sa.Con(outRaster2 < 1, 0, arcpy.sa.Con(outRaster2 > 9998, 9999, outRaster2)) # Temporary file
outRaster3.save(os.path.join(ws_scr, "aju_{}".format(raster1)))
print (outRaster3 + " ajusted")
mean = arcpy.GetRasterProperties_management(outRaster3, "MEAN")
std = arcpy.GetRasterProperties_management(outRaster3, "STD")
outRaster4 = ((outRaster3-mean)/std) # Temporary file
outRaster4.save(os.path.join(ws_scr, "std_{}".format(raster1)))
print (outRaster4 + " standardized")
ras_max = arcpy.GetRasterProperties_management(outRaster4, "MAXIMUM")
ras_min = arcpy.GetRasterProperties_management(outRaster4, "MINIMUM")
outRaster = ((outRaster4-ras_min)/(ras_max-ras_min)) # Final result
outRaster.save(os.path.join(ws_out, "norm_{}".format(raster1)))
print (outRaster + " processed")
arcpy.CheckInExtension("spatial")
if __name__ == '__main__':
main()
In my code, the ProjectRaster_management worked and the projected raster dataset(outRaster1) was generated. In the next step, the error occoured at "ExtractByMask".
In your code, neither the "project" step worked. The error was:
Traceback (most recent call last):
File "<string>", line 48, in <module>
File "<string>", line 20, in main
File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 9090, in ProjectRaster
raise e
File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 9087, in ProjectRaster
retval = convertArcObjectToPythonObject(gp.ProjectRaster_management(*gp_fixargs((in_raster, out_raster, out_coor_system, resampling_type, cell_size, geographic_transform, Registration_Point, in_coor_system), True)))
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 000445: Extension is invalid for the output raster format.
Failed to execute (ProjectRaster).
In line 19 is ".join(raster1)" or .format(raster1)?
Good catch, that should be format as you mentioned.
I replaced the line 19 and executed the code in both ArcGis 10.5 and ArcGis Pro. Before execute the code, I deleted all Temp files. Follow the prints and errors:
ArcGis 10.5: After run the "project" step, the system adds the rep_ndvi_2002_dec_10.tif to the ArcMap´s TOC. In the Temp workspace, two tifs were generated: "rep_ndvi_2002_dec_10.tif" and "mask_ndvi_2002_dec_10.tif", so I think that the "ExtractByMask" is working.
Processing: ndvi_2002_dec_10.tif
D:\Modis250\Temp\rep_ndvi_2002_dec_10.tif projected
Runtime error
Traceback (most recent call last):
File "<string>", line 48, in <module>
File "<string>", line 26, in main
File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Functions.py", line 4329, in Plus
in_raster_or_constant2)
File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Utils.py", line 53, in swapper
result = wrapper(*args, **kwargs)
File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Functions.py", line 4326, in Wrapper
["Plus", in_raster_or_constant1, in_raster_or_constant2])
RuntimeError: ERROR 000732: Input Raster: Dataset masked does not exist or is not supported
ArcGisPro: In the Temp workspace, two tifs were generated: "rep_ndvi_2002_dec_10.tif" and "mask_ndvi_2002_dec_10.tif", so I think that the "ExtractByMask" is working.
Processing: ndvi_2002_dec_10.tif
D:\Modis250\Temp\rep_ndvi_2002_dec_10.tif projected
Traceback (most recent call last):
File "<string>", line 48, in <module>
File "<string>", line 26, in main
File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Functions.py", line 4353, in Plus
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 4350, in Wrapper
["Plus", in_raster_or_constant1, in_raster_or_constant2])
RuntimeError: ERROR 000732: Input Raster: Dataset masked does not exist or is not supported
Running from prompt, same error:
C:\Program Files\ArcGIS\Pro\bin\Python\Scripts>propy.bat D:\Modis250\Python\Script2_2.py
Processing: ndvi_2002_dec_10.tif
D:\Modis250\Temp\rep_ndvi_2002_dec_10.tif projected
Traceback (most recent call last):
File "D:\Modis250\Python\Script2_2.py", line 48, in <module>
main()
File "D:\Modis250\Python\Script2_2.py", line 26, in main
print (outRaster2 + " masked")
File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\sa\Functions.py", line 4353, in Plus
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 4350, in Wrapper
["Plus", in_raster_or_constant1, in_raster_or_constant2])
RuntimeError: ERROR 000732: Input Raster: Dataset masked does not exist or is not supported
C:\Program Files\ArcGIS\Pro\bin\Python\Scripts>