Latest Contributions by LuizVianna

‎06-15-2018

In line 19 is ". join(raster1)" or .format(raster1)?

‎06-15-2018

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

‎06-15-2018

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()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

‎06-14-2018

Script 1 working now! Thank you! Now, can you please help me with the script 2? Pre-processing raster datasets: (Python script 2) Project the renamed raster datasets in "X:Modis250\NDVI_10dias" as ("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 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)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

‎06-14-2018

Problem was in line 7: Originals -> Originais The script runns, however does´t overwrite the output files: If there are more than one tif per decade, the higher julian day should represent the decade. Traceback (most recent call last): File "<string>", line 37, in <module> File "<string>", line 22, in main FileExistsError: [WinError 183] Não é possível criar um arquivo já existente: 'D:\\Modis250\\Originais\\MODIS.ndvi.2002017.yL6000.BOKU.tif' -> 'D:\\Modis250\\NDVI_10dias\\ndvi_2002_dec_2.tif' Does the code considered the two first files as belonging to decade two?

‎06-14-2018

Same error for your code: Traceback (most recent call last): File "<string>", line 37, in <module> File "<string>", line 13, in main TypeError: 'NoneType' object is not iterable

‎06-14-2018

Error running at ArcGisPro: Traceback (most recent call last): File "<string>", line 40, in <module> File "<string>", line 13, in main TypeError: 'NoneType' object is not iterable in code: def main(): import arcpy import os arcpy.env.overwriteOutput = True ws_in = r'D:\Modis250\Originals' ws_out = r'D:\Modis250\NDVI_10dias' arcpy.env.workspace = ws_in lst_ras = arcpy.ListRasters("*.tif") for ras_name in sorted(lst_ras): in_name_path = os.path.join(ws_in, ras_name) j = int(ras_name[15:18]) d = GetDecadeFromJulian(j) out_name = "ndvi_{}_dec_{}.tif".format(j, d) out_name_path = os.path.join(ws_out, out_name) os.rename(in_name_path, out_name_path) #def GetDecadeFromJulian(j): #d = (j+1 - ((j-1) % 10)) / 10 + 1 #return str(d) if d < 37 else '36' def GetDecadeFromJulian(year, j): import datetime date = datetime.datetime(year, 1, 1) + datetime.timedelta(j) m = (date.month - 1) * 3 if date.day <= 10: d = 1 elif date.day <= 20: d = 2 else: d = 3 return m + d if __name__ == '__main__': main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

‎06-14-2018

Question: how often do you plan to run this script? Does the input data change?‍ The first process will be runned one time for all the past years files (from 2000 to today). After that, the script will be runned every ten days period (decade). The new input (original file "MODIS.ndvi.YYYYjjj.yL6000.BOKU.tif") will be pre-processed and the result will be added to the list of standardized and scaled (0-100) raster datasets, to calculate the anomaly for the last decade in current year. It is a continous process. As the name of the original file have the "YYYYjjj" format, is it possible to rename the files considering the leap years, according to the table below? Decade Date julian day in common year julian day in leap year 1 01/01 to 01/10 10 10 2 01/11 to 01/20 20 20 3 01/21 to 01/31 31 31 4 02/01 to 02/10 41 41 5 02/11 to 02/20 51 51 6 02/21 to 02/28-29 59 60 7 03/01 to 03/10 69 70 8 03/11 to 03/20 79 80 9 03/21 to 03/31 90 91 10 04/01 to 04/10 100 101 11 04/11 to 04/20 110 111 12 04/21 to 04/30 120 121 13 05/01 to 05/10 130 131 14 05/11 to 05/20 140 141 15 05/21 to 05/31 151 152 16 06/01 to 06/10 161 162 17 06/11 to 06/20 171 172 18 06/21 to 06/30 181 182 19 07/01 to 07/10 191 192 20 07/11 to 07/20 201 202 21 07/21 to 07/31 212 213 22 08/01 to 08/10 222 223 23 08/11 to 08/20 232 233 24 08/21 to 08/31 243 244 25 09/01 to 09/10 253 254 26 09/11 to 09/20 263 264 27 09/21 to 09/30 273 274 28 10/01 to 10/10 283 284 29 10/11 to 10/20 293 294 30 10/21 to 10/31 304 305 31 11/01 to 11/10 314 315 32 11/11 to 11/20 324 325 33 11/21 to 11/30 334 335 34 12/01 to 12/10 344 345 35 12/11 to 12/20 354 355 36 12/21 to 12/31 365 366

‎06-14-2018

The script 1 worked, however the output name is wrong: The names should be "ndvi_yyyy_dec_d.tif" where "yyyy" is the year (from the input file name (MODIS.ndvi.YYYYjjj.yL6000.BOKU.tif ) and "d" the decade. Why the ".2" after the decade number in the outputs? If there are more then one input file for decade, de code shoud overwrite the first.

‎06-14-2018

https://community.esri.com/message/778938-python-routine-to-calculate-modis-ndvi-anomaly

‎06-14-2018

Xander Bakker‌ 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()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

‎06-13-2018

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: Select and download NDVI filtered raster datasets; Organizing files and folders: (Python script 1) Store downloaded original raster datasets at X:\Modis250\Originals 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). Transfer renamed files from X:\Modis250\Originals to X:\Modis250\NDVI10dias Pre-processing raster datasets: (Python script 2) Project original raster datasets from Lat/Lon WGS 1984 to Lat/Lon SIRGAS2000 Overwrite spurious data: NDVI = (if NDVI <1, 0 and if NDVI >=9999, 9999) Standardize raster datasets: NDVI = (NDVI – MEAN(NDVI))/STD(NDVI) Scale raster datasets to 0 – 100: NDVI = ((NDVI-MIN(NDVI))/(MAX(NDVI)-MIN(NDVI)) Mask standardized raster datasets by the area of interest 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 = 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!

‎06-13-2018

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

‎06-13-2018

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()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

‎06-13-2018

My code: I´ve tried line 27 and lines 28-29. Both returned errors. # 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))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Online Status	Offline
Date Last Visited	‎11-11-2020 02:23 AM

My Ideas

Latest Contributions by LuizVianna

Re: Python routine to calculate MODIS NDVI anomaly

Re: Python routine to calculate MODIS NDVI anomaly

Re: Python routine to calculate MODIS NDVI anomaly

Re: Python routine to calculate MODIS NDVI anomaly

Re: Python routine to calculate MODIS NDVI anomaly

Re: Python routine to calculate MODIS NDVI anomaly

Re: Python routine to calculate MODIS NDVI anomaly

Re: Python routine to calculate MODIS NDVI anomaly

Re: Python routine to calculate MODIS NDVI anomaly

Re: Python script to calculate NDVI anomalies.

Python routine to calculate MODIS NDVI anomaly

Re: Python script to calculate NDVI anomalies.

Re: Python script to calculate NDVI anomalies.

Re: Python script to calculate NDVI anomalies.

Re: Python script to calculate NDVI anomalies.

Re: How to calculate multiple equations (a+bX+cY+d...

Re: Python routine to calculate MODIS NDVI anomaly

How to calculate multiple equations (a+bX+cY+dZ), ...

Re: Python routine to calculate MODIS NDVI anomaly

Re: Python script to calculate NDVI anomalies.