Latest Contributions by BIJOYGAYEN

Hello Xander Bakker‌ sir have you check my problem? ... View more

    Thank you  Xander Bakker‌ sir, for your suggestion.I could not able to write your previously mentioned function. I think, Right now it will work smoothly. And this is my modified code. cosSenZ = arcpy.sa.Cos(senzen * math.pi / 180) cosSolZ = arcpy.sa.Cos(solzen * math.pi / 180) sinSenZ = arcpy.sa.Cos(senzen * math.pi / 180) sinSolZ = arcpy.sa.Cos(solzen * math.pi / 180) RelAzm = arcpy.sa.Abs(senazm - solazm) index = arcpy.sa.where(RelAzm > 180.0) RelAzm[index] = 360.0- RelAzm[index] index = arcpy.sa.where(RelAzm <= 180.0) RelAzm[index] = 180.0- RelAzm[index] cosRelA = arcpy.sa.Cos(RelAzm * math.pi / 180) SctAgl = arcpy.sa.ACos((-cosSolZ * cosSenZ) + ((sinSolZ * sinSenZ) * cosRelA)) OutRaster = os.path.join(outdir,'MOD022222.{0}.tif'.format("SctAgl")) SctAgl.save(OutRaster) # # RCR for TOA_BAND_1 and TOA_BAND_2 # TOA_BAND_1 B1_ROD = (0.00864 + (0.0000065)) * (0.67)**(-(3.916 + (0.074 * 0.67)+(0.05/0.67))) Pr1 = (3.0/16.0*math.pi) * (1 + (arcpy.sa.Cos(SctAgl) * arcpy.sa.Cos(SctAgl))) wr1 = 1.0 RayRef_B1 = (wr1 * B1_ROD * Pr1 ) / (4.0 * (cosSolZ * cosSenZ)) R_TOA_B1 = toa_b1 - RayRef_B1 OutRaster = os.path.join(outdir,'MOD022222.{0}.tif'.format("TOA_B1")) R_TOA_B1.save(OutRaster)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ But meanwhile, a new problem arises in where function. When I use np.where‍ function , np.where‍‍ RelAzm = arcpy.sa.Abs(senazm - solazm) index = np.where(RelAzm > 180.0) RelAzm[index] = 360.0- RelAzm[index] index = np.where(RelAzm <= 180.0) RelAzm[index] = 180.0- RelAzm[index] cosRelA = arcpy.sa.Cos(RelAzm * math.pi / 180)‍‍‍‍‍‍  it shows me this error. File "F:\DB_test_data\python_script\FINAL_DB.py", line 306, in <module> RelAzm[index] = 360.0 - RelAzm[index] TypeError: 'Raster' object has no attribute '__getitem__'‍‍‍‍‍‍  In other hands when I use arcpy.sa.where function  RelAzm = arcpy.sa.Abs(senazm - solazm) index = arcpy.sa.where(RelAzm > 180.0) RelAzm[index] = 360.0- RelAzm[index] index = arcpy.sa.where(RelAzm <= 180.0) RelAzm[index] = 180.0- RelAzm[index] cosRelA = arcpy.sa.Cos(RelAzm * math.pi / 180)‍‍‍‍‍‍ it shows this error.  Traceback (most recent call last): File "F:\DB_test_data\python_script\FINAL_DB.py", line 305, in <module> index = arcpy.sa.where(RelAzm > 180.0) AttributeError: 'module' object has no attribute 'where'‍‍‍‍ please, sir, look at this matter. ... View more

Have there any alternative way to address this issue? This conversion is very complicated so I am asking an alternative way to address this issue. ... View more

Hello Xander Bakker‌ sir, sorry for the late reply. I am confused to write the code using your mentioned function. Can you write one example code with my code? ... View more

Hello Xander Bakker‌ sir, I rechecked my code as per your suggestion, But when I try to save the output SctAgl  and toa_b1 after some processes this code gave me some error in output (the output gave me blank result with unwanted value range in two images) SctAgl: Toa_B1: Second, when I run this code it shows some unwanted messages, I cannot understand why these messages occur F:\DB_test_data\python_script\FINAL_DB.py:298: RuntimeWarning: overflow encountered in multiply cosSenZ = np.cos(senzen_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:298: RuntimeWarning: invalid value encountered in cos cosSenZ = np.cos(senzen_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:300: RuntimeWarning: overflow encountered in multiply cosSolZ = np.cos(cosSolZ_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:300: RuntimeWarning: invalid value encountered in cos cosSolZ = np.cos(cosSolZ_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:302: RuntimeWarning: overflow encountered in multiply sinSenZ = np.sin(sinSenZ_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:302: RuntimeWarning: invalid value encountered in sin sinSenZ = np.sin(sinSenZ_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:304: RuntimeWarning: overflow encountered in multiply sinSolZ = np.sin(sinSolZ_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:304: RuntimeWarning: invalid value encountered in sin sinSolZ = np.sin(sinSolZ_na * np.pi / 180)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ Third things, when I try to run this code with row data (whole India) after some processing it shows memory issue    sinSenZ_na = arcpy.RasterToNumPyArray(senzen) File "C:\Program Files (x86)\ArcGIS\Desktop10.6\ArcPy\arcpy\__init__.py", line 2275, in RasterToNumPyArray return _RasterToNumPyArray(*args, **kwargs) MemoryError‍‍‍‍ There is no option, without conversion (Raster To Numpy and then Numpy To Raster) to write the code for raster data analysis directly use cos, sin, subtract, plus, minus, division function, etc. I am totally confused to write this code. Give any solution regarding this issue. Dan Patterson‌ Joe Borgione‌Joshua Bixby‌Luke Webb‌ My_full_code: #------------------------------------------------------------------------------- # Name: module1 # Purpose: # # Author: HONEY # # Created: 18-06-2019 # Copyright: (c) HONEY 2019 # Licence: <your licence> #------------------------------------------------------------------------------- import datetime start = datetime.datetime.now() print 'start run: %s\n' % (start) import arcpy ,os ,sys,csv,errno from arcpy import env from arcpy.sa import * import datetime import re import glob import itertools import math import numpy as np arcpy.CheckOutExtension('Spatial') arcpy.env.overwriteOutput = True cellsize = "F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip\\MOD02HKM_A2017001_0530_NDVI_AA.img" #Call cloud mask images from directory d1="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" CLDMASK = glob.glob(d1 + os.sep + "*.Aerosol_Cldmask_Land_Ocean-Aerosol_Cldmask_Land_Ocean.tif") CLDMASK.sort() if len(CLDMASK) == 0: print 'Could not open the CLDMASK raster files' sys.exit(1) else: print 'The CLDMASK raster files was opened successfully' #Call AOD images from directory d2 = r"F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" AOD = glob.glob(d2 + os.sep + "*Corrected_Optical_Depth_Land_2-Corrected_Optical_Depth_Land.tif") AOD.sort() if len(AOD) == 0: print 'Could not open the AOD raster files' sys.exit(1) else: print 'The AOD raster files was opened successfully' #Call MOD02HKM_TOA-B1_TOA-B2 images from directory d3="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" TOA_B1 = glob.glob(d3 + os.sep + "*EV_500_RefSB_1-EV_500_RefSB.tif") TOA_B1.sort() if len(TOA_B1) == 0: print 'Could not open the TOA_B1 raster files' sys.exit(1) else: print 'The TOA_B1 raster files was opened successfully' TOA_B2 = glob.glob(d3 + os.sep + "*EV_500_RefSB_2-EV_500_RefSB.tif") TOA_B2.sort() if len(TOA_B2) == 0: print 'Could not open the TOA_B2 raster files' sys.exit(1) else: print 'The TOA_B2 raster files was opened successfully' #Call MOD03 (Height,SenZen,SenAzm,SolZen,SolAzm) images from directory d4="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" Height = glob.glob(d4 + os.sep + "*.Height-Height.tif") Height.sort() if len(Height) == 0: print 'Could not open the Height raster files' sys.exit(1) else: print 'The Height raster files was opened successfully' SenZen = glob.glob(d4 + os.sep + "*.SensorZenith-SensorZenith.tif") SenZen.sort() if len(SenZen) == 0: print 'Could not open the SenZen raster files' sys.exit(1) else: print 'The SenZen raster files was opened successfully' SenAzm = glob.glob(d4 + os.sep + "*.SensorAzimuth-SensorAzimuth.tif") SenAzm.sort() if len(SenAzm) == 0: print 'Could not open the SenAzm raster files' sys.exit(1) else: print 'The SenAzm raster files was opened successfully' SolZen = glob.glob(d4 + os.sep + "*.SolarZenith-SolarZenith.tif") SolZen.sort() if len(SolZen) == 0: print 'Could not open the SolZen raster files' sys.exit(1) else: print 'The SolZen raster files was opened successfully' SolAzm = glob.glob(d4 + os.sep + "*.SolarAzimuth-SolarAzimuth.tif") SolAzm.sort() if len(SolAzm) == 0: print 'Could not open the SolAzm raster files' sys.exit(1) else: print 'The SolAzm raster files was opened successfully' #Call MOD09 images from directory d5="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" Sur_B1 = glob.glob(d5 + os.sep + "*.500m_Surface_Reflectance_Band_1-500m_Surface_Reflectance_Band_1.tif") Sur_B1.sort() if len(Sur_B1) == 0: print 'Could not open the Sur_B1 raster files' sys.exit(1) else: print 'The Sur_B1 raster files was opened successfully' Sur_B2 = glob.glob(d5 + os.sep + "*.500m_Surface_Reflectance_Band_2-500m_Surface_Reflectance_Band_2.tif") Sur_B2.sort() if len(Sur_B2) == 0: print 'Could not open the Sur_B2 raster files' sys.exit(1) else: print 'The Sur_B2 raster files was opened successfully' Sur_B3 = glob.glob(d5 + os.sep + "*.500m_Surface_Reflectance_Band_3-500m_Surface_Reflectance_Band_3.tif") Sur_B3.sort() if len(Sur_B3) == 0: print 'Could not open the Sur_B3 raster files' sys.exit(1) else: print 'The Sur_B3 raster files was opened successfully' Land_cover="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip\\LAND_COVER_TYPE_1_Grid_2D.img" if len(Land_cover) == 0: print 'Could not open the Land_cover raster files' sys.exit(1) else: print 'The Land_cover raster files was opened successfully' outdir="F:\\DB_test_data\\TEST_RAY\\TEST1\\c\\d\\" cellsize = "F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip\\MOD02HKM_A2017001_0530_NDVI_AA.img" for i in range(len(CLDMASK)): AA=CLDMASK[i].split("\\")[5][0:114] print ("processing_a:"+ AA) BB=AOD[i].split("\\")[9][0:120] print ("processing_b:"+ BB) CC=TOA_B1[i].split("\\")[5][0:88] print ("processing_c:"+ CC) DD=TOA_B2[i].split("\\")[5][0:88] print ("processing_d:"+ DD) EE=Height[i].split("\\")[5][0:71] print ("processing_e:"+ EE) FF=SenZen[i].split("\\")[5][0:83] print ("processing_f:"+ FF) GG=SenAzm[i].split("\\")[5][0:85] print ("processing_g:"+ GG) HH=SolZen[i].split("\\")[5][0:81] print ("processing_h:"+ HH) II=SolAzm[i].split("\\")[5][0:83] print ("processing_i:"+ II) JJ=Sur_B1[i].split("\\")[5][0:121] print ("processing_j:"+ JJ) KK=Sur_B2[i].split("\\")[5][0:121] print ("processing_k:"+ KK) LL=Sur_B3[i].split("\\")[5][0:121] print ("processing_l:"+ LL) ####################################################################### #Cloud Mask scalling Scale_factor1 = float(1.0) add_offset1 = float(0.0) whereClause = "VALUE <= 0" # Execute SetNull outSetNull = SetNull(CLDMASK[i], CLDMASK[i], whereClause) cldmask=(outSetNull-add_offset1)*Scale_factor1 #AOD scalling Scale_factor2 = float(0.0010000000474974513) add_offset2 = float(0.0) whereClause= "VALUE = -9999" # Execute SetNull outSetNull2 = SetNull(AOD[i],AOD[i], whereClause) aod=(outSetNull2-add_offset2)*Scale_factor2 aod_1 = Con((aod >= 0.0) & (aod <= 0.1),1) #TOA_B1 Scalling Scale_factor3 = float(0.000053811826) add_offset3 = float(-0.0) whereClause= "VALUE = 65535" # Execute SetNull outSetNull3 = SetNull(TOA_B1[i],TOA_B1[i], whereClause) toa_b1=(outSetNull3-add_offset3)*Scale_factor3 #TOA_B2 Scalling Scale_factor4 = float(0.00003255546) add_offset4 = float(-0.0) whereClause= "VALUE = 65535" # Execute SetNull outSetNull4 = SetNull(TOA_B2[i],TOA_B2[i], whereClause) toa_b2=(outSetNull4-add_offset4)*Scale_factor4 #Height Scale_factor5 = float(1.0) add_offset5 = float(-0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull5 = SetNull(Height[i],Height[i], whereClause) height=(outSetNull5-add_offset5)*Scale_factor5 #Sensor Zenith angle Scale_factor6 = float(0.01) add_offset6 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull6 = SetNull(SenZen[i],SenZen[i], whereClause) senzen=(outSetNull6-add_offset6)*Scale_factor6 #Sensor azimuth angle Scale_factor7 = float(0.01) add_offset7 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull7 = SetNull(SenAzm[i],SenAzm[i], whereClause) senazm=(outSetNull7-add_offset7)*Scale_factor7 #solar Zenith angle Scale_factor8 = float(0.01) add_offset8 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull8 = SetNull(SolZen[i],SolZen[i], whereClause) solzen=(outSetNull8-add_offset8)*Scale_factor8 #solar azimuth angle Scale_factor9 = float(0.01) add_offset9 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull9 = SetNull(SolAzm[i],SolAzm[i], whereClause) solazm=(outSetNull9-add_offset9)*Scale_factor9 # OutRaster = os.path.join(outdir,'MOD04_L22.{0}.tif'.format("solazm")) # solazm.save(OutRaster) #(MODO9)surface reflectance band 1 Scale_factor_10 = float(0.0001) add_offset_10 = float(0.0) whereClause = "VALUE = -28672" # Execute SetNull outSetNull10 = SetNull(Sur_B1[i],Sur_B1[i], whereClause) sur_b1=(outSetNull10-add_offset_10)*Scale_factor_10 #(MODO9)surface reflectance band 2 Scale_factor_11 = float(0.0001) add_offset_11 = float(0.0) whereClause = "VALUE = -28672" # Execute SetNull outSetNull11 = SetNull(Sur_B2[i],Sur_B2[i], whereClause) sur_b2=(outSetNull11-add_offset_11)*Scale_factor_11 #(MODO9)surface reflectance band 3 Scale_factor_12 = float(0.0001) add_offset_12 = float(0.0) whereClause = "VALUE = -28672" # Execute SetNull outSetNull12 = SetNull(Sur_B3[i],Sur_B3[i], whereClause) sur_b3=(outSetNull12-add_offset_12)*Scale_factor_12 # Process: Extract by Cloud_Mask tempEnvironment0 = arcpy.env.cellSize arcpy.env.cellSize = "MAXOF" toa_bb1=arcpy.sa.ExtractByMask(toa_b1, cldmask) toa_bb2=arcpy.sa.ExtractByMask(toa_b2, cldmask) heightt=arcpy.sa.ExtractByMask(height, cldmask) senzenn=arcpy.sa.ExtractByMask(senzen, cldmask) senazmm=arcpy.sa.ExtractByMask(senazm, cldmask) solzenn=arcpy.sa.ExtractByMask(solzen, cldmask) solazmm=arcpy.sa.ExtractByMask(solazm, cldmask) sur_b11=arcpy.sa.ExtractByMask(sur_b1, cldmask) sur_b22=arcpy.sa.ExtractByMask(sur_b2, cldmask) sur_b33=arcpy.sa.ExtractByMask(sur_b3, cldmask) arcpy.env.cellSize = tempEnvironment0 # # Process: Extract by Mask using AOD less than 0.1 value tempEnvironment0 = arcpy.env.cellSize arcpy.env.cellSize = cellsize toa_b1=arcpy.sa.ExtractByMask(toa_bb1, aod_1) toa_b2=arcpy.sa.ExtractByMask(toa_bb2, aod_1) height=arcpy.sa.ExtractByMask(heightt, aod_1) senzen=arcpy.sa.ExtractByMask(senzenn, aod_1) senazm=arcpy.sa.ExtractByMask(senazmm, aod_1) solzen=arcpy.sa.ExtractByMask(solzenn, aod_1) solazm=arcpy.sa.ExtractByMask(solazmm, aod_1) sur_b1=arcpy.sa.ExtractByMask(sur_b11, aod_1) sur_b2=arcpy.sa.ExtractByMask(sur_b22, aod_1) sur_b3=arcpy.sa.ExtractByMask(sur_b33, aod_1) arcpy.env.cellSize = tempEnvironment0 ############################################################################# #Calculate the scattering angle from Resample MOD03 products senzen_na = arcpy.RasterToNumPyArray(senzen) cosSenZ = np.cos(senzen_na * np.pi / 180) cosSolZ_na = arcpy.RasterToNumPyArray(solzen) cosSolZ = np.cos(cosSolZ_na * np.pi / 180) sinSenZ_na = arcpy.RasterToNumPyArray(senzen) sinSenZ = np.sin(sinSenZ_na * np.pi / 180) sinSolZ_na = arcpy.RasterToNumPyArray(solzen) sinSolZ = np.sin(sinSolZ_na * np.pi / 180) senazm_na = arcpy.RasterToNumPyArray(senazm) solazm_na = arcpy.RasterToNumPyArray(solazm) RelAzm = np.abs((senazm_na ) - (solazm_na)) index = np.where(RelAzm > 180.0) RelAzm[index] = 360.0- RelAzm[index] index = np.where(RelAzm <= 180.0) RelAzm[index] = 180.0- RelAzm[index] cosRelA = np.cos(np.radians(RelAzm)) SctAgl = np.arccos((-cosSolZ * cosSenZ) + ((sinSolZ * sinSenZ) * cosRelA)) dscc=arcpy.Describe(senazm) srr=dscc.SpatialReference extt=dscc.Extent lll=arcpy.Point(extt.XMin,extt.YMin) noDataValuee=dscc.noDataValue myArraySumm = SctAgl.sum(1) myArraySumm.shape = (SctAgl.shape[0],1) myArrayPercc = (SctAgl * 1.0)/ myArraySumm newRasterr = arcpy.NumPyArrayToRaster(myArrayPercc,lll,dscc.meanCellWidth,dscc.meanCellHeight,noDataValuee) arcpy.DefineProjection_management(newRasterr, srr) OutRaster = os.path.join(outdir,'MOD022222.{0}.tif'.format("SctAgl")) newRasterr.save(OutRaster) # # RCR for TOA_BAND_1 and TOA_BAND_2 # TOA_BAND_1 dsc=arcpy.Describe(toa_b1) sr=dsc.SpatialReference ext=dsc.Extent ll=arcpy.Point(ext.XMin,ext.YMin) noDataValue=dsc.noDataValue toa_b1_na = arcpy.RasterToNumPyArray(toa_b1) B1_ROD = (0.00864 + (0.0000065)) * (0.67)**(-(3.916 + (0.074 * 0.67)+ (0.05/0.67))) Pr1 = (3.0/16.0*np.pi) * (1 + (np.cos(SctAgl) * np.cos(SctAgl))) wr1 = 1.0 RayRef_B1 = (wr1 * B1_ROD * Pr1 ) / (4.0 * (cosSolZ * cosSenZ)) myArray = toa_b1_na - RayRef_B1 myArraySum = myArray.sum(1) myArraySum.shape = (myArray.shape[0],1) myArrayPerc = (myArray * 1.0)/ myArraySum newRaster = arcpy.NumPyArrayToRaster(myArrayPerc,ll,dsc.meanCellWidth,dsc.meanCellHeight,noDataValue) arcpy.DefineProjection_management(newRaster, sr) OutRaster = os.path.join(outdir,'MOD022222.{0}.tif'.format("TOA_B1")) # toa_b2.save(OutRaster) newRaster.save(OutRaster) print arcpy.GetMessages() print 'finished run: %s\n\n' % (datetime.datetime.now() - start) ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ... View more

Hi Xander Bakker‌ sir, have you check my problem? ... View more

Hi Xander Bakker, sir Thank you for the suggestion.    I trying to modify the code as you suggested but i don't understand, why this error is occurring.  Please check this error. F:\DB_test_data\python_script\FINAL_DB.py:298: RuntimeWarning: overflow encountered in multiply cosSenZ = np.cos(senzen_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:298: RuntimeWarning: invalid value encountered in cos cosSenZ = np.cos(senzen_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:300: RuntimeWarning: overflow encountered in multiply cosSolZ = np.cos(cosSolZ_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:300: RuntimeWarning: invalid value encountered in cos cosSolZ = np.cos(cosSolZ_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:302: RuntimeWarning: overflow encountered in multiply sinSenZ = np.sin(sinSenZ_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:302: RuntimeWarning: invalid value encountered in sin sinSenZ = np.sin(sinSenZ_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:304: RuntimeWarning: overflow encountered in multiply sinSolZ = np.sin(sinSolZ_na * np.pi / 180) F:\DB_test_data\python_script\FINAL_DB.py:304: RuntimeWarning: invalid value encountered in sin sinSolZ = np.sin(sinSolZ_na * np.pi / 180) Traceback (most recent call last): File "F:\DB_test_data\python_script\FINAL_DB.py", line 320, in <module> SctAgl = np.acos((-cosSolZ * cosSenZ) + ((sinSolZ * sinSenZ) * cosRelA)) AttributeError: 'module' object has no attribute 'acos'‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ My Modified code: #------------------------------------------------------------------------------- # Name: module1 # Purpose: # # Author: HONEY # # Created: 18-06-2019 # Copyright: (c) HONEY 2019 # Licence: <your licence> #------------------------------------------------------------------------------- import datetime start = datetime.datetime.now() print 'start run: %s\n' % (start) import arcpy ,os ,sys,csv,errno from arcpy import env from arcpy.sa import * import datetime import re import glob import itertools import math import numpy as np arcpy.CheckOutExtension('Spatial') arcpy.env.overwriteOutput = True cellsize = "F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip\\MOD02HKM_A2017001_0530_NDVI_AA.img" #Call cloud mask images from directory d1="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" CLDMASK = glob.glob(d1 + os.sep + "*.Aerosol_Cldmask_Land_Ocean-Aerosol_Cldmask_Land_Ocean.tif") CLDMASK.sort() if len(CLDMASK) == 0: print 'Could not open the CLDMASK raster files' sys.exit(1) else: print 'The CLDMASK raster files was opened successfully' #Call AOD images from directory d2 = r"F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" AOD = glob.glob(d2 + os.sep + "*Corrected_Optical_Depth_Land_2-Corrected_Optical_Depth_Land.tif") AOD.sort() if len(AOD) == 0: print 'Could not open the AOD raster files' sys.exit(1) else: print 'The AOD raster files was opened successfully' #Call MOD02HKM_TOA-B1_TOA-B2 images from directory d3="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" TOA_B1 = glob.glob(d3 + os.sep + "*EV_500_RefSB_1-EV_500_RefSB.tif") TOA_B1.sort() if len(TOA_B1) == 0: print 'Could not open the TOA_B1 raster files' sys.exit(1) else: print 'The TOA_B1 raster files was opened successfully' TOA_B2 = glob.glob(d3 + os.sep + "*EV_500_RefSB_2-EV_500_RefSB.tif") TOA_B2.sort() if len(TOA_B2) == 0: print 'Could not open the TOA_B2 raster files' sys.exit(1) else: print 'The TOA_B2 raster files was opened successfully' #Call MOD03 (Height,SenZen,SenAzm,SolZen,SolAzm) images from directory d4="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" Height = glob.glob(d4 + os.sep + "*.Height-Height.tif") Height.sort() if len(Height) == 0: print 'Could not open the Height raster files' sys.exit(1) else: print 'The Height raster files was opened successfully' SenZen = glob.glob(d4 + os.sep + "*.SensorZenith-SensorZenith.tif") SenZen.sort() if len(SenZen) == 0: print 'Could not open the SenZen raster files' sys.exit(1) else: print 'The SenZen raster files was opened successfully' SenAzm = glob.glob(d4 + os.sep + "*.SensorAzimuth-SensorAzimuth.tif") SenAzm.sort() if len(SenAzm) == 0: print 'Could not open the SenAzm raster files' sys.exit(1) else: print 'The SenAzm raster files was opened successfully' SolZen = glob.glob(d4 + os.sep + "*.SolarZenith-SolarZenith.tif") SolZen.sort() if len(SolZen) == 0: print 'Could not open the SolZen raster files' sys.exit(1) else: print 'The SolZen raster files was opened successfully' SolAzm = glob.glob(d4 + os.sep + "*.SolarAzimuth-SolarAzimuth.tif") SolAzm.sort() if len(SolAzm) == 0: print 'Could not open the SolAzm raster files' sys.exit(1) else: print 'The SolAzm raster files was opened successfully' #Call MOD09 images from directory d5="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" Sur_B1 = glob.glob(d5 + os.sep + "*.500m_Surface_Reflectance_Band_1-500m_Surface_Reflectance_Band_1.tif") Sur_B1.sort() if len(Sur_B1) == 0: print 'Could not open the Sur_B1 raster files' sys.exit(1) else: print 'The Sur_B1 raster files was opened successfully' Sur_B2 = glob.glob(d5 + os.sep + "*.500m_Surface_Reflectance_Band_2-500m_Surface_Reflectance_Band_2.tif") Sur_B2.sort() if len(Sur_B2) == 0: print 'Could not open the Sur_B2 raster files' sys.exit(1) else: print 'The Sur_B2 raster files was opened successfully' Sur_B3 = glob.glob(d5 + os.sep + "*.500m_Surface_Reflectance_Band_3-500m_Surface_Reflectance_Band_3.tif") Sur_B3.sort() if len(Sur_B3) == 0: print 'Could not open the Sur_B3 raster files' sys.exit(1) else: print 'The Sur_B3 raster files was opened successfully' Land_cover="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip\\LAND_COVER_TYPE_1_Grid_2D.img" if len(Land_cover) == 0: print 'Could not open the Land_cover raster files' sys.exit(1) else: print 'The Land_cover raster files was opened successfully' outdir="F:\\DB_test_data\\TEST_RAY\\TEST1\\c\\d\\" cellsize = "F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip\\MOD02HKM_A2017001_0530_NDVI_AA.img" for i in range(len(CLDMASK)): AA=CLDMASK[i].split("\\")[5][0:114] print ("processing_a:"+ AA) BB=AOD[i].split("\\")[9][0:120] print ("processing_b:"+ BB) CC=TOA_B1[i].split("\\")[5][0:88] print ("processing_c:"+ CC) DD=TOA_B2[i].split("\\")[5][0:88] print ("processing_d:"+ DD) EE=Height[i].split("\\")[5][0:71] print ("processing_e:"+ EE) FF=SenZen[i].split("\\")[5][0:83] print ("processing_f:"+ FF) GG=SenAzm[i].split("\\")[5][0:85] print ("processing_g:"+ GG) HH=SolZen[i].split("\\")[5][0:81] print ("processing_h:"+ HH) II=SolAzm[i].split("\\")[5][0:83] print ("processing_i:"+ II) JJ=Sur_B1[i].split("\\")[5][0:121] print ("processing_j:"+ JJ) KK=Sur_B2[i].split("\\")[5][0:121] print ("processing_k:"+ KK) LL=Sur_B3[i].split("\\")[5][0:121] print ("processing_l:"+ LL) ####################################################################### #Cloud Mask scalling Scale_factor1 = float(1.0) add_offset1 = float(0.0) whereClause = "VALUE <= 0" # Execute SetNull outSetNull = SetNull(CLDMASK[i], CLDMASK[i], whereClause) cldmask=(outSetNull-add_offset1)*Scale_factor1 #AOD scalling Scale_factor2 = float(0.0010000000474974513) add_offset2 = float(0.0) whereClause= "VALUE = -9999" # Execute SetNull outSetNull2 = SetNull(AOD[i],AOD[i], whereClause) aod=(outSetNull2-add_offset2)*Scale_factor2 aod_1 = Con((aod >= 0.0) & (aod <= 0.1),1) #TOA_B1 Scalling Scale_factor3 = float(0.000053811826) add_offset3 = float(-0.0) whereClause= "VALUE = 65535" # Execute SetNull outSetNull3 = SetNull(TOA_B1[i],TOA_B1[i], whereClause) toa_b1=(outSetNull3-add_offset3)*Scale_factor3 #TOA_B2 Scalling Scale_factor4 = float(0.00003255546) add_offset4 = float(-0.0) whereClause= "VALUE = 65535" # Execute SetNull outSetNull4 = SetNull(TOA_B2[i],TOA_B2[i], whereClause) toa_b2=(outSetNull4-add_offset4)*Scale_factor4 #Height Scale_factor5 = float(1.0) add_offset5 = float(-0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull5 = SetNull(Height[i],Height[i], whereClause) height=(outSetNull5-add_offset5)*Scale_factor5 #Sensor Zenith angle Scale_factor6 = float(0.01) add_offset6 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull6 = SetNull(SenZen[i],SenZen[i], whereClause) senzen=(outSetNull6-add_offset6)*Scale_factor6 #Sensor azimuth angle Scale_factor7 = float(0.01) add_offset7 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull7 = SetNull(SenAzm[i],SenAzm[i], whereClause) senazm=(outSetNull7-add_offset7)*Scale_factor7 #solar Zenith angle Scale_factor8 = float(0.01) add_offset8 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull8 = SetNull(SolZen[i],SolZen[i], whereClause) solzen=(outSetNull8-add_offset8)*Scale_factor8 #solar azimuth angle Scale_factor9 = float(0.01) add_offset9 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull9 = SetNull(SolAzm[i],SolAzm[i], whereClause) solazm=(outSetNull9-add_offset9)*Scale_factor9 # OutRaster = os.path.join(outdir,'MOD04_L22.{0}.tif'.format("solazm")) # solazm.save(OutRaster) #(MODO9)surface reflectance band 1 Scale_factor_10 = float(0.0001) add_offset_10 = float(0.0) whereClause = "VALUE = -28672" # Execute SetNull outSetNull10 = SetNull(Sur_B1[i],Sur_B1[i], whereClause) sur_b1=(outSetNull10-add_offset_10)*Scale_factor_10 #(MODO9)surface reflectance band 2 Scale_factor_11 = float(0.0001) add_offset_11 = float(0.0) whereClause = "VALUE = -28672" # Execute SetNull outSetNull11 = SetNull(Sur_B2[i],Sur_B2[i], whereClause) sur_b2=(outSetNull11-add_offset_11)*Scale_factor_11 #(MODO9)surface reflectance band 3 Scale_factor_12 = float(0.0001) add_offset_12 = float(0.0) whereClause = "VALUE = -28672" # Execute SetNull outSetNull12 = SetNull(Sur_B3[i],Sur_B3[i], whereClause) sur_b3=(outSetNull12-add_offset_12)*Scale_factor_12 # Process: Extract by Cloud_Mask tempEnvironment0 = arcpy.env.cellSize arcpy.env.cellSize = "MAXOF" toa_bb1=arcpy.sa.ExtractByMask(toa_b1, cldmask) toa_bb2=arcpy.sa.ExtractByMask(toa_b2, cldmask) heightt=arcpy.sa.ExtractByMask(height, cldmask) senzenn=arcpy.sa.ExtractByMask(senzen, cldmask) senazmm=arcpy.sa.ExtractByMask(senazm, cldmask) solzenn=arcpy.sa.ExtractByMask(solzen, cldmask) solazmm=arcpy.sa.ExtractByMask(solazm, cldmask) sur_b11=arcpy.sa.ExtractByMask(sur_b1, cldmask) sur_b22=arcpy.sa.ExtractByMask(sur_b2, cldmask) sur_b33=arcpy.sa.ExtractByMask(sur_b3, cldmask) arcpy.env.cellSize = tempEnvironment0 # # Process: Extract by Mask using AOD less than 0.1 value tempEnvironment0 = arcpy.env.cellSize arcpy.env.cellSize = cellsize toa_b1=arcpy.sa.ExtractByMask(toa_bb1, aod_1) toa_b2=arcpy.sa.ExtractByMask(toa_bb2, aod_1) height=arcpy.sa.ExtractByMask(heightt, aod_1) senzen=arcpy.sa.ExtractByMask(senzenn, aod_1) senazm=arcpy.sa.ExtractByMask(senazmm, aod_1) solzen=arcpy.sa.ExtractByMask(solzenn, aod_1) solazm=arcpy.sa.ExtractByMask(solazmm, aod_1) sur_b1=arcpy.sa.ExtractByMask(sur_b11, aod_1) sur_b2=arcpy.sa.ExtractByMask(sur_b22, aod_1) sur_b3=arcpy.sa.ExtractByMask(sur_b33, aod_1) arcpy.env.cellSize = tempEnvironment0 ############################################################################# #Calculate the scattering angle from Resample MOD03 products senzen_na = arcpy.RasterToNumPyArray(senzen) cosSenZ = np.cos(senzen_na * np.pi / 180) cosSolZ_na = arcpy.RasterToNumPyArray(solzen) cosSolZ = np.cos(cosSolZ_na * np.pi / 180) sinSenZ_na = arcpy.RasterToNumPyArray(senzen) sinSenZ = np.sin(sinSenZ_na * np.pi / 180) sinSolZ_na = arcpy.RasterToNumPyArray(solzen) sinSolZ = np.sin(sinSolZ_na * np.pi / 180) senazm_na = arcpy.RasterToNumPyArray(senazm) solazm_na = arcpy.RasterToNumPyArray(solazm) RelAzm = np.abs((senazm_na ) - (solazm_na)) index = np.where(RelAzm > 180.0) RelAzm[index] = 360.0- RelAzm[index] index = np.where(RelAzm <= 180.0) RelAzm[index] = 180.0- RelAzm[index] cosRe = np.cos(np.radians(RelAzm)) SctAgl = np.acos((-cosSolZ * cosSenZ) + ((sinSolZ * sinSenZ) * cosRelA)) # # RCR for TOA_BAND_1 and TOA_BAND_2 # TOA_BAND_1 B1_ROD = (0.00864 + (0.0000065)) * (0.67)**(-(3.916 + (0.074 * 0.67)+ (0.05/0.67))) Pr1 = (3.0/16.0*np.pi) * (1 + (np.cos(SctAgl) * np.cos(SctAgl))) wr1 = 1.0 RayRef_B1 = (wr1 * B1_ROD * Pr1 ) / (4.0 * (cosSolZ * cosSenZ)) toa_b2_na = arcpy.RasterToNumPyArray(toa_b2) RCR_B1 = toa_b1 - RayRef_B1 #TOA_BAND_2 B2_ROD = (0.00864 + (0.0000065)) * (0.86)**(-(3.916 + (0.074 * 0.86)+ (0.05/0.86))) Pr2 = (3.0/16.0*np.pi) * (1 + (np.cos(SctAgl) * np.cos(SctAgl))) wr2 = 1.0 RayRef_B2 = (wr2 * B2_ROD * Pr2 ) / (4.0 * (cosSolZ * cosSenZ)) toa_b2_na = arcpy.RasterToNumPyArray(toa_b2) RCR_B2 = toa_b2 - RayRef_B2 print arcpy.GetMessages() print 'finished run: %s\n\n' % (datetime.datetime.now() - start) ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ... View more

OK Xander Bakker‌ sir, You can take your time. ... View more

Xander Bakker‌ Sir, Did you check my code and attached TEST DATA? ... View more

Xander Bakker My code is: You can test my code using my "test_data" that I attached . #------------------------------------------------------------------------------- # Name: module1 # Purpose: # # Author: HONEY # # Created: 18-06-2019 # Copyright: (c) HONEY 2019 # Licence: <your licence> #------------------------------------------------------------------------------- import datetime start = datetime.datetime.now() print 'start run: %s\n' % (start) import arcpy ,os ,sys,csv,errno from arcpy import env from arcpy.sa import * import datetime import re import numpy as np import glob import itertools arcpy.CheckOutExtension('Spatial') arcpy.env.overwriteOutput = True cellsize = "F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip\\MOD02HKM_A2017001_0530_NDVI_AA.img" #Call cloud mask images from directory d1="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" CLDMASK = glob.glob(d1 + os.sep + "*.Aerosol_Cldmask_Land_Ocean-Aerosol_Cldmask_Land_Ocean.tif") CLDMASK.sort() if len(CLDMASK) == 0: print 'Could not open the CLDMASK raster files' sys.exit(1) else: print 'The CLDMASK raster files was opened successfully' #Call AOD images from directory d2 = r"F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" AOD = glob.glob(d2 + os.sep + "*Corrected_Optical_Depth_Land_2-Corrected_Optical_Depth_Land.tif") AOD.sort() if len(AOD) == 0: print 'Could not open the AOD raster files' sys.exit(1) else: print 'The AOD raster files was opened successfully' #Call MOD02HKM_TOA-B1_TOA-B2 images from directory d3="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" TOA_B1 = glob.glob(d3 + os.sep + "*EV_500_RefSB_1-EV_500_RefSB.tif") TOA_B1.sort() if len(TOA_B1) == 0: print 'Could not open the TOA_B1 raster files' sys.exit(1) else: print 'The TOA_B1 raster files was opened successfully' TOA_B2 = glob.glob(d3 + os.sep + "*EV_500_RefSB_2-EV_500_RefSB.tif") TOA_B2.sort() if len(TOA_B2) == 0: print 'Could not open the TOA_B2 raster files' sys.exit(1) else: print 'The TOA_B2 raster files was opened successfully' #Call MOD03 (Height,SenZen,SenAzm,SolZen,SolAzm) images from directory d4="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" Height = glob.glob(d4 + os.sep + "*.Height-Height.tif") Height.sort() if len(Height) == 0: print 'Could not open the Height raster files' sys.exit(1) else: print 'The Height raster files was opened successfully' SenZen = glob.glob(d4 + os.sep + "*.SensorZenith-SensorZenith.tif") SenZen.sort() if len(SenZen) == 0: print 'Could not open the SenZen raster files' sys.exit(1) else: print 'The SenZen raster files was opened successfully' SenAzm = glob.glob(d4 + os.sep + "*.SensorAzimuth-SensorAzimuth.tif") SenAzm.sort() if len(SenAzm) == 0: print 'Could not open the SenAzm raster files' sys.exit(1) else: print 'The SenAzm raster files was opened successfully' SolZen = glob.glob(d4 + os.sep + "*.SolarZenith-SolarZenith.tif") SolZen.sort() if len(SolZen) == 0: print 'Could not open the SolZen raster files' sys.exit(1) else: print 'The SolZen raster files was opened successfully' SolAzm = glob.glob(d4 + os.sep + "*.SolarAzimuth-SolarAzimuth.tif") SolAzm.sort() if len(SolAzm) == 0: print 'Could not open the SolAzm raster files' sys.exit(1) else: print 'The SolAzm raster files was opened successfully' #Call MOD09 images from directory d5="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" Sur_B1 = glob.glob(d5 + os.sep + "*.500m_Surface_Reflectance_Band_1-500m_Surface_Reflectance_Band_1.tif") Sur_B1.sort() if len(Sur_B1) == 0: print 'Could not open the Sur_B1 raster files' sys.exit(1) else: print 'The Sur_B1 raster files was opened successfully' Sur_B2 = glob.glob(d5 + os.sep + "*.500m_Surface_Reflectance_Band_2-500m_Surface_Reflectance_Band_2.tif") Sur_B2.sort() if len(Sur_B2) == 0: print 'Could not open the Sur_B2 raster files' sys.exit(1) else: print 'The Sur_B2 raster files was opened successfully' Sur_B3 = glob.glob(d5 + os.sep + "*.500m_Surface_Reflectance_Band_3-500m_Surface_Reflectance_Band_3.tif") Sur_B3.sort() if len(Sur_B3) == 0: print 'Could not open the Sur_B3 raster files' sys.exit(1) else: print 'The Sur_B3 raster files was opened successfully' Land_cover="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip\\LAND_COVER_TYPE_1_Grid_2D.img" if len(Land_cover) == 0: print 'Could not open the Land_cover raster files' sys.exit(1) else: print 'The Land_cover raster files was opened successfully' outdir="F:\\DB_test_data\\TEST_RAY\\TEST1\\c\\d\\" cellsize = "F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip\\MOD02HKM_A2017001_0530_NDVI_AA.img" for a,b,c,d,e,f,g,h,i,j,k,l, in zip (CLDMASK,AOD,TOA_B1,TOA_B2,Height,SenZen,SenAzm,SolZen,SolAzm,Sur_B1,Sur_B2,Sur_B3): #arcpy.AddMessage("processing:{}".format(b.split('\\')[4][9:30])) AA=a.split("\\")[5][0:114] print ("processing_a:"+ AA) BB=b.split("\\")[9][0:120] print ("processing_b:"+ BB) CC=c.split("\\")[5][0:88] print ("processing_c:"+ CC) DD=d.split("\\")[5][0:88] print ("processing_d:"+ DD) EE=e.split("\\")[5][0:71] print ("processing_e:"+ EE) FF=f.split("\\")[5][0:83] print ("processing_f:"+ FF) GG=g.split("\\")[5][0:85] print ("processing_g:"+ GG) HH=h.split("\\")[5][0:81] print ("processing_h:"+ HH) II=i.split("\\")[5][0:83] print ("processing_i:"+ II) JJ=j.split("\\")[5][0:121] print ("processing_j:"+ JJ) KK=k.split("\\")[5][0:121] print ("processing_k:"+ KK) LL=l.split("\\")[5][0:121] print ("processing_l:"+ LL) ####################################################################### #Cloud Mask scalling Scale_factor1 = float(1.0) add_offset1 = float(0.0) whereClause = "VALUE <= 0" # Execute SetNull outSetNull = SetNull(a, a, whereClause) cldmask=(outSetNull-add_offset1)*Scale_factor1 #AOD scalling Scale_factor2 = float(0.0010000000474974513) add_offset2 = float(0.0) whereClause= "VALUE = -9999" # Execute SetNull outSetNull2 = SetNull(b,b, whereClause) aod=(outSetNull2-add_offset2)*Scale_factor2 aod_1 = Con((aod >= 0.0) & (aod <= 0.1),1) #TOA_B1 Scalling Scale_factor3 = float(0.000053811826) add_offset3 = float(-0.0) whereClause= "VALUE = 65535" # Execute SetNull outSetNull3 = SetNull(c,c, whereClause) toa_b1=(outSetNull3-add_offset3)*Scale_factor3 #TOA_B2 Scalling Scale_factor4 = float(0.00003255546) add_offset4 = float(-0.0) whereClause= "VALUE = 65535" # Execute SetNull outSetNull4 = SetNull(d,d, whereClause) toa_b2=(outSetNull4-add_offset4)*Scale_factor4 #Height Scale_factor5 = float(1.0) add_offset5 = float(-0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull5 = SetNull(e,e, whereClause) height=(outSetNull5-add_offset5)*Scale_factor5 #Sensor Zenith angle Scale_factor6 = float(0.01) add_offset6 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull6 = SetNull(f,f, whereClause) senzen=(outSetNull6-add_offset6)*Scale_factor6 #Sensor azimuth angle Scale_factor7 = float(0.01) add_offset7 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull7 = SetNull(g,g, whereClause) senazm=(outSetNull7-add_offset7)*Scale_factor7 #solar Zenith angle Scale_factor8 = float(0.01) add_offset8 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull8 = SetNull(h,h, whereClause) solzen=(outSetNull8-add_offset8)*Scale_factor8 #solar azimuth angle Scale_factor9 = float(0.01) add_offset9 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull9 = SetNull(i,i, whereClause) solazm=(outSetNull9-add_offset9)*Scale_factor9 # OutRaster = os.path.join(outdir,'MOD04_L22.{0}.tif'.format("solazm")) # solazm.save(OutRaster) #(MODO9)surface reflectance band 1 Scale_factor_10 = float(0.0001) add_offset_10 = float(0.0) whereClause = "VALUE = -28672" # Execute SetNull outSetNull10 = SetNull(j,j, whereClause) sur_b1=(outSetNull10-add_offset_10)*Scale_factor_10 #(MODO9)surface reflectance band 2 Scale_factor_11 = float(0.0001) add_offset_11 = float(0.0) whereClause = "VALUE = -28672" # Execute SetNull outSetNull11 = SetNull(k,k, whereClause) sur_b2=(outSetNull11-add_offset_11)*Scale_factor_11 #(MODO9)surface reflectance band 3 Scale_factor_12 = float(0.0001) add_offset_12 = float(0.0) whereClause = "VALUE = -28672" # Execute SetNull outSetNull12 = SetNull(l,l, whereClause) sur_b3=(outSetNull12-add_offset_12)*Scale_factor_12 # Process: Extract by Cloud_Mask tempEnvironment0 = arcpy.env.cellSize arcpy.env.cellSize = "MAXOF" toa_bb1=arcpy.sa.ExtractByMask(toa_b1, cldmask) toa_bb2=arcpy.sa.ExtractByMask(toa_b2, cldmask) heightt=arcpy.sa.ExtractByMask(height, cldmask) senzenn=arcpy.sa.ExtractByMask(senzen, cldmask) senazmm=arcpy.sa.ExtractByMask(senazm, cldmask) solzenn=arcpy.sa.ExtractByMask(solzen, cldmask) solazmm=arcpy.sa.ExtractByMask(solazm, cldmask) sur_b11=arcpy.sa.ExtractByMask(sur_b1, cldmask) sur_b22=arcpy.sa.ExtractByMask(sur_b2, cldmask) sur_b33=arcpy.sa.ExtractByMask(sur_b3, cldmask) arcpy.env.cellSize = tempEnvironment0 # # Process: Extract by Mask using AOD less than 0.1 value tempEnvironment0 = arcpy.env.cellSize arcpy.env.cellSize = cellsize toa_b1=arcpy.sa.ExtractByMask(toa_bb1, aod_1) toa_b2=arcpy.sa.ExtractByMask(toa_bb2, aod_1) height=arcpy.sa.ExtractByMask(heightt, aod_1) senzen=arcpy.sa.ExtractByMask(senzenn, aod_1) senazm=arcpy.sa.ExtractByMask(senazmm, aod_1) solzen=arcpy.sa.ExtractByMask(solzenn, aod_1) solazm=arcpy.sa.ExtractByMask(solazmm, aod_1) sur_b1=arcpy.sa.ExtractByMask(sur_b11, aod_1) sur_b2=arcpy.sa.ExtractByMask(sur_b22, aod_1) sur_b3=arcpy.sa.ExtractByMask(sur_b33, aod_1) arcpy.env.cellSize = tempEnvironment0 ############################################################################# #Calculate the scattering angle from Resample MOD03 products cosSenZ = np.cos(senzen * np.pi / 180) cosSolZ = np.cos(solzen * np.pi / 180) sinSenZ = np.sin(senzen * np.pi / 180) sinSolZ = np.sin(solzen * np.pi / 180) RelAzm = np.abs((senazm ) - (solazm)) index = np.where(RelAzm > 180.0) RelAzm[index] = 360.0- RelAzm[index] index = np.where(RelAzm <= 180.0) RelAzm[index] = 180.0- RelAzm[index] cosRe = np.cos(np.radians(RelAzm)) SctAgl = np.acos((-cosSolZ * cosSenZ) + ((sinSolZ * sinSenZ) * cosRelA)) # RCR for TOA_BAND_1 and TOA_BAND_2 # TOA_BAND_1 B1_ROD = (0.00864 + (0.0000065)) * (0.67)**(-(3.916 + (0.074 * 0.67)+ (0.05/0.67))) Pr1 = (3.0/16.0*np.pi) * (1 + (np.cos(SctAgl) * np.cos(SctAgl))) wr1 = 1.0 RayRef_B1 = (wr1 * B1_ROD * Pr1 ) / (4.0 * (cosSolZ * cosSenZ)) RCR_B1 = toa_b1 - RayRef_B1 OutRaster = os.path.join(outdir,'MOD02HKM_L22.{0}.tif'.format("RAY_CORRECT_TOA_B1")) RCR_B1.save(OutRaster) #TOA_BAND_2 B2_ROD = (0.00864 + (0.0000065)) * (0.86)**(-(3.916 + (0.074 * 0.86)+ (0.05/0.86))) Pr2 = (3.0/16.0*np.pi) * (1 + (np.cos(SctAgl) * np.cos(SctAgl))) wr2 = 1.0 RayRef_B2 = (wr2 * B2_ROD * Pr2 ) / (4.0 * (cosSolZ * cosSenZ)) RCR_B2 = toa_b2 - RayRef_B2 OutRaster = os.path.join(outdir,'MOD02HKM_L22.{0}.tif'.format("RAY_CORRECT_TOA_B2")) RCR_B2.save(OutRaster) # Calculating NDVI float raster ndvi_raster = Divide(Float(Raster(RCR_B2) - Raster(RCR_B1)), Float(Raster(RCR_B2) + Raster(RCR_B1))) OutRaster = os.path.join(outdir,'MOD0_L22.{0}.tif'.format("NDVI")) ndvi_raster.save(OutRaster) print arcpy.GetMessages() print 'finished run: %s\n\n' % (datetime.datetime.now() - start) ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ... View more

Xander Bakker‌ , sir I wrote all the "np.cos" ,"np.pi" and "np.sin" in lower case but it is now giving another error. Traceback (most recent call last): File "F:\DB_test_data\python_script\FINAL_DB.py", line 295, in <module> cosSenZ = np.cos(senzen * np.pi / 180) AttributeError: 'Raster' object has no attribute 'cos'‍‍‍‍ code: cosSenZ = np.cos(senzen * np.pi / 180) cosSolZ = np.cos(solzen * np.pi / 180) sinSenZ = np.sin(senzen * np.pi / 180) sinSolZ = np.sin(solzen * np.pi / 180) RelAzm = np.abs((senazm ) - (solazm)) index = np.where(RelAzm > 180.0) RelAzm[index] = 360.0- RelAzm[index] index = np.where(RelAzm <= 180.0) RelAzm[index] = 180.0- RelAzm[index] cosRe = np.cos(np.radians(RelAzm)) SctAgl = np.acos((-cosSolZ * cosSenZ) + ((sinSolZ * sinSenZ) * cosRelA)) # RCR for TOA_BAND_1 and TOA_BAND_2 # TOA_BAND_1 B1_ROD = (0.00864 + (0.0000065)) * (0.67)**(-(3.916 + (0.074 * 0.67)+ (0.05/0.67))) Pr1 = (3.0/16.0*np.pi) * (1 + (np.cos(SctAgl) * np.cos(SctAgl))) wr1 = 1.0 RayRef_B1 = (wr1 * B1_ROD * Pr1 ) / (4.0 * (cosSolZ * cosSenZ)) RCR_B1 = toa_b1 - RayRef_B1 OutRaster = os.path.join(outdir,'MOD02HKM_L22.{0}.tif'.format("RAY_CORRECT_TOA_B1")) RCR_B1.save(OutRaster) #TOA_BAND_2 B2_ROD = (0.00864 + (0.0000065)) * (0.86)**(-(3.916 + (0.074 * 0.86)+ (0.05/0.86))) Pr2 = (3.0/16.0*np.pi) * (1 + (np.cos(SctAgl) * np.cos(SctAgl))) wr2 = 1.0 RayRef_B2 = (wr2 * B2_ROD * Pr2 ) / (4.0 * (cosSolZ * cosSenZ)) RCR_B2 = toa_b2 - RayRef_B2 OutRaster = os.path.join(outdir,'MOD02HKM_L22.{0}.tif'.format("RAY_CORRECT_TOA_B2")) RCR_B2.save(OutRaster) # Calculating NDVI float raster ndvi_raster = Divide(Float(Raster(RCR_B2) - Raster(RCR_B1)), Float(Raster(RCR_B2) + Raster(RCR_B1))) OutRaster = os.path.join(outdir,'MOD0_L22.{0}.tif'.format("NDVI")) ndvi_raster.save(OutRaster)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ... View more

Thank you Xander Bakker, sir, for your suggestion. Yes is the need to avoid "In-memory" function when we analysis Big size data.I wrote the code a bit differently to remove the "in-memory" function, and it worked for me. Scale_factor1 = float(1.0) add_offset1 = float(0.0) whereClause = "VALUE <= 0" # Execute SetNull outSetNull = SetNull(a, a, whereClause) cldmask=(outSetNull-add_offset1)*Scale_factor1 OutRaster = os.path.join(outdir,'MOD04_L22.{0}.tif'.format("CLDMASK")) cldmask.save(OutRaster)‍‍‍‍‍‍‍‍ I have another problem regarding the trigonometric calculation in raster data. I am trying to write simple trigonometric calculation(code line 295-330) using numpy but it gave me one error. Traceback (most recent call last): File "F:\DB_test_data\python_script\FINAL_DB.py", line 295, in <module> CosSenZ = np.Cos(senzen * np.pi / 180) AttributeError: 'module' object has no attribute 'Cos'‍‍‍‍ MY_CODE( trigonometric calculation(code line 295-330)) #------------------------------------------------------------------------------- # Name: module1 # Purpose: # # Author: HONEY # # Created: 18-06-2019 # Copyright: (c) HONEY 2019 # Licence: <your licence> #------------------------------------------------------------------------------- import datetime start = datetime.datetime.now() print 'start run: %s\n' % (start) import arcpy ,os ,sys,csv,errno from arcpy import env from arcpy.sa import * import datetime import re import numpy as np import glob import itertools arcpy.CheckOutExtension('Spatial') arcpy.env.overwriteOutput = True cellsize = "F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip\\MOD02HKM_A2017001_0530_NDVI_AA.img" #Call cloud mask images from directory d1="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" CLDMASK = glob.glob(d1 + os.sep + "*.Aerosol_Cldmask_Land_Ocean-Aerosol_Cldmask_Land_Ocean.tif") CLDMASK.sort() if len(CLDMASK) == 0: print 'Could not open the CLDMASK raster files' sys.exit(1) else: print 'The CLDMASK raster files was opened successfully' #Call AOD images from directory d2 = r"F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" AOD = glob.glob(d2 + os.sep + "*Corrected_Optical_Depth_Land_2-Corrected_Optical_Depth_Land.tif") AOD.sort() if len(AOD) == 0: print 'Could not open the AOD raster files' sys.exit(1) else: print 'The AOD raster files was opened successfully' #Call MOD02HKM_TOA-B1_TOA-B2 images from directory d3="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" TOA_B1 = glob.glob(d3 + os.sep + "*EV_500_RefSB_1-EV_500_RefSB.tif") TOA_B1.sort() if len(TOA_B1) == 0: print 'Could not open the TOA_B1 raster files' sys.exit(1) else: print 'The TOA_B1 raster files was opened successfully' TOA_B2 = glob.glob(d3 + os.sep + "*EV_500_RefSB_2-EV_500_RefSB.tif") TOA_B2.sort() if len(TOA_B2) == 0: print 'Could not open the TOA_B2 raster files' sys.exit(1) else: print 'The TOA_B2 raster files was opened successfully' #Call MOD03 (Height,SenZen,SenAzm,SolZen,SolAzm) images from directory d4="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" Height = glob.glob(d4 + os.sep + "*.Height-Height.tif") Height.sort() if len(Height) == 0: print 'Could not open the Height raster files' sys.exit(1) else: print 'The Height raster files was opened successfully' SenZen = glob.glob(d4 + os.sep + "*.SensorZenith-SensorZenith.tif") SenZen.sort() if len(SenZen) == 0: print 'Could not open the SenZen raster files' sys.exit(1) else: print 'The SenZen raster files was opened successfully' SenAzm = glob.glob(d4 + os.sep + "*.SensorAzimuth-SensorAzimuth.tif") SenAzm.sort() if len(SenAzm) == 0: print 'Could not open the SenAzm raster files' sys.exit(1) else: print 'The SenAzm raster files was opened successfully' SolZen = glob.glob(d4 + os.sep + "*.SolarZenith-SolarZenith.tif") SolZen.sort() if len(SolZen) == 0: print 'Could not open the SolZen raster files' sys.exit(1) else: print 'The SolZen raster files was opened successfully' SolAzm = glob.glob(d4 + os.sep + "*.SolarAzimuth-SolarAzimuth.tif") SolAzm.sort() if len(SolAzm) == 0: print 'Could not open the SolAzm raster files' sys.exit(1) else: print 'The SolAzm raster files was opened successfully' #Call MOD09 images from directory d5="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip" Sur_B1 = glob.glob(d5 + os.sep + "*.500m_Surface_Reflectance_Band_1-500m_Surface_Reflectance_Band_1.tif") Sur_B1.sort() if len(Sur_B1) == 0: print 'Could not open the Sur_B1 raster files' sys.exit(1) else: print 'The Sur_B1 raster files was opened successfully' Sur_B2 = glob.glob(d5 + os.sep + "*.500m_Surface_Reflectance_Band_2-500m_Surface_Reflectance_Band_2.tif") Sur_B2.sort() if len(Sur_B2) == 0: print 'Could not open the Sur_B2 raster files' sys.exit(1) else: print 'The Sur_B2 raster files was opened successfully' Sur_B3 = glob.glob(d5 + os.sep + "*.500m_Surface_Reflectance_Band_3-500m_Surface_Reflectance_Band_3.tif") Sur_B3.sort() if len(Sur_B3) == 0: print 'Could not open the Sur_B3 raster files' sys.exit(1) else: print 'The Sur_B3 raster files was opened successfully' Land_cover="F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip\\LAND_COVER_TYPE_1_Grid_2D.img" if len(Land_cover) == 0: print 'Could not open the Land_cover raster files' sys.exit(1) else: print 'The Land_cover raster files was opened successfully' outdir="F:\\DB_test_data\\TEST_RAY\\TEST1\\c\\d\\" cellsize = "F:\\DB_test_data\\VAR2\\TEST_DATA\\Clip\\MOD02HKM_A2017001_0530_NDVI_AA.img" for a,b,c,d,e,f,g,h,i,j,k,l, in zip (CLDMASK,AOD,TOA_B1,TOA_B2,Height,SenZen,SenAzm,SolZen,SolAzm,Sur_B1,Sur_B2,Sur_B3): #arcpy.AddMessage("processing:{}".format(b.split('\\')[4][9:30])) AA=a.split("\\")[5][0:114] print ("processing_a:"+ AA) BB=b.split("\\")[9][0:120] print ("processing_b:"+ BB) CC=c.split("\\")[5][0:88] print ("processing_c:"+ CC) DD=d.split("\\")[5][0:88] print ("processing_d:"+ DD) EE=e.split("\\")[5][0:71] print ("processing_e:"+ EE) FF=f.split("\\")[5][0:83] print ("processing_f:"+ FF) GG=g.split("\\")[5][0:85] print ("processing_g:"+ GG) HH=h.split("\\")[5][0:81] print ("processing_h:"+ HH) II=i.split("\\")[5][0:83] print ("processing_i:"+ II) JJ=j.split("\\")[5][0:121] print ("processing_j:"+ JJ) KK=k.split("\\")[5][0:121] print ("processing_k:"+ KK) LL=l.split("\\")[5][0:121] print ("processing_l:"+ LL) ####################################################################### #Cloud Mask scalling Scale_factor1 = float(1.0) add_offset1 = float(0.0) whereClause = "VALUE <= 0" # Execute SetNull outSetNull = SetNull(a, a, whereClause) cldmask=(outSetNull-add_offset1)*Scale_factor1 #AOD scalling Scale_factor2 = float(0.0010000000474974513) add_offset2 = float(0.0) whereClause= "VALUE = -9999" # Execute SetNull outSetNull2 = SetNull(b,b, whereClause) aod=(outSetNull2-add_offset2)*Scale_factor2 aod_1 = Con((aod >= 0.0) & (aod <= 0.1),1) #TOA_B1 Scalling Scale_factor3 = float(0.000053811826) add_offset3 = float(-0.0) whereClause= "VALUE = 65535" # Execute SetNull outSetNull3 = SetNull(c,c, whereClause) toa_b1=(outSetNull3-add_offset3)*Scale_factor3 #TOA_B2 Scalling Scale_factor4 = float(0.00003255546) add_offset4 = float(-0.0) whereClause= "VALUE = 65535" # Execute SetNull outSetNull4 = SetNull(d,d, whereClause) toa_b2=(outSetNull4-add_offset4)*Scale_factor4 #Height Scale_factor5 = float(1.0) add_offset5 = float(-0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull5 = SetNull(e,e, whereClause) height=(outSetNull5-add_offset5)*Scale_factor5 #Sensor Zenith angle Scale_factor6 = float(0.01) add_offset6 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull6 = SetNull(f,f, whereClause) senzen=(outSetNull6-add_offset6)*Scale_factor6 #Sensor azimuth angle Scale_factor7 = float(0.01) add_offset7 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull7 = SetNull(g,g, whereClause) senazm=(outSetNull7-add_offset7)*Scale_factor7 #solar Zenith angle Scale_factor8 = float(0.01) add_offset8 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull8 = SetNull(h,h, whereClause) solzen=(outSetNull8-add_offset8)*Scale_factor8 #solar azimuth angle Scale_factor9 = float(0.01) add_offset9 = float(0.0) whereClause = "VALUE = -32767" # Execute SetNull outSetNull9 = SetNull(i,i, whereClause) solazm=(outSetNull9-add_offset9)*Scale_factor9 # OutRaster = os.path.join(outdir,'MOD04_L22.{0}.tif'.format("solazm")) # solazm.save(OutRaster) #(MODO9)surface reflectance band 1 Scale_factor_10 = float(0.0001) add_offset_10 = float(0.0) whereClause = "VALUE = -28672" # Execute SetNull outSetNull10 = SetNull(j,j, whereClause) sur_b1=(outSetNull10-add_offset_10)*Scale_factor_10 #(MODO9)surface reflectance band 2 Scale_factor_11 = float(0.0001) add_offset_11 = float(0.0) whereClause = "VALUE = -28672" # Execute SetNull outSetNull11 = SetNull(k,k, whereClause) sur_b2=(outSetNull11-add_offset_11)*Scale_factor_11 #(MODO9)surface reflectance band 3 Scale_factor_12 = float(0.0001) add_offset_12 = float(0.0) whereClause = "VALUE = -28672" # Execute SetNull outSetNull12 = SetNull(l,l, whereClause) sur_b3=(outSetNull12-add_offset_12)*Scale_factor_12 # Process: Extract by Cloud_Mask tempEnvironment0 = arcpy.env.cellSize arcpy.env.cellSize = "MAXOF" toa_bb1=arcpy.sa.ExtractByMask(toa_b1, cldmask) toa_bb2=arcpy.sa.ExtractByMask(toa_b2, cldmask) heightt=arcpy.sa.ExtractByMask(height, cldmask) senzenn=arcpy.sa.ExtractByMask(senzen, cldmask) senazmm=arcpy.sa.ExtractByMask(senazm, cldmask) solzenn=arcpy.sa.ExtractByMask(solzen, cldmask) solazmm=arcpy.sa.ExtractByMask(solazm, cldmask) sur_b11=arcpy.sa.ExtractByMask(sur_b1, cldmask) sur_b22=arcpy.sa.ExtractByMask(sur_b2, cldmask) sur_b33=arcpy.sa.ExtractByMask(sur_b3, cldmask) arcpy.env.cellSize = tempEnvironment0 # # Process: Extract by Mask using AOD less than 0.1 value tempEnvironment0 = arcpy.env.cellSize arcpy.env.cellSize = cellsize toa_b1=arcpy.sa.ExtractByMask(toa_bb1, aod_1) toa_b2=arcpy.sa.ExtractByMask(toa_bb2, aod_1) height=arcpy.sa.ExtractByMask(heightt, aod_1) senzen=arcpy.sa.ExtractByMask(senzenn, aod_1) senazm=arcpy.sa.ExtractByMask(senazmm, aod_1) solzen=arcpy.sa.ExtractByMask(solzenn, aod_1) solazm=arcpy.sa.ExtractByMask(solazmm, aod_1) sur_b1=arcpy.sa.ExtractByMask(sur_b11, aod_1) sur_b2=arcpy.sa.ExtractByMask(sur_b22, aod_1) sur_b3=arcpy.sa.ExtractByMask(sur_b33, aod_1) arcpy.env.cellSize = tempEnvironment0 ############################################################################# #Calculate the scattering angle from Resample MOD03 products CosSenZ = np.Cos(senzen * np.pi / 180) CosSolZ = np.Cos(solzen * np.pi / 180) SinSenZ = np.Sin(senzen * np.pi / 180) SinSolZ = np.Sin(solzen * np.pi / 180) RelAzm = np.abs((senazm ) - (solazm)) index = np.where(RelAzm > 180.0) RelAzm[index] = 360.0- RelAzm[index] index = np.where(RelAzm <= 180.0) RelAzm[index] = 180.0- RelAzm[index] CosRe = np.Cos(np.radians(RelAzm)) SctAgl = np.acos((-CosSolZ * CosSenZ) + ((SinSolZ * SinSenZ) * CosRelA)) # RCR for TOA_BAND_1 and TOA_BAND_2 # TOA_BAND_1 B1_ROD = (0.00864 + (0.0000065)) * (0.67)**(-(3.916 + (0.074 * 0.67)+ (0.05/0.67))) Pr1 = (3.0/16.0*np.pi) * (1 + (np.cos(SctAgl) * np.cos(SctAgl))) wr1 = 1.0 RayRef_B1 = (wr1 * B1_ROD * Pr1 ) / (4.0 * (CosSolZ * CosSenZ)) RCR_B1 = toa_b1 - RayRef_B1 OutRaster = os.path.join(outdir,'MOD02HKM_L22.{0}.tif'.format("RAY_CORRECT_TOA_B1")) RCR_B1.save(OutRaster) #TOA_BAND_2 B2_ROD = (0.00864 + (0.0000065)) * (0.86)**(-(3.916 + (0.074 * 0.86)+ (0.05/0.86))) Pr2 = (3.0/16.0*np.pi) * (1 + (np.cos(SctAgl) * np.cos(SctAgl))) wr2 = 1.0 RayRef_B2 = (wr2 * B2_ROD * Pr2 ) / (4.0 * (CosSolZ * CosSenZ)) RCR_B2 = toa_b2 - RayRef_B2 OutRaster = os.path.join(outdir,'MOD02HKM_L22.{0}.tif'.format("RAY_CORRECT_TOA_B2")) RCR_B2.save(OutRaster) # Calculating NDVI float raster ndvi_raster = Divide(Float(Raster(RCR_B2) - Raster(RCR_B1)), Float(Raster(RCR_B2) + Raster(RCR_B1))) OutRaster = os.path.join(outdir,'MOD0_L22.{0}.tif'.format("NDVI")) ndvi_raster.save(OutRaster) print arcpy.GetMessages() print 'finished run: %s\n\n' % (datetime.datetime.now() - start) ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ How do I resolve this matter or have any other way to write this code in the arcpy environment? Dan Patterson‌Joshua Bixby‌Joe Borgione‌Luke Webb‌ ... View more

# #Cloud Mask scalling Scale_factor1 = float(1.0) add_offset1 = float(0.0) setnull1 =arcpy.gp.SetNull_sa(a,a, "in_memory/dat1", "\"Value\" <= 0") ras1=arcpy.Raster(setnull1) cldmask=(ras1-add_offset1)*Scale_factor1 OutRaster1 = os.path.join(outdir,'{0}.img'.format(AA)) arcpy.CopyRaster_management(cldmask,OutRaster1) arcpy.Delete_management("in_memory/dat1") arcpy.Delete_management(ras1)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ In this function how to avoid the in_memory workspace? can you write any example code regarding this code? Xander Bakker‌ ... View more

How can I write this code without using in-memory workspace in this context? One more thing, here I write a function(set-null and some algebra function) for individual files and I write this function multiple time for multiples files(code line 162-270). can I use this function in a single for loop for Individual file?  Xander Bakker‌ ... View more