IDL实现风云三号D星1km数据定标及几何校正
发布日期:2021-05-14 03:07:03 浏览次数:19 分类:精选文章

本文共 5808 字,大约阅读时间需要 19 分钟。

������

  • ������������D���1KM������
  • ������GEO1K������������

������

  • ������250m���������������������������������������������������������������������������������������������1KM���������������
  • ������������GUI���������������������������������������������������������������������������
  • ���������������������������������������������������������������������������
  • ������������������������������������������������������������������������������������������������������������
  • 2020-07-15���������������������������������������������������������������������

������������

  • ������������������IDL���������������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
    1. ������������������������������������������������������
      • ������������FY3D_date_TOA.dat���������19������������TOA������
      • ������������FY3D_date_BT.dat���������5������������������������

      ������������

      ������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 fibClose, 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) endfor 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] endfor 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] EMVI_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  endfunction 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  endfunction 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
    上一篇:pygame教你从0到1一步步实现点到点的智能追踪系统(其一)
    下一篇:东方大观:使用IDL程序制作风云四号的时序动图

    发表评论

    最新留言

    做的很好,不错不错
    [***.243.131.199]2025年04月11日 21时38分29秒