一、介绍
MODIS HDF包含经纬度数据集，因此可以进行几何校正。系统使用的方法是先读取HDF中的经纬度数据集，然后建立51*51的GCP控制点（见图1），通过对GCP点的投影转换来对角度数据集和辐射定标后的科学数据集进行几何校正。在几何校正的过程中将角度数据重采样到1000*1000分辨率，以便后来进行气溶胶反演的时候，按行列号查找角度数据。
二、IDL源码
源码下载：http://pan.baidu.com/s/1dDcXOEP
;+
;:Description:
; Describe the procedure.
;
;:Author: tiands@esrichina.com.cn
;
;:Date: 2013-9-5 11:10:27
;-
;;=======================================================================
;;; 改程序实现对辐射校正后的科学数据几何校正以及对角度数据校正，输出控制点文件
;;;======================================================================
PRO MODIS_REGISTER,filename,fs_name,outname_angle,outname_data,GCP1,GCP2
ENVI, /RESTORE_BASE_SAVE_FILES
ENVI_BATCH_INIT, LOG_FILE = 'batch_log.txt'
COMPILE_OPT idl2
; filename ='H:\MODIS\modis数据\1km\MOD021KM.A2012126.0245.005.2012128183247.hdf'
; fs_name='C:\Users\tiandeshan\Desktop\idl\fushejiaozheng.img'
; outname_angle='C:\Users\tiandeshan\Desktop\idl\angledata.img'
; outname_data='C:\Users\tiandeshan\Desktop\idl\scaledata.img'
; GCP1='C:\Users\tiandeshan\Desktop\idl\GCP1.dat'
; GCP2='C:\Users\tiandeshan\Desktop\idl\GCP2.dat'
;经纬度读取
Lon=READ_DATASET(filename,'Longitude')
Lat=READ_DATASET(filename,'Latitude')
;重采样经纬度
ENVI_WRITE_ENVI_FILE,Lon[*,*],r_fid=fid_Lon,out_name='c:/temp'
FIDlon= MODIS_RESIZE(fid_Lon)
ENVI_FILE_QUERY, FIDlon,dims=dimslon
datalon = ENVI_GET_DATA(fid=FIDlon, dims=dimslon, pos=0)
ENVI_WRITE_ENVI_FILE,Lat[*,*],r_fid=fid_Lat,out_name='c:/temp'
FIDlat= MODIS_RESIZE(fid_Lat)
ENVI_FILE_QUERY, FIDlat,dims=dimslat
datalat = ENVI_GET_DATA(fid=FIDlat, dims=dimslat, pos=0)
;;;=======================================================================
;;; GCP控制点建立
;;;======================================================================
OPENW,lun,GCP1,/Get_lun
PRINTF,lun,' gcp[0,l] ','gcp[1,l] ','gcp[2,l] ','gcp[2,l] ','n ','m '
length=0L
;定义控制点为51*51
gcp= DBLARR(4,2601)
rstr = ["Input File :" , "Output File :"]
ENVI_REPORT_INIT,rstr,title="output the GCP", base=base,/INTERRUPT
FOR i=3.5d,2030,40.53d DO BEGIN
FOR j=3.5d,1354,27.01d DO BEGIN
gcp[0,length]=datalon[j,i]
gcp[1,length]=datalat[j,i]
gcp[2,length]=j
gcp[3,length]=i
PRINTF,lun,gcp[0,length],gcp[1,length],gcp[2,length],gcp[3,length]
length=length+1L
ENDFOR
ENVI_REPORT_STAT, base, i*1354, 2601
ENDFOR
ENVI_REPORT_INIT, base=base, /finish
FREE_LUN,lun
;;;=======================================================================
;;; GCP控制点投影转换
;;;======================================================================
OPENW,lun,GCP2,/Get_lun
PRINTF,lun,'ogcp[0,l]','ogcp[1,l] ','ogcp[2,l]','ogcp[2,l]'
iproj= ENVI_PROJ_CREATE(/geographic)
units = envi_translate_projection_units('Meters')
datum = 'WGS-84'
oproj= ENVI_PROJ_CREATE(/utm, zone=51,datum=datum,units=units)
ENVI_CONVERT_PROJECTION_COORDINATES, gcp[0,*], gcp[1,*],iproj,outx, outy, oproj
; print,gcp[0,*]
ogcp= DBLARR(4,2601)
length=0L
FOR i=3.5d,2030,40.53d DO BEGIN
FOR j=3.5d,1354,27.01d DO BEGIN
ogcp[0,length]=outx[[(i-3.5)/40.53]*51+[(j-3.5)/27.01]]
ogcp[1,length]=outy[[(i-3.5)/40.53]*51+[(j-3.5)/27.01]]
ogcp[2,length]=gcp[2,length]
ogcp[3,length]=gcp[3,length]
PRINTF,lun,ogcp[0,length],ogcp[1,length],ogcp[2,length],ogcp[3,length]
length=length+1L
ENDFOR
ENDFOR
FREE_LUN,lun
;;;=======================================================================
;;; 开始校正
;;;======================================================================
;角度校正
OUT_BNAME=['Warp(SensorZenith)','Warp(SensorAzimuth)',$
'Warp(SolarZenith)','Warp(SolarAzimuth)']
anglefid=ANGLE_DATA(filename)
ENVI_FILE_QUERY,anglefid,nb=nb_angle,dims=dims_angle
REGISTER,anglefid,nb_angle,dims_angle,ogcp,oproj,outname_angle,OUT_BNAME
;科学数据校正
OUT_BNAME=['Warp(1KM Reflectance (band1) [250 Aggr])','Warp(1KM Reflectance (band2) [250 Aggr])',$
'Warp(1KM Reflectance (band3) [500 Aggr])','Warp(1KM Reflectance (band4) [500 Aggr])','Warp(1KM Reflectance (band6) [500 Aggr])',$
'Warp(1KM Reflectance (band7) [500 Aggr])','Warp(1KM Reflectance (band8))','Warp(1KM Reflectance (band19))','Warp(1KM Reflectance (band26))',$
'Warp(1KM Emissive (band29))','Warp(1KM Emissive (band31))','Warp(1KM Emissive (band32))']
;数据读取
IF (ENVI_IS_GDB(fs_name)) THEN BEGIN
ENVI_OPEN_GDB,fs_name,r_fid=data_fid
ENDIF ELSE BEGIN
ENVI_OPEN_FILE,fs_name,r_fid=data_fid
ENDELSE
;envi_open_file, fs_name, r_fid=data_fid
ENVI_FILE_QUERY, data_fid, dims=dims_data, nb=nb_data
REGISTER,data_fid,nb_data,dims_data,ogcp,oproj,outname_data,OUT_BNAME
ENVI_BATCH_EXIT
END
FUNCTION ANGLE_DATA,filename
COMPILE_OPT idl2
;角度信息读取
sensorzenith=READ_DATASET(filename,'SensorZenith')
sensorazimuth=READ_DATASET(filename,'SensorAzimuth')
solarzenith=READ_DATASET(filename,'SolarZenith')
solarazimuth=READ_DATASET(filename,'SolarAzimuth')
;;;=======================================================================
;;; 角度信息合成
;;;======================================================================
ENVI_WRITE_ENVI_FILE,sensorzenith[*,*],r_fid=fid_sez,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,sensorazimuth[*,*],r_fid=fid_sea,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,solarzenith[*,*],r_fid=fid_soz,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,solarazimuth[*,*],r_fid=fid_soa,/IN_MEMORY
ENVI_FILE_QUERY,fid_sez,dims=sezdims
fid_sez_c=ANGLE_CALCULATE(fid_sez)
fid_sea_c=ANGLE_CALCULATE(fid_sea)
fid_soz_c=ANGLE_CALCULATE(fid_soz)
fid_soa_c=ANGLE_CALCULATE(fid_soa)
;将四个有用信息的图像合成一个图像，便于后续的几何校正工作
ENVI_DOIT, 'cf_doit',fid=[fid_sez_c,fid_sea_c,fid_soz_c,fid_soa_c], $
pos=[0,0,0,0], dims=sezdims, $
remove=0,r_fid=info_fid ,/in_memory
;将四个单独的几何信息文件释放
ENVI_FILE_MNG, id=fid_sez, /remove
ENVI_FILE_MNG, id=fid_sea, /remove
ENVI_FILE_MNG, id=fid_soz, /remove
ENVI_FILE_MNG, id=fid_soa, /remove
angle_fid= MODIS_RESIZE(info_fid)
RETURN,angle_fid
END
FUNCTION ANGLE_CALCULATE,fid
COMPILE_OPT idl2
ENVI_FILE_QUERY, fid, dims=dims
t_fid = [fid]
pos = [0]
;表达式
exp='b1*0.0100'
ENVI_DOIT, 'math_doit', $
fid=t_fid, pos=pos, dims=dims,$
exp=exp,r_fid=r_fid, /IN_MEMORY
ENVI_FILE_MNG, id=fid, /remove
RETURN,r_fid
END
FUNCTION READ_DATASET,filename,dataset
SD_id = HDF_SD_START(filename)
index = HDF_SD_NAMETOINDEX(SD_id,dataset)
sds_id=HDF_SD_SELECT(SD_id,index)
HDF_SD_GETINFO,sds_id,DIMS=dim
HDF_SD_GETDATA,sds_id,data
HDF_SD_ENDACCESS,sds_id
RETURN,data
END
FUNCTION MODIS_RESIZE,fid
COMPILE_OPT IDL2
; Open the input file
IF (fid EQ -1) THEN BEGIN
RETURN,-1
ENDIF
ENVI_FILE_QUERY, fid, dims=dims, nb=nb
pos = LINDGEN(nb)
out_name = 'testimg'
;把数据集重采样成1354*2030
ENVI_DOIT, 'resize_doit', $
fid=fid, pos=pos, dims=dims, $
interp=1, rfact=[.20015,.2], $
out_name=out_name, r_fid=r_fid
RETURN,r_fid
END
PRO REGISTER,fid,nb,dims,ogcp,oproj,out_name,OUT_BNAME
COMPILE_OPT idl2
pos = LINDGEN(nb)
pixel_size = [1000.00,1000.00]
ENVI_DOIT, 'envi_register_doit', w_fid=fid, w_pos=pos,w_dims=dims,method=7,$
out_name=out_name,pts=ogcp, proj=oproj,r_fid=r_fid3,pixel_size=pixel_size,OUT_BNAME=OUT_BNAME
ENVI_FILE_MNG, id=r_fid3, /remove
ENVI_FILE_MNG, id=fid, /remove
END
分享：