How To Automatic Calculate SEBAL Model " Landsat 8.0 "
Run in IDL:
pro SEBALAutomatic
; Programmer : Dr. Hamid Rahimi
; Ph.D of Geography (Climatic Hazards)
COMPILE_OPT idl2
e = ENVI(/headless)
IF e EQ !NULL THEN e = ENVI(/headless)
; Landsat Data Import
CD, 'D:\LC08_L1TP_167036_20190821_20190903_01_T1'
File = DIALOG_PICKFILE(TITLE='Import Landsat MTL. File !')
; Landsat File to Raster
Raster = e.OpenRaster(File)
; Spectral Subset
MS = Raster[0]
TIR = Raster[3]
; Spatial Subset Multispectral Data
ULX = 0.0
ULY = 0.0
LRX = 3796
LRY = 3870
Sub_Rect = [ULX,ULY,LRX,LRY]
Out_MS_Sub = DIALOG_PICKFILE(/WRITE, TITLE='Out_MS_Sub !') MS_Sub = ENVITask('SubsetRaster')
MS_Sub.Input_Raster = MS
MS_Sub.Sub_Rect = Sub_Rect
MS_Sub.Output_Raster_Uri = Out_MS_Sub
MS_Sub.Execute
; Spatial Subset Thermal Data
Out_TIR_Sub = DIALOG_PICKFILE(/WRITE, TITLE='Out_TIR_Sub !') TIR_Sub = ENVITask('SubsetRaster')
TIR_Sub.Input_Raster = TIR
TIR_Sub.Sub_Rect = Sub_Rect
TIR_Sub.Output_Raster_Uri = Out_TIR_Sub
TIR_Sub.Execute
; Radiometric Calibration Multispectral Data ; Radiarion Calibration MultiSpectral
Out_MS_Rad = DIALOG_PICKFILE(/WRITE, TITLE='Out_MS_Rad !')
MS_Rad = ENVITask('RadiometricCalibration')
MS_Rad.Input_Raster = MS_Sub.Output_Raster
MS_Rad.Calibration_Type = 'Radiance'
MS_Rad.Output_Raster_Uri = Out_MS_Rad
MS_Rad.Execute
; Reflectance Calibration MultiSpectral
Out_MS_Ref_TOA = DIALOG_PICKFILE(/WRITE, TITLE='Out_MS_Ref_TOA !')
MS_Ref_TOA = ENVITask('RadiometricCalibration')
MS_Ref_TOA.Input_Raster = MS_Sub.Output_Raster
MS_Ref_TOA.Calibration_Type = 'Top-Of-Atmosphere Reflectance'
MS_Ref_TOA.Output_Raster_Uri = Out_MS_Ref_TOA
MS_Ref_TOA.Execute
; Radiometric Calibration Thermal Data ; Radiarion Calibration Thermal
Out_TIR_Rad = DIALOG_PICKFILE(/WRITE, TITLE='Out_TIR_Rad !')
TIR_Rad = ENVITask('RadiometricCalibration')
TIR_Rad.Input_Raster = TIR_Sub.Output_Raster
TIR_Rad.Calibration_Type = 'Radiance'
TIR_Rad.Output_Raster_Uri = Out_TIR_Rad
TIR_Rad.Execute
; Brightness Temperature Thermal
Out_TIR_BT = DIALOG_PICKFILE(/WRITE, TITLE='Out_TIR_BT !')
TIR_BT = ENVITask('RadiometricCalibration')
TIR_BT.Input_Raster = TIR_Sub.Output_Raster
TIR_BT.Calibration_Type = 'Brightness Temperature'
TIR_BT.Output_Raster_Uri = Out_TIR_BT
TIR_BT.Execute
; Atmospheric Correction Multispectral Data ; QUAC MUltispectral Out_MS_Ref_Surf = DIALOG_PICKFILE(/WRITE, TITLE='Out_MS_Ref_Surf !') MS_Ref_Surf = ENVITask('QUAC')
MS_Ref_Surf.Input_Raster = MS_Rad.Output_Raster
MS_Ref_Surf.Sensor = 'Landsat TM/ETM/OLI'
MS_Ref_Surf.Output_Raster_Uri = Out_MS_Ref_Surf
MS_Ref_Surf.Execute
; Spectral Indices : NDVI and LAI
Out_MS_VEG_Ind = DIALOG_PICKFILE(/WRITE, TITLE='Out_MS_VEG_Ind !') MS_VEG_Ind = ENVITask('SpectralIndices')
MS_VEG_Ind.Input_Raster = MS_Ref_Surf.Output_Raster
MS_VEG_Ind.Index = ['Normalized Difference Vegetation Index','Leaf Area Index']
MS_VEG_Ind.Output_Raster_Uri = Out_MS_VEG_Ind
MS_VEG_Ind.Execute
; Emissivity
MS_VEG_Ind = MS_VEG_Ind.Output_Raster
NDVI = MS_VEG_Ind.GetData(Bands = [0])
NDVImin = MIN(NDVI)
NDVImax = MAX(NDVI)
PV = (NDVI-NDVImin)/(NDVImax-NDVImin)
Emissivity = 0.004 * PV + 0.986
Emissivity[WHERE(NDVI LT 0.2,/NULL)] = 0.97
Emissivity[WHERE(NDVI GT 0.5,/NULL)] = 0.99
; Land Surface Temperature
Out_LST_Kel = DIALOG_PICKFILE(/WRITE, TITLE='Out_LST_Kel !')
TIR_BT = TIR_BT.Output_Raster
TIR_BT = TIR_BT.GetData(Bands = [0])
TIR_Rad = TIR_Rad.Output_Raster
TIR_Rad = TIR_Rad.GetData(Bands = [0])
LST_Kel = ( TIR_BT ) / ( 1 + ( TIR_Rad * (TIR_BT / 14380.0 ) * ALOG( Emissivity ) ) )
LST_Cel = LST_Kel - 273.15
Cor = MS.SpatialRef
LSTEx= ENVIRaster(LST_Kel,SPATIALREF=Cor,URI=Out_LST_Kel) LSTEx.SAVE
; Albedo
Out_Albedo_Surf = DIALOG_PICKFILE(/WRITE, TITLE='Out_Albedo_Surf !')
ESUN_1 = 1719.0
ESUN_2 = 1787.0
ESUN_3 = 1746.0
ESUN_4 = 1536.0
ESUN_5 = 997.0
ESUN_6 = 811.0
ESUN_7 = 75.0
ESUN_SUM = ESUN_1 + ESUN_2 + ESUN_3 + ESUN_4 + ESUN_5 + ESUN_6 + ESUN_7
W_1 = ESUN_1 / ESUN_SUM
W_2 = ESUN_2 / ESUN_SUM
W_3 = ESUN_3 / ESUN_SUM
W_4 = ESUN_4 / ESUN_SUM
W_5 = ESUN_5 / ESUN_SUM
W_6 = ESUN_6 / ESUN_SUM
W_7 = ESUN_7 / ESUN_SUM
MS_Ref_TOA = MS_Ref_TOA.Output_Raster
Ref_TOA_1 = MS_Ref_TOA.GetData(Bands = [0])
Ref_TOA_2 = MS_Ref_TOA.GetData(Bands = [1])
Ref_TOA_3 = MS_Ref_TOA.GetData(Bands = [2])
Ref_TOA_4 = MS_Ref_TOA.GetData(Bands = [3])
Ref_TOA_5 = MS_Ref_TOA.GetData(Bands = [4])
Ref_TOA_6 = MS_Ref_TOA.GetData(Bands = [5])
Ref_TOA_7 = MS_Ref_TOA.GetData(Bands = [6])
Albedo_TOA = (Ref_TOA_1 * W_1) + (Ref_TOA_2 * W_2) + (Ref_TOA_3 * W_3) + (Ref_TOA_4 * W_4) + (Ref_TOA_5 * W_5) + (Ref_TOA_6 * W_6) + (Ref_TOA_7 * W_7)
Elevation = 0.0
READ, Elevation, PROMPT= 'Elevation : '
T_sw = 0.75 + (2 * 0.00001) * Elevation
Albedo_Surf = (Albedo_TOA - 0.03 ) / (T_sw ^ 2.0 )
Cor = MS.SpatialRef
Albedo_SurfEx = ENVIRaster(Albedo_Surf,SPATIALREF=Cor,URI=Out_Albedo_Surf)
Albedo_SurfEx.SAVE
; Shortwave Radiation Incoming (RS_Income)
Sun_Elevation = 0.0
READ, Sun_Elevation, PROMPT= 'Sun_Elevation : '
G_sc = 1367.0
Teta = 90 - Sun_Elevation
Teta_d = (Teta * !PI)/180.0
Earth_Sun_Distance = 0.0
READ, Earth_Sun_Distance , PROMPT= 'Earth_Sun_Distance : '
dr = 1 / ( Earth_Sun_Distance ^ 2.0 )
RS_Income = G_sc * COS(Teta_d) * dr * T_sw
; Logwave Radiation Outgoing (RL_Outgoing)
LAI = MS_VEG_Ind.Getdata(Bands = [1])
Emissivity0 = 0.95 + ( 0.01 * LAI )
Emissivity0[WHERE(LAI GE 3.0,/NULL)] = 0.98
Stefan_Boltzman_Cons = 5.67 * 0.00000001
RL_Outgoing = Emissivity0 * Stefan_Boltzman_Cons * ( LST_Kel ^ 4.0 )
; Longwave Radiation Incoming
EmissivityA = 0.85 * ( (-ALOG(T_sw) ) ^ 0.09 )
T_ins = 0.0
READ, T_ins, PROMPT= 'T_ins : '
RL_Incoming = EmissivityA * Stefan_Boltzman_Cons * ( T_ins ^ 4.0 )
; Net Radiation : Rn
Out_Rn = DIALOG_PICKFILE(/WRITE, TITLE='Out_Rn !')
Rn = ( 1 - Albedo_Surf ) * RS_Income + RL_Incoming - RL_Outgoing - ( 1 - Emissivity0 ) * RL_Incoming
Cor = MS.SpatialRef
RnEx = ENVIRaster(Rn,SPATIALREF=Cor,URI=Out_Rn)
RnEx.SAVE
; Soil Heat Flux : G
Out_G = DIALOG_PICKFILE(/WRITE, TITLE='Out_G !')
G_Rn = ( LST_Cel / Albedo_Surf ) * ( ((0.0038 * Albedo_Surf) + (0.0074 * (Albedo_Surf ^ 2.0)) * (1-(0.98 * (NDVI ^ 4.0) ) ) ) )
G = Rn * G_Rn
Cor = MS.SpatialRef
GEx = ENVIRaster(G,SPATIALREF=Cor,URI=Out_G)
GEx.SAVE
; Sensible Heat Flux : H : rah
K = 0.41
Zx = 10.0
Ux = 0.0
READ, Ux , PROMPT= 'Wind_Speed : '
Zom = 0.018 * LAI
U_Star_10m = ( K * Ux ) / ( ALOG( Zx / Zom ) )
U200 = U_Star_10m * ( (ALOG( 200.0 / Zom ) ) / K )
U_Star_200m = ( K * U200 ) / ( ALOG( 200.0 / Zom ) )
Z1 = 0.1
Z2 = 2.0
Rah = ( ALOG(Z2 / Z1) ) / ( U_Star_200m * K )
; ET Refernce : Penman Montith
T_ins_Cel = 0.0
T_min = 0.0
T_max = 0.0
Julian = 0.0
Lat = 0.0
Lon = 0.0
Albedo = 0.0
Humidity = 0.0
Time = 0.0 Lz = 51.43 Sunshine = 0.0 READ, T_ins_Cel , PROMPT= 'T_ins_cel : ' READ, T_min , PROMPT= 'T_min : ' READ, T_max , PROMPT= 'T_max : ' READ, Julian , PROMPT= 'Julian_Date : ' READ, Lat , PROMPT= 'Lat : ' READ, Lon , PROMPT= 'Lon : ' READ, Albedo , PROMPT= 'Albedo : ' READ, Humidity , PROMPT= 'Humidity : ' READ, Time , PROMPT= 'Time : ' READ, Sunshine , PROMPT= 'Sunshine : ' phi = (lat*!PI)/180.0 ; latitude (radian) delta1 = (2503.0*exp((17.27*T_ins_Cel)/(T_ins_Cel+237.3)))/(T_ins_Cel+237.3)^2.0 delta2 = 0.409*sin(((2*!PI*julian)/365.0)-1.39) ; used for hour angle omega_s = ACOS(-TAN(phi)*TAN(delta2)) ; hour angle n = (24.0*omega_s)/(!PI) ; maximum sunshine b = (2*!PI*(julian-81.0))/364.0 ; coefficient of Sc sc = 0.1645*sin(2*b)-0.1255*cos(b)-0.025*sin(b) ; corrective coefficient for solar hour angle w = (!PI/12.0)*((time+(lz-lon)/15.0)+sc-12.0) ; solar hourly angle at half of periode w1 = w-(!PI/24.0) w2 = w+(!PI/24.0) ra = ((12.0*4.92*0.96)/!PI)*((w2-w1)*sin(phi)*sin(delta2)+cos(phi)*cos(delta2)*(sin(w2)-sin(w1))) ; ra is hyper earth radiation rs = (0.25+0.5*(sunshine/n))*ra ; solar radiation (shortwave) rn_s = (1-albedo)*rs ; solar pure radiation (shortwave) e0_t = 0.6108*exp((17.27*T_ins_Cel)/(T_ins_Cel+237.3)) ; temperature dew point e_a = (humidity/100.0)*e0_t ; actual pressure rs_o = (0.75+2*(0.00001)*elevation)*ra ; solar radiation in a clear sky fcd = 1.35*(rs/rs_o)-0.35 ; a part of rn_l computation rn_l = (2.042*0.0000000001)*fcd*(0.34-0.14*(e_a^0.5))*(T_ins_Cel+273.15)^4.0 ; solar pure radiation (longwave) rn_ref = (rn_s)-(rn_l) ; net radiation g_ref = 0.04*rn_ref ; soil heat flux p = 101.3*((293-0.0065*elevation)/293.0)^5.26 ; atmosperic pressure gamma = 0.000665*p ; scycrometer constatnt u2 = (Ux*4.87)/(alog(67.8*10.0-5.42)) e0_t_min = 0.6108*exp((17.27*t_min)/(t_min+237.3)) e0_t_max = 0.6108*exp((17.27*t_max)/(t_max+237.3)) e_s = (e0_t_max+e0_t_min)/2.0 ; saturated vepour pressure cn = 66.0 cd = 0.25 et_reference = ((0.408*delta1*(rn_ref-g_ref)+(gamma*cn)*u2*(e_s-e_a))/(T_ins_Cel+273.15))/$ (delta1+gamma*(1+cd*u2)) ; Latent Heat Flux : Landa Landa = (2.501 - ( 0.00236 * LST_Cel) ) * (10 ^ 6.0) ; H cold and Hot pixels Rn_Cold = 0.0 G_Cold = 0.0 Rn_Hot = 0.0 G_Hot = 0.0 READ, Rn_Cold , PROMPT= 'Rn_Cold : ' READ, G_Cold , PROMPT= 'G_Cold : ' READ, Rn_Hot , PROMPT= 'Rn_Hot : ' READ, G_Hot , PROMPT= 'G_Hot : ' cp = 1004.0 LST_Cold = 0.0 LST_Hot = 0.0 READ, LST_Cold, PROMPT= 'LST_Cold : ' READ, LST_Hot, PROMPT= 'LST_Hot : ' MinElevation = 0.0 READ, MinElevation, PROMPT= 'Min_Elevation :' Gr = 9.81 Temp_Rah = Rah Temp_U_Star = U_Star_200m FOR i = 1, 10 do begin H_Cold = ( Rn_Cold - G_Cold ) - 1.05 * 1000.0 * ( (et_reference * 0.001) / (3600.0)) * Landa H_Hot = Rn_Hot - G_Hot ; DT Computation Air_Pressure = 101.3 * (( (293 - 0.0065 * Elevation)/293.0 ) ^ 5.26) Rho = (1000.0 * Air_Pressure) / (287.0 * 1.01 * LST_Kel) DT_Cold = ( H_Cold * Temp_Rah ) / (Rho * cp) DT_Hot = ( H_Hot * Temp_Rah ) / (Rho * cp) Dif_DT = DT_Hot - DT_Cold Dif_LST = LST_Hot - LST_Cold a = Dif_DT / Dif_LST b = (-a) * LST_Hot + DT_Hot LST_DEM = LST_Kel + 0.0065 * (Elevation - MinElevation) DT_Total = a * LST_DEM + b ; Sensible Heat Flux (H) : Not Corrected H = ( rho * cp ) * ( DT_Total / Temp_Rah ) print, 'Loop' ; L : Monine Abokhov L = (-( ( Rho * cp * (Temp_U_Star ^ 3.0) * LST_Kel ) / ( K * Gr * H ) )) ; Sai Parameters Sai_200m = FLTARR(SIZE(L,/DIMENSIONS)) Sai_2m = Sai_200m Sai_01m = Sai_200m Unstability = WHERE(L LT 0,/NULL, count) IF count NE 0 THEN BEGIN X_200m = (1.0 - 16.0 * (200.0 / L[Unstability]) ) ^ 0.25 Sai_200m[Unstability] = 2 * ALOG( (1 + X_200m) / (2.0) ) + ALOG( (1 + (X_200m ^ 2.0)) / 2.0 ) - 2 * ATAN(X_200m) + 0.5 * !PI X_2m = (1.0 - 16.0 * (2.0 / L[Unstability]) ) ^ 0.25 Sai_2m[Unstability] = 2 * ALOG( (1 + (X_2m ^ 2.0)) / 2.0 ) X_01m = (1.0 - 16.0 * (0.1 / L[Unstability]) ) ^ 0.25 Sai_01m[Unstability] = 2 * ALOG( (1 + (X_01m ^ 2.0)) / 2.0 ) ENDIF Stability = WHERE(L GT 0,/NULL, count) IF count NE 0 THEN BEGIN Sai_200m[Stability] = (-5) * ( 2.0 / L[Stability] ) Sai_2m[Stability] = (-5) * ( 2.0 / L[Stability] ) Sai_01m[Stability] = (-5) * ( 0.1 / L[Stability] ) ENDIF ; Rah Corrected U_Star_200m_Cor = (U200 * K) / (ALOG( 200.0 / Zom ) - Sai_200m) Rah_Cor = ( ALOG(Z2/Z1) - Sai_2m + Sai_01m ) / (U_Star_200m_Cor * K) Temp_Rah = Rah_Cor Temp_U_Star = U_Star_200m_Cor ENDFOR Out_ET_insEx = DIALOG_PICKFILE(/WRITE, TITLE='Out_ET_insEx !') LET = Rn - G - H ET_ins = ( 3600.0 * LET ) / ( (2.501 - 0.00236 * (LST_Kel - 273.15)) * (10 ^ 6.0) ) Cor = MS.SpatialRef ET_insEx = ENVIRaster(ET_ins,SPATIALREF=Cor,URI=Out_ET_insEx) ET_insEx.SAVEEND