Latest Contributions by LuizVianna

Python routine to calculate MODIS NDVI anomaly: Xander Bakker‌ modis‌ndvi calculation ndvi 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 def main():      import arcpy      import os      import shutil      arcpy.env.overwriteOutput = True      ws_in = r'D:\Modis250\Originais'      ws_out = r'D:\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)           shutil.move(in_name_path, out_name_path)           print (out_name + " renamed") def GetDecadeFromJulian(year, j):      import datetime      date = datetime.datetime(year-1, 12, 31) + 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()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ 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 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")      lst_ras.sort()     # Pre-processing original rasters      for raster1 in lst_ras:           print ("Processing: " + raster1)                   #outRaster1= Original raster projected from CGS-WGS1964 to CGS-SIRGAS-2000                    outRaster1 = os.path.join(ws_scr, "rep_{}".format(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 ("rep_{}".format(raster1) + " projected")           #outRaster2 = outRaster1 masked by Area of Interest (AOI)           outRaster2 = arcpy.sa.ExtractByMask(outRaster1, mask)  # Temporary file           outRaster2.save(os.path.join(ws_scr, "mask_{}".format(raster1)))           print ("mask_{}".format(raster1) + " masked")           #outRaster3 = Overwrited spurious data values from outRaster2           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 ("aju_{}".format(raster1) + " ajusted")           #outRaster4 = Standardized outRaster3           #ras_mean = arcpy.GetRasterProperties_management(outRaster3, "MEAN")           #ras_std = arcpy.GetRasterProperties_management(outRaster3, "STD")           ras_mean = float(arcpy.GetRasterProperties_management(outRaster3, "MEAN").getOutput(0))           ras_std = float(arcpy.GetRasterProperties_management(outRaster3, "STD").getOutput(0))           outRaster4 = ((outRaster3-ras_mean)/ras_std)  # Temporary file           outRaster4.save(os.path.join(ws_scr, "ras_std_{}".format(raster1)))           print ("ras_std_{}".format(raster1) + " standardized")           #outRaster = Final raster scaled to 0 -100           #ras_max = arcpy.GetRasterProperties_management(outRaster4, "MAXIMUM")           #ras_min = arcpy.GetRasterProperties_management(outRaster4, "MINIMUM")           ras_max = float(arcpy.GetRasterProperties_management(outRaster4, "MAXIMUM").getOutput(0))           ras_min = ras_min = float(arcpy.GetRasterProperties_management(outRaster4, "MINIMUM").getOutput(0))           outRaster = ((outRaster4-ras_min)/(ras_max-ras_min) * 100)  # Final result           outRaster.save(os.path.join(ws_out, "norm_{}".format(raster1)))           print ("norm_{}".format(raster1) + " processed")      arcpy.CheckInExtension("spatial") if __name__ == '__main__':     main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ 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 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 ("Raster list: {}".format(lst_ndvi_ras))      # get list of years      decades = GetListOfDecades(lst_ndvi_ras)      decades.sort()      #print ("Decades list: {}".format(decades))      # loop through each year      for decade in decades:           # Get rasters for given year           lst_ndvi_decade = GetListOfRasterForGivenDecade(lst_ndvi_ras, decade)           print("Decade by year list: {}".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 ("Mean raster: {}".format(out_name_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 ("Processing: {}".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_norm) - arcpy.Raster(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 ("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()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ Batch file to run on prompt: @ECHO OFF REM Runs NDVI anomaly scripts ECHO Backup original files xcopy D:\Modis250\Originais\*.tif D:\Modis250\Originais_bkp /a /c /i /y /z PAUSE ECHO Selecting decade and renaming original files to "ndvi_YYYY_dec_dd.tif" standard. C:\Python27\ArcGIS10.5\Python.exe D:\Modis250\Python\Script1.py ECHO NDVI files renamed, hit enter to pre-processing data. PAUSE ECHO Projecting, removing spurious data, standardizing and scaling ndvi decade files to 0-100. C:\Python27\ArcGIS10.5\Python.exe D:\Modis250\Python\Script2.py ECHO NDVI pre-processed, hit enter to calculate anomalies. PAUSE ECHO Calculating decade anomalies. C:\Python27\ArcGIS10.5\Python.exe D:\Modis250\Python\Script3.py Del D:\Temp /s Del D:\NDVI_10dias /s ECHO Done. PAUSE‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ... View more