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

BIJOYGAYEN · ‎07-23-2019

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 Patterson‌Xander Bakker‌Joe Borgione‌Joshua Bixby‌Luke Webb‌

CharlesPreppernau · ‎07-23-2019

Hello,

It looks like the issue may be that np.cos, a NumPy function, is being called on an ArcPy Raster object. My understanding is that ArcPy rasters need to be converted to Numpy arrays with arcpy.RasterToNumPyArray() before this will work. You can convert back with arcpy.NumPyArrayToRaster().

Hope this helps,

-Charles

BIJOYGAYEN · ‎07-23-2019

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.

CharlesPreppernau · ‎07-23-2019

I don't know of any alternatives; I've always converted rasters to arrays. I find it helpful to use something like this pattern:

inputRaster = arcpy.Raster(inputPath)

cX = inputRaster.meanCellWidth

cY = inputRaster.meanCellHeight

crs=inputRaster.spatialReference

anchorX = inputRaster.extent.XMin

anchorY = inputRaster.extent.YMin

corner = arcpy.Point(anchorX, anchorY)

inputArray = arcpy.RasterToNumPyArray(inputRaster)

outputArray = DoAllYourNumpyStuff(inputArray)

outputRaster = arcpy.NumPyArrayToRaster(outputArray,corner,cX,cY)

arcpy.DefineProjection_managment(outputRaster,crs)

XanderBakker · ‎07-23-2019

Hi BIJOY GAYEN ,

In your other thread (https://community.esri.com/thread/235537-runtimeerror-error-999998-unexpected-error-in-set-null-func...) I already suggested to convert the rasters to numpy arrays as required to be able to use numpy. The errors that you have seen are probably caused by invalid values in your array. This could be caused by the NoData values in your rasters.

As Charles Preppernau mentioned, you can't use an ArcGIS Raster object directly with numpy. You have to use RasterToNumPyArray to convert your raster to a numpy array and than use it. This you already did in the other thread, but than you ran into problems with invalid values for the numpy functions.

I also suggested to try and use the SA function. Can you try and change the lines:

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

... into something like this:

    cosSenZ = arcpy.sa.Cos(senzen * math.pi / 180.0)
    cosSolZ = arcpy.sa.Cos(solzen * math.pi / 180.0)
    sinSenZ = arcpy.sa.Sin(senzen * math.pi / 180.0)
    sinSolZ = arcpy.sa.Sin(solzen * math.pi / 180.0)‍‍‍‍

Does that work?

BIJOYGAYEN · ‎07-24-2019

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.

BIJOYGAYEN · ‎07-25-2019

Hello Xander Bakker‌ sir have you check my problem?

XanderBakker · ‎07-25-2019

Hi BIJOY GAYEN ,

I haven't tested it,but there are a couple of things I can see that won't work.

You have a line of code where you are using arcpy.sa.where. This does not exists (as far as I know). You should use a Con statement instead. See http://desktop.arcgis.com/en/arcmap/latest/tools/spatial-analyst-toolbox/con-.htm
You create a raster called index, but you cannot use this raster as an index. This would be another Con statement to determine what to do with those pixels that meet the condition.
You are still using raster object in numpy function. When you want to use numpy, please convert the raster to numpy array first or find similar functionality available in arcpy.sa.