|
POST
|
In line 19 is ". join(raster1)" or .format(raster1)?
... View more
06-15-2018
11:32 AM
|
0
|
28
|
1818
|
|
POST
|
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).
... View more
06-15-2018
11:26 AM
|
0
|
0
|
1818
|
|
POST
|
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()
... View more
06-15-2018
07:06 AM
|
0
|
31
|
1818
|
|
POST
|
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)
... View more
06-14-2018
01:17 PM
|
0
|
33
|
1818
|
|
POST
|
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?
... View more
06-14-2018
12:13 PM
|
0
|
36
|
3274
|
|
POST
|
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
... View more
06-14-2018
12:07 PM
|
1
|
37
|
3274
|
|
POST
|
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()
... View more
06-14-2018
11:48 AM
|
0
|
39
|
3275
|
|
POST
|
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
... View more
06-14-2018
10:41 AM
|
0
|
41
|
3274
|
|
POST
|
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.
... View more
06-14-2018
07:47 AM
|
0
|
0
|
3274
|
|
POST
|
https://community.esri.com/message/778938-python-routine-to-calculate-modis-ndvi-anomaly
... View more
06-14-2018
05:38 AM
|
1
|
0
|
808
|
|
POST
|
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()
... View more
06-14-2018
05:38 AM
|
0
|
44
|
10471
|
|
POST
|
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!
... View more
06-13-2018
12:42 PM
|
1
|
2
|
1870
|
|
POST
|
The code you provided. Exactly as above. I only inserted the "()" at lines with "print".
... View more
06-13-2018
11:37 AM
|
0
|
5
|
1870
|
|
POST
|
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()
... View more
06-13-2018
06:16 AM
|
0
|
7
|
1870
|
|
POST
|
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))
... View more
06-13-2018
05:33 AM
|
0
|
1
|
1600
|
| Title | Kudos | Posted |
|---|---|---|
| 1 | 09-18-2018 10:16 AM | |
| 1 | 06-19-2018 01:10 PM | |
| 1 | 09-03-2018 01:23 PM | |
| 1 | 06-14-2018 12:07 PM | |
| 1 | 06-14-2018 05:38 AM |
| Online Status |
Offline
|
| Date Last Visited |
11-11-2020
02:23 AM
|