HANIMANI1

Raster Data analysis with trigonometry function in arcpy. Giving Attribute Error: 'Raster'

Discussion created by HANIMANI1 on Jul 23, 2019
Latest reply on Jul 25, 2019 by xander_bakker

Previously, I discussed this problem in this platform, but I am unable to short out this problem as per discussion directive.Right now, I am specifying the persisting problem for further solutions.

I am trying to write simple trigonometric calculation (code line 297-324) using numpy but it gave me one error.

Traceback (most recent call last):
 File "F:\DB_test_data\python_script\FINAL_DB.py", line 297, in <module>
cosSenZ = np.cos(senzen * np.pi / 180)
AttributeError: 'Raster' object has no attribute 'cos'

Here uploading MY_CODE( trigonometric calculation(code line 297-324)) and test data for further checking. Please look at this matter.

#-------------------------------------------------------------------------------
# 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
    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]
    cosRelA = np.cos(np.radians(RelAzm))

    SctAgl =  np.arccos((-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*np.pi) * (1 + (np.cos(SctAgl) * np.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)




print arcpy.GetMessages()
print 'finished run: %s\n\n' % (datetime.datetime.now() - start)

Dan PattersonXander BakkerJoe BorgioneJoshua BixbyLuke Webb

Attachments

Outcomes