IDL 解析葵花8Himawari-8标准数据（HSD），辐射定标、重投影、裁剪_葵花数据投影程序

代码解析全云盘FLDK全波段（16个）的数据并进行定标、重投影和裁剪，重投影用的GLT生成方法参考之前的博客https://blog.csdn.net/qq_33339770/article/details/102957857
输入的是包含全部DAT数据的文件夹
输出的是与NC相同范围的2个影像（1000m和2000m各一景），其中Band3的500m分辨率影像被我采样到1000m

PRO Himawari8
 COMPILE_OPT idl2
  start=systime(1)
  inputfolder='***\HS_H08_20190916_0300_FLDK'
  ;获取影像时间
  basename=file_basename(inputfolder)
  ;获取1km和2km定标后影像数据
  imgdata=H8_Preprocess(inputfolder=inputfolder) 
  R10_data=imgdata['R10']
  R20_data=imgdata['R20']
  ;销毁哈希表
  OBJ_DESTROY,imgdata
  ;应用GTL重投影并裁剪
  outR10=GLTproject(inputdata=R10_data,results=results,timename=timename)
  outR20=GLTproject(inputdata=R20_data,results=results,timename=timename)
  print,systime(1)-start
    
END
;对合成的影像应用GLT几何校正并裁减成NC的范围[80,60,200,-60]
Function GLTproject,inputdata=inputdata,results=results,timename=timename
  COMPILE_OPT idl2
  e=envi(/headless)
  
  inputraster=ENVIRaster(inputdata,URI=e.GetTemporaryFilename(CLEANUP_ON_EXIT='True'))
  inputraster.Save

  bands=inputraster.NBANDS
  cols=inputraster.NCOLUMNS
  if cols eq 11000 then begin
    res='1000'
    gltfile='***\GLT_1000M.dat'
  endif else if cols eq 5500 then begin
    res='2000'
    gltfile='***\GLT_2000M.dat'
  endif
  
  gltraster=e.OpenRaster(gltfile)
  glt_id=ENVIRastertoFID(gltraster)
  ENVI_FILE_QUERY,glt_id,dims=glt_dims
  ;应用GLT校正
  fid=ENVIRasterToFID(inputraster)
  out_name=e.GetTemporaryFilename(CLEANUP_ON_EXIT='True')
  ENVI_DOIT,'ENVI_GEOREF_FROM_GLT_DOIT',FID=fid,GLT_DIMS=glt_dims,GLT_FID=glt_id,OUT_NAME=out_name,pos=lindgen(bands),r_fid=out_fid

  ;裁剪输出
  outraster=ENVIFIDtoRaster(out_fid)
  SpatialRef = outraster.SPATIALREF
  SpatialRef.ConvertLonLatToMap, 80.0, 60.0, MapX, MapY
  SpatialRef.ConvertLonLatToMap, 200.0, -60.0, MapX2, MapY2
  subraster = outraster.Subset(SPATIALREF=SpatialRef,SUB_RECT=[MapX, MapY2, MapX2, MapY])
  subraster.Export,'******.dat','ENVI'
     
  gltraster.Close
  inputraster.Close
  outraster.Close
  subraster.Close

  return,outtiff
END

;--解析数据并将500m重采样至1000m,合成两个数据集--
Function H8_Preprocess,inputfolder=inputfolder
  COMPILE_OPT idl2

  R10_data=make_array(11000,11000,4,/FLOAT)
  R20_data=make_array(5500,5500,12,/FLOAT)
  bandnames=['01','02','03','04','05','06','07','08','09','10','11','12','13','14','15','16']
  ;按波段查找文件
  for i=0,15 do begin
    band_files=file_search(inputfolder,'*B'+bandnames[i]+'_FLDK*.DAT',/test_regular)
    ;必须找到全部10个部分的数据
    if N_ELEMENTS(band_files) EQ 10 THEN BEGIN
      band_data=[]
      seg=['01','02','03','04','05','06','07','08','09','10']
      for j=0,9 do begin
        infile=band_files[j]
        segment=read_hsd(inputfile=infile[0])
        imgdata=segment['pixels']
        imgtime=segment['time'];后续制作大气校正需要的太阳方位角、天顶角使用，代码后续补充
        sun_pos=segment['sun_pos']
        band_data=[[band_data],[imgdata]]
        ;销毁哈希表
        OBJ_DESTROY,segment
      endfor
      if i lt 2 then begin
        R10_data[*,*,i]=band_data
      endif else if i eq 2 then begin
        e=envi(/headless);把B03五百米分辨率重采样
        band_raster=ENVIRaster(band_data,URI=e.GetTemporaryFilename(CLEANUP_ON_EXIT='True'))
        band_raster.Save
        sampleRaster=ENVIResampleRaster(band_raster, DIMENSIONS=[11000,11000]);, METHOD='Bilinear')
        band_data=sampleRaster.GetData()
        sampleRaster.Close
        band_raster.Close
        e.Close
        R10_data[*,*,i]=band_data
      endif else if i lt 4 then begin
        R10_data[*,*,i]=band_data
      endif else begin
        R20_data[*,*,i-4]=band_data
      endelse

    endif
  endfor
  
  outdata=hash('R10',R10_data,'R20',R20_data)
  return,outdata
END

;读取HSD数据并定标为albedo/BT
Function read_hsd,inputfile=inputfile
  COMPILE_OPT idl2
  
  spli=strsplit(file_basename(inputfile),'_')
  ;获取分辨率
  resolution=fix(strmid(file_basename(inputfile),spli[6]+1,2))
  if resolution eq 5 then begin
    cols=22000
    rows=2200
  endif else if resolution eq 10 then begin
    cols=11000
    rows=1100
  endif else if resolution eq 20 then begin
    cols=5500
    rows=550
  endif
  
  openr,h8_lun,inputfile,/get_lun
  ;获取Observation start time
  point_lun,h8_lun,46
  imgtime=dblarr(1);R8
  readu,h8_lun,imgtime
  ;获取Total header length
  point_lun,h8_lun,70
  header_length=uintarr(1);I4
  readu,h8_lun,header_length
  ;获取影像
  point_lun,h8_lun,header_length
  imgdata=uintarr(cols,rows)
  readu,h8_lun,imgdata 
  ;获取Sun's position
  point_lun,h8_lun,510
  sun_pos=dblarr(3);R8
  readu,h8_lun,sun_pos
  ;获取Band number
  point_lun,h8_lun,601
  band=uintarr(1);I2
  readu,h8_lun,band  
  ;获取Gain
  point_lun,h8_lun,617
  Gain=dblarr(1);R8
  readu,h8_lun,Gain
  ;获取Offset
  point_lun,h8_lun,625
  Offset=dblarr(1);R8
  readu,h8_lun,Offset
  ;计算radiance
  radiance=imgdata*Gain[0]+Offset[0]
  
  if band le 6 then begin;计算反射率
    ;获取radiance to albedo
    point_lun,h8_lun,633
    cc=dblarr(1);R8
    readu,h8_lun,cc
    
    ;计算反射率
    albedo=radiance*cc[0]   
    ;返回值outdata
    outdata=albedo
  endif else begin;计算亮温
    ;获取Central wave length
    point_lun,h8_lun,603
    wv=dblarr(1);R8
    readu,h8_lun,wv
    ;获取radiance to brightness temperature(c0)
    point_lun,h8_lun,633
    c0=dblarr(1);R8
    readu,h8_lun,c0
    ;获取radiance to brightness temperature(c1)
    point_lun,h8_lun,641
    c1=dblarr(1);R8
    readu,h8_lun,c1
    ;获取radiance to brightness temperature(c2)
    point_lun,h8_lun,649
    c2=dblarr(1);R8
    readu,h8_lun,c2
    ;获取Speed of light (c)
    point_lun,h8_lun,681
    c=dblarr(1);R8
    readu,h8_lun,c    
    ;获取Planck constant (h)
    point_lun,h8_lun,689
    h=dblarr(1);R8
    readu,h8_lun,h
    ;获取Boltzmann constant(k)
    point_lun,h8_lun,697
    k=dblarr(1);R8
    readu,h8_lun,k
    
    ;计算亮温
    wv=wv[0]*1e-6
    rad=radiance*1e6
    Te=h[0]*c[0]/k[0]/wv[0]/(ALOG(2*h[0]*c[0]^2/(wv[0]^5*rad)+1))
    BT=c0[0]+c1[0]*Te+c2[0]*Te^2
    
    ;返回值 outdata
    outdata=BT
  endelse
  free_lun,h8_lun 
  ;返回：albedo/BT,时间,太阳坐标
  out=hash('pixels',outdata,'time',imgtime,'sun_pos',sun_pos)
  return,out;data
  
END



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
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217