by
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:
- Select and download NDVI filtered raster datasets;
- Store downloaded original raster datasets at X:\Modis250\Originais\
- Organizing files and folders: (Python script 1)
- Rename original raster datasets (MODIS.ndvi.YYYYjjj.yL6000.BOKU.tif) to filename standard “ndvi_YYYY_dec_dd.tif”, where “YYYY” is the year, “jjj” is the Julian day and “dd” is the decade (ten days period).
- Transfer renamed files from X:\Modis250\Originais to X:\Modis250\NDVI10dias
- Pre-processing raster datasets: (Python script 2)
- Project renamed raster datasets ("ndvi_YYYY_dec_dd.tif") from Lat/Lon WGS 1984 to Lat/Lon SIRGAS2000 using "SAD_1969_To_WGS_1984_14 + SAD_1969_To_SIRGAS_2000_1" geographic transformation.
- Overwrite spurious data values: NDVI = Con("%NDVI%"<1,0,Con("%NDVI%">9998,9999,"%NDVI%"))
- Standardize raster datasets: NDVI = (NDVI – Mean(NDVI))/Std(NDVI)
- Scale raster datasets to 0 – 100: NDVI = ((NDVI-Min(NDVI))/(Max(NDVI)-Min(NDVI)) * 100
- Mask standardized raster datasets by the area of interest (X:\Modis250\SRTM\mde_sc_90.tif) and output to X:\Modis250\NDVI_10dias_norm
- Processing data (Python script 3)
- 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,…)
- Calculate NDVI anomaly for each year by decade: anom_yyyy_dec_dd.tif = ndvi_YYYY_dec_dd.tif - mean_dec_dd.tif
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()
