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
加到出图前就行啦!