FY-3D如何预处理，做NDVI出图

FY-3D相关数据及文件下载位置：风云卫星遥感数据服务网
以MERSI传感器为例。文档中有它的定标方法。

定标方法：

大致是通过原来先求dn值再求反射率数据。
其中Slope、Intersept都是以数据的一个属性出现的：

对应的Cal0、1、2等在另外的数据集中：Calibration中

对应的数据如下：

按照列分别为式中的Cal0、Cal1、Cal2

一般常用的波段应该是前四个波段，就下面这几个：

对应的Cal系数就是0-3行的，具体的对应也可以在前面的帮助文档中看到。
那么我们看下怎么读取FY3D数据并给他做个简单的NDVI处理吧！
首先，该数据后缀是.HDF, 按理来说在IDL语言中用hdf4可以打开，但是我在打开过程中发现这样根本读取不了，在抱着试试看的心理，我用hdf5的方式给他打开了。下面附上我IDL处理的代码

;author:Lf Zhang
;2021-05-25
function get_hdf5_data,hd_name,filename;filename是数据集的具体路径
  hd_id=h5f_open(hd_name)
  sd_id=h5d_open(hd_id,filename)
  data=h5d_read(sd_id)
  return,data
  h5d_close,sd_id
  h5f_close,hd_id
end
function get_hdf5_att_data,hd_name,filename,sds_name;filename是数据集的具体路径
  hd_id=h5f_open(hd_name)
  sd_id=h5d_open(hd_id,filename)
  sds_id=h5a_open_name(sd_id,sds_name)
  data=h5a_read(sds_id)
  return,data
  h5d_close,sd_id
  h5f_close,hd_id
end
pro fy3d_glt
  dir="H:\fy3D\"
  file_list=file_search(dir,'*250M_MS.HDF',count=file_n)
  for file_i=0,file_n-1 do begin
    band_nir=get_hdf5_data(file_list[file_i],'/Data/EV_250_RefSB_b4/')
    band_red=get_hdf5_data(file_list[file_i],'/Data/EV_250_RefSB_b3/')
    nir_slope=get_hdf5_att_data(file_list[file_i],'/Data/EV_250_RefSB_b4/','Slope')
    red_slope=get_hdf5_att_data(file_list[file_i],'/Data/EV_250_RefSB_b3/','Slope')
    nir_intercept=get_hdf5_att_data(file_list[file_i],'/Data/EV_250_RefSB_b4/','Intercept')
    red_intercept=get_hdf5_att_data(file_list[file_i],'/Data/EV_250_RefSB_b3/','Intercept')
    cal=get_hdf5_data(file_list[file_i],'/Calibration/VIS_Cal_Coeff/')
    lat=get_hdf5_data(file_list[file_i],'/Geolocation/Latitude/')
    lon=get_hdf5_data(file_list[file_i],'/Geolocation/Longitude/')
    size_data=size(band_nir)
    size_col=size_data[1]
    size_row=size_data[2]
    lat=congrid(lat,size_col,size_row,/interp);此处由于它的经纬度与数据格式不一样，如果后边做GLT是需要先给他插值成一样大小的
    lon=congrid(lon,size_col,size_row,/interp)
    lon_min=min(lon)
    lat_max=max(lat)
    cal_0=cal[0,*]
    cal_1=cal[1,*]
    cal_2=cal[2,*]
    ;print,cal_0
    nir_cal_0=cal_0[3]
    red_cal_0=cal_0[2]
    nir_cal_1=cal_1[3]
    red_cal_1=cal_1[2]
    nir_cal_2=cal_2[3]
    red_cal_2=cal_2[2]
    dn_nir=((band_nir ge 0) and (band_nir le 4095))*band_nir*nir_slope[0]+nir_intercept[0]
    ref_nir=nir_cal_2*dn_nir^2.0+nir_cal_1*dn_nir+nir_cal_0
    
    dn_red=((band_red ge 0) and (band_red le 4095))*band_red*red_slope[0]+red_intercept[0]
    ref_red=red_cal_2*dn_red^2.0+red_cal_1*dn_red+red_cal_0
    
    ndvi=float(ref_nir-ref_red)/(ref_nir+ref_red)
    write_tiff,dir+file_basename(file_list[file_i],'.HDF')+'_ndvi.tiff',ndvi_result,/float
endfor
end

在上面我没有重投影，如果需要可以自己再输出一下经纬度，然后再ENVI中BUILD GLT操作就可以啦，然后未投影的效果图如下：

GLT后的效果如下：

因为可能有有的数值超过范围，可以自己在输出之前限制下数值范围，不在正常值范围的直接舍去：

ndvit=((ndvi le 1)and(ndvi ge -1))*ndvi

加到出图前就行啦！