topic Re: Python routine to calculate MODIS NDVI anomaly in ArcGIS Pro Questions

Python routine to calculate MODIS NDVI anomaly

LuizVianna — Sat, 11 Dec 2021 12:17:55 GMT

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;
1. Store downloaded original raster datasets at X:\Modis250\Originais\
2. Organizing files and folders: (Python script 1)
3. 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).
4. Transfer renamed files from X:\Modis250\Originais to X:\Modis250\NDVI10dias
Pre-processing raster datasets: (Python script 2)
1. 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.
2. Overwrite spurious data values: NDVI = Con("%NDVI%"<1,0,Con("%NDVI%">9998,9999,"%NDVI%"))
3. Standardize raster datasets: NDVI = (NDVI – Mean(NDVI))/Std(NDVI)
4. Scale raster datasets to 0 – 100: NDVI = ((NDVI-Min(NDVI))/(Max(NDVI)-Min(NDVI)) * 100
5. 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)
1. 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,…)
2. 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()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Python routine to calculate MODIS NDVI anomaly

XanderBakker — Sat, 11 Dec 2021 12:17:58 GMT

Let's try and do this in parts.

Question 1: What is the best option: to run separated scripts or create a unique code?

You can work with a single script that does all 3 steps, but it might be better to break it up in parts(as you did). This way you can check the intermediate results and run parts separately and also reuse parts for other workflows.

Question 2: How to translate the model builder to python script (Python script 2)?

You can export the Model to a Python script, but that probably results in some really ugly code. It is better to rewrite the model into a new script. This also allows for some additional checks that might not be easy to do inside a model.

Question 3: If the best option is to run separated, how to make it in a logical sequence (script 1 -> script 2 -> script 3)

Depends. As I mentioned above, it may be easier to run parts separately and do some (manual) checks between the steps to be sure things go right. When you have tested thoroughly, you can call the scripts from a "master" script and run the entire process. Question: how often do you plan to run this script? Does the input data change?

I think that the first script could be rewritten like this (untested):

def main():
    import arcpy
    import os

    arcpy.env.overwriteOutput = True

    ws_in = 'O:\Modis250\Originals'
    ws_out = 'O:\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'

if __name__ == '__main__':
    main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Python routine to calculate MODIS NDVI anomaly

LuizVianna — Thu, 14 Jun 2018 14:47:38 GMT

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.

Re: Python routine to calculate MODIS NDVI anomaly

LuizVianna — Thu, 14 Jun 2018 17:41:25 GMT

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

Re: Python routine to calculate MODIS NDVI anomaly

XanderBakker — Sat, 11 Dec 2021 12:18:00 GMT

This would require changing the function to get the decade like this:

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

Both year and Julian (j) should be provided as numeric value.

Re: Python routine to calculate MODIS NDVI anomaly

LuizVianna — Sat, 11 Dec 2021 12:18:03 GMT

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

Re: Python routine to calculate MODIS NDVI anomaly

XanderBakker — Sat, 11 Dec 2021 12:18:06 GMT

It seems you are missing some code:

def main():
    import arcpy
    import os

    arcpy.env.overwriteOutput = True

    ws_in = 'O:\Modis250\Originals'
    ws_out = 'O:\Modis250\NDVI_10dias'

    arcpy.env.workspace = ws_in
    lst_ras = arcpy.ListRasters("*.tif")

    for ras_name in sorted(lst_ras):
        # MODIS.ndvi.YYYYjjj.yL6000.BOKU.tif
         in_name_path = os.path.join(ws_in, ras_name)
         y = ras_name[11:15]
         j = int(ras_name[15:18])
         d = GetDecadeFromJulian(int(y), j)

         out_name = "ndvi_{}_dec_{}.tif".format(y, d)
         out_name_path = os.path.join(ws_out, out_name)
         os.rename(in_name_path, out_name_path)

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

Re: Python routine to calculate MODIS NDVI anomaly

LuizVianna — Thu, 14 Jun 2018 19:07:49 GMT

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

Re: Python routine to calculate MODIS NDVI anomaly

LuizVianna — Thu, 14 Jun 2018 19:13:08 GMT

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?

Re: Python routine to calculate MODIS NDVI anomaly

XanderBakker — Thu, 14 Jun 2018 19:19:11 GMT

You're using "os.rename" and on Windows if the destination exists, it will throw an error.

Re: Python routine to calculate MODIS NDVI anomaly

XanderBakker — Thu, 14 Jun 2018 19:21:15 GMT

Perhaps "shutil.move()" could do the trick. See python - Force Overwrite in Os.Rename - Stack Overflow

Re: Python routine to calculate MODIS NDVI anomaly

LuizVianna — Sat, 11 Dec 2021 12:18:09 GMT

Script 1 working now! Thank you!

Now, can you please help me with the script 2?

Pre-processing raster datasets: (Python script 2)
1. 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.
2. Overwrite spurious data values: NDVI = Con("%NDVI%"<1,0,Con("%NDVI%">9998,9999,"%NDVI%"))
3. Standardize raster datasets: NDVI = (NDVI – Mean(NDVI))/Std(NDVI)
4. Scale raster datasets to 0 – 100: NDVI = ((NDVI-Min(NDVI))/(Max(NDVI)-Min(NDVI)) * 100
5. 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)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Python routine to calculate MODIS NDVI anomaly

XanderBakker — Fri, 15 Jun 2018 12:31:51 GMT

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.

Re: Python routine to calculate MODIS NDVI anomaly

LuizVianna — Sat, 11 Dec 2021 12:18:12 GMT

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

Re: Python routine to calculate MODIS NDVI anomaly

XanderBakker — Sat, 11 Dec 2021 12:18:15 GMT

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

Re: Python routine to calculate MODIS NDVI anomaly

LuizVianna — Fri, 15 Jun 2018 18:26:04 GMT

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

Re: Python routine to calculate MODIS NDVI anomaly

LuizVianna — Fri, 15 Jun 2018 18:32:25 GMT

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

Re: Python routine to calculate MODIS NDVI anomaly

XanderBakker — Fri, 15 Jun 2018 18:57:41 GMT

Good catch, that should be format as you mentioned.

Re: Python routine to calculate MODIS NDVI anomaly

LuizVianna — Fri, 15 Jun 2018 19:22:10 GMT

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

Re: Python routine to calculate MODIS NDVI anomaly

LuizVianna — Fri, 15 Jun 2018 20:15:44 GMT

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>