;;=======================================================================
;;; 改程序实现对辐射校正后的科学数据几何校正以及对角度数据校正，输出控制点文件
;;;======================================================================
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
