-
Notifications
You must be signed in to change notification settings - Fork 9
/
Copy pathtest_albedo.pro
101 lines (93 loc) · 6.65 KB
/
test_albedo.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
;+
;author:shaodonghang
;date:2016-4-22
;-
pro TEST_Albedo
;根据MOD09GA反射率和角度计算反照率
;打开反射率及角度数据
fn_reflectance=dialog_pickfile(title='选择反射率数据', get_path=work_dir)
cd, work_dir
;fn_angle=dialog_pickfile(title='选择角度数据')
;读入数据
envi_open_file, fn_reflectance, r_fid=fid_reflectance
;envi_open_file, fn_angle, r_fid=fid_angle
envi_file_query,fid_reflectance,ns=ns,nl=nl,nb=nb,dims=dims,$
data_type=data_type,interleave=interleave,offset=offset
map_info=envi_get_map_info(fid=fid_reflectance)
;for i=0, nb-1 do begin
;读取反射率数据
ref_b1=envi_get_data(fid=fid_reflectance,dims=dims,pos=0)
; ref_b2=envi_get_data(fid=fid_reflectance,dims=dims,pos=1)
; ref_b3=envi_get_data(fid=fid_reflectance,dims=dims,pos=2)
; ref_b4=envi_get_data(fid=fid_reflectance,dims=dims,pos=3)
; ref_b5=envi_get_data(fid=fid_reflectance,dims=dims,pos=4)
; ref_b6=envi_get_data(fid=fid_reflectance,dims=dims,pos=5)
; ref_b7=envi_get_data(fid=fid_reflectance,dims=dims,pos=6)
;读取角度数据
angle1=envi_get_data(fid=fid_reflectance,dims=dims,pos=7)
angle2=envi_get_data(fid=fid_reflectance,dims=dims,pos=8)
angle3=envi_get_data(fid=fid_reflectance,dims=dims,pos=9)
;调用ART_BRDF函数计算Albedo
; Albedo_result_b1=ART_BRDF(angle1,angle2,angle3,ref_b1)
; Albedo_result_b2=ART_BRDF(angle1,angle2,angle3,ref_b2)
; Albedo_result_b3=ART_BRDF(angle1,angle2,angle3,ref_b3)
; Albedo_result_b4=ART_BRDF(angle1,angle2,angle3,ref_b4)
; Albedo_result_b5=ART_BRDF(angle1,angle2,angle3,ref_b5)
; Albedo_result_b6=ART_BRDF(angle1,angle2,angle3,ref_b6)
; Albedo_result_b7=ART_BRDF(angle1,angle2,angle3,ref_b7)
A=1.247
B=1.186
C=5.157
;
;
u0=cos(angle1*!DTOR)
u=cos(angle2*!DTOR)
;
;cos(a*3.14/180)
s0=sin(angle1*!DTOR)
s=sin(angle2*!DTOR)
;
; sita = acos(-u*u0+s*s0*cos(Relative*!DTOR))
; P = 11.1*exp(-0.087*sita*180.0/!DPI)+1.1*exp(-0.014*sita*180.00/!DPI)
; R0 = ((A+B*(u+u0)+C*u*u0+P)*1.0)/(4.0*(u+u0))
; print,R0,P,sita,u,u0
; ku=3*(1+2*u)/7.00
; ku0=3*(1+2*u0)/7.00
; v=ku*ku0/R0
; print,ku,ku0,v,1/v
;
angle4=acos(-u*u0+s*s0*cos((180.0-angle3)*!DTOR))
;从P开始，结果出现问题
P=11.1*exp(-0.087*angle4*180.0/!DPI)+1.1*exp(-0.014*angle4*180.00/!DPI)
;print,P
R0 = ((A+B*(u+u0)+C*u*u0+P)*1.0)/(4.0*(u+u0))
z1=max(R0)
z2=min(R0)
z3=float(mean(R0))
print,z1,z2,z3
uu0=3*(1+2*u0)/7.00
uu=3*(1+2*u)/7.00
;print,uu0,uu
;print,R0
;R1=0.97570502
;R1=1.0400882
R2=float(R0)
Rs=(ref_b1/R2)^(R2/(uu0*uu)) ;Rs为积雪光谱白空反照率
Rp=Rs^uu0 ;Rp为积雪光谱黑空反照率
f=0.3 ;f为天空散射光因子
result=(1-f)*Rp+f*Rs ;result位积雪光谱反照率，f为天空散射光因子，取地表反照率产品MOD43中的0.3
;PRINT,result
;RETURN,result
;endfor
;窄波段反照率向宽波段反照率转换;;窄波段反照率向宽波段反照率转换，A=-0.0093+0.157a1+0.2789a2+0.3829a3+0.1131a5+0.069a7
;ALBEDO=-0.0093+0.157*ref_b1+0.2789*ref_b2+0.3829*ref_b3+0.1131*ref_b5+0.069*ref_b7
;保存反照率结果
; file=FILE_BASENAME(fn)
; filetime=strpos(file,'A+1',/reverse_search)
; filename=strmid(file,filetime,9)
; out_name = imgpath+filename+'_Albedo'+StrTrim(i+1,2)+'.tif'
o_fn=dialog_pickfile(title='反照率结果保存为')
ENVI_WRITE_ENVI_FILE, result, OUT_NAME=o_fn, /NO_COPY, $
NS=NS, NL=NL, NB=1, DATA_TYPE=4, INTERLEAVE=INTERLEAVE, $
OFFSET=OFFSET, MAP_INFO=MAP_INFO
END