IDL实现风云三号D星1km数据定标及几何校正

数据

• 风云三号D星1KM数据
• 对应GEO1K定位数据

说明

• 因为250m的定位数据缺少角度信息，不能定标至大气层顶反射率，所以只做了1KM
• 暂时没做GUI和批量功能，等到明年年初再看情况做
• 关于发射率，官方说明是不需要定标，只需要通过slope和intercept线性转换就行
• 关于亮温的计算，单位给我整懵了，所以直接参考了下面这篇硕士论文里的计算公式
• 2020-07-15更改：感谢易智瑞的杜老师能够发现代码中的定标问题，现已修改。

基于多源卫星数据的秸秆焚烧监测研究

使用

1. 编译
2. 命令行键入如下代码，其中fnDat为数据文件，fnLoc为定位文件

fnDat = 'F:\FY3D_MERSI_GBAL_L1_date_1000M_MS.HDF'
fnDat = 'F:\FY3D_MERSI_GBAL_L1_date_GEO1K_MS.HDF'
preFY_1K, fnDat, fnLoc

3. 数据所在文件夹生成两个文件，形如FY3D_date_TOA.dat为定标的19谱段TOA数据，形如FY3D_date_BT.dat为5谱段的亮温数据。

源码

pro preFY_1K, fnDat, fnLoc
  compile_opt idl2, hidden
  tic
  DLM_LOAD, 'HDF5', 'XML', 'MAP_PE', 'NATIVE'
  e = ENVI(/h)

  dirname = FILE_DIRNAME(fnDat)
  date = (FILE_BASENAME(fnDat)).Extract('20[1-4][0-9][0,1][0-9][0-3][0-9]_[0-9][0-9][0-9][0-9]')
  fnRef = dirname + PATH_SEP() + 'FY3D_' + date + '_TOA.dat'
  fnBT = dirname + PATH_SEP() + 'FY3D_' + date + '_BT.dat'
  
  fdid = H5F_OPEN(fnDat)
    rawRefAggr = h5GetData(fdid, '/Data/EV_250_Aggr.1KM_RefSB')
    rawRef = h5GetData(fdid, '/Data/EV_1KM_RefSB')
    rawEmsAggr = h5GetData(fdid, '/Data/EV_250_Aggr.1KM_Emissive')
    rawEms = h5GetData(fdid, '/Data/EV_1KM_Emissive')
    calCoef = h5GetData(fdid, '/Calibration/VIS_Cal_Coeff')

    ESUN = h5GetAttr(fdid, 'Solar_Irradiance')
    dst = (h5GetAttr(fdid, 'EarthSun Distance Ratio'))[0]
    tbbA = h5GetAttr(fdid, 'TBB_Trans_Coefficient_A')
    tbbB = h5GetAttr(fdid, 'TBB_Trans_Coefficient_B')
  H5F_CLOSE, fdid

  flid = H5F_OPEN(fnLoc)
    lat = h5GetData(flid, '/Geolocation/Latitude')
    lon = h5GetData(flid, '/Geolocation/Longitude')
    slrZth = h5GetData(flid, '/Geolocation/SolarZenith') * 0.01 * !pi / 180E
  H5F_CLOSE, flid

  refChnl = MAKE_ARRAY([2048, 2000, 19], type = 5)
  refChnl[*, *, 0 : 3] = rawRefAggr
  refChnl[*, *, 4 : 18] = rawRef

  emsChnl = MAKE_ARRAY([2048, 2000, 6], type = 5)
  emsChnl[*, *, 0 : 3] = rawEms
  emsChnl[*, *, 4 : 5] = rawEmsAggr

  k = (dst ^ 2)
  for i = 0, 18 do refChnl[*, *, i] = k * $
    (calCoef[0, i] + calCoef[1, i] * refChnl[*, *, i] $
    + calCoef[2, i] * refChnl[*, *, i] ^ 2) $
    /  COS(slrZth)
  
  wn = [2634.359, 2471.654, 1382.621, 1168.182, 933.364, 836.941]
  for i = 0, 5 do emsChnl[*, *, i] = $
    tbbA[i] * (1.43878 * wn[i] / ALOG(1 + $
    (1.191E-5 * wn[i] ^ 3 / emsChnl[*, *, i]))) + tbbB[i]
  
  refTempFn = e.GetTemporaryFilename()
  ref = ENVIRaster(refChnl, $
    uri = refTempFn, interleave = 'BSQ')
  ref.Save
  ref_id = ENVIRasterToFID(ref)
  ENVI_FILE_QUERY,ref_id, dims = ref_dims
  
  emsTempFn = e.GetTemporaryFilename()
  ems = ENVIRaster(emsChnl, $
    uri = emsTempFn, interleave = 'BSQ')
  ems.Save
  ems_id = ENVIRasterToFID(ems)
  
  latTempFn = e.GetTemporaryFilename()
  latRaster = ENVIRaster(lat, uri = latTempFn)
  latRaster.Save
  lat_id = ENVIRasterToFID(latRaster)
  
  lonTempFn = e.GetTemporaryFilename()
  lonRaster = ENVIRaster(lon, uri = lonTempFn)
  lonRaster.Save
  lon_id = ENVIRasterToFID(lonRaster)
  
  proj = ENVI_PROJ_CREATE(/GEOGRAPHIC)
  gltTempFn = e.GetTemporaryFilename()
  ENVI_DOIT, 'ENVI_GLT_DOIT', DIMS = ref_dims, $
    I_PROJ = proj,  O_PROJ = proj, $
    OUT_NAME = gltTempFn, $
    R_FID = glt_id, ROTATION = 0 , $
    X_FID = lon_id, X_POS = [0], $
    Y_FID = lat_id, Y_POS = [0]
  
  ENVI_DOIT, 'ENVI_GEOREF_FROM_GLT_DOIT', $
    BACKGROUND = 0E, FID = ref_id , $
    GLT_FID = glt_id, R_FID = outR_fid, $
    OUT_NAME = fnRef, POS = [0 : 18]
  
  ENVI_DOIT, 'ENVI_GEOREF_FROM_GLT_DOIT', $
    BACKGROUND = 0E, FID = ems_id , $
    GLT_FID = glt_id, R_FID = outB_fid, $
    OUT_NAME = fnBT, POS = [0 : 5]
  
  CD, FILE_DIRNAME(e.GetTemporaryFilename())
  ref.Close & ems.Close & latRaster.Close & lonRaster.Close
  FILE_DELETE, refTempFn, emsTempFn, latTempFn, lonTempFn, $
    (FILE_BASENAME(refTempFn)).Replace('dat', 'hdr'), $
    (FILE_BASENAME(emsTempFn)).Replace('dat', 'hdr'), $
    (FILE_BASENAME(latTempFn)).Replace('dat', 'hdr'), $
    (FILE_BASENAME(lonTempFn)).Replace('dat', 'hdr'), $
    /ALLOW_NONEXISTENT, /QUIET
  ENVI_FILE_MNG, id = glt_id, /d, /r
  
  refRaster = e.OpenRaster(fnRef)
  refRaster.Metadata.UpdateItem, 'band names', $
    ['Ref_1', 'Ref_2', 'Ref_3(Aggr_1)', 'Ref_4',' Ref_5(Aggr_2)', $
    'Ref_6', 'Ref_7(Aggr_3)', 'Ref_8', 'Ref_9', 'Ref_10', $
    'Ref_11', 'Ref_12(Aggr_4)', 'Ref_13', 'Ref_14', $
    'Ref_15', 'Ref_16', 'Ref_17', 'Ref_18', 'Ref_19']
  refRaster.Metadata.UpdateItem, 'wavelength units', 'Micrometers'
  refRaster.Metadata.AddItem, 'wavelength', $
    [0.412, 0.443, 0.47, 0.49, 0.55, 0.555, 0.65, 0.67, $
    0.709, 0.746, 0.865, 0.865, 0.905, 0.936, 0.94, $
    1.24, 1.38, 1.64, 2.13]
  refRaster.WriteMetadata
  refRaster.Close
  
  btRaster = e.OpenRaster(fnBT)
  btRaster.Metadata.UpdateItem, 'band names', $
    ['Ems_1', 'Ems_2', 'Ems_3', 'Ems_4', $
    'Ems_5(Aggr_1)', 'Ems_6(Aggr_2)']
  btRaster.Metadata.UpdateItem, 'wavelength units', 'Micrometers'
  btRaster.Metadata.AddItem, 'wavelength', $
    [3.8, 4.05, 7.2, 8.55, 10.8, 12]
  btRaster.WriteMetadata
  btRaster.Close
  
  ENVI_FILE_MNG, id = outR_fid, /r
  ENVI_FILE_MNG, id = outB_fid, /r
  
  e.Close
  toc
end

function h5GetData, fid, str
  compile_opt idl2, hidden

  str_id = H5D_OPEN(fid, str)
    slope = FLOAT(h5GetAttr(str_id, 'Slope'))
    intercept = FLOAT(h5GetAttr(str_id, 'Intercept'))
    data = FLOAT(H5D_READ(str_id))
    foreach _slope, slope, index do $
      data[*, *, index] = _slope * data[*, *, index] + intercept[index]
  H5D_CLOSE, str_id

  RETURN, data
end

function h5GetAttr, fid, str
  compile_opt idl2, hidden

  str_id = H5A_OPEN_NAME(fid, str)
    attr = H5A_READ(str_id)
  H5A_CLOSE, str_id

  RETURN, attr
end