本文以MOD07数据为例,利用IDL程序进行几何校正,拼接和剪切处理。
1、Inputs.pro
calcu,2008144,2008233,'G:\tangronglin\mod07\','MOD07_L2.A'
2、calcu.pro
pro calcu,startday,endday,path,prefix
a=systime(1)
daystart=long(startday)
dayend=long(endday)
while(daystart le dayend) do begin
;suffix='Retrieved_Temperature_Profile'
suffix='Retrieved_Moisture_Profile'
arg1=strtrim(string(daystart),2)
time=100
while (time le 800) do begin
timestr=strtrim(string(time),2)
file=file_search(path+prefix+arg1+'.0'+timestr+'.005'+'.*.hdf',COUNT=count1)
if count1 ne 0 then begin
MODIS_georeference,path,arg1,prefix,timestr,suffix,file[0]
endif
time=time+1
endwhile
;-----------------------------------------------------------
file=file_search(path+prefix+arg1+'.*_'+suffix+'.tif',COUNT=count1);to be modified
;-----------------------------------------------------------
if count1 ne 0 then begin
fill_value=-32768
example_resize_doit,arg1,path,prefix,suffix,fill_value
endif
daystart=daystart+1
endwhile
b=systime(1)
print,b-a
end
stop
3、MODIS_georeference.pro
pro MODIS_georeference,path,arg1,prefix,timestr,suffix,file
envi, /restore_base_save_files
envi_batch_init, log_file='c:/batch.txt'
SD_id = HDF_SD_START(file)
index = HDF_SD_NAMETOINDEX(SD_id, suffix)
sds_id=HDF_SD_SELECT(SD_id,index)
HDF_SD_GETINFO,sds_id,DIMS=dim
ns=dim[0]
nl=dim[1]
print,ns,nl
hdf_sd_getdata,sds_id,b
HDF_SD_END, SD_id
write_tiff,'refb0.tif',b[*,*,16],/float
;;==========================================================
; map_temp=bytarr(ns,nl)
; print,ns,nl
rstr = ["Input File :" , "Output File :"]
envi_report_init,rstr,title="map", base=base,/INTERRUPT
envi_report_init, base=base, /finish
ENVI_WRITE_ENVI_FILE,b[*,*,16],r_fid=fid_map,out_name='c:/temp'
;;==========================================================
; ;read coordinate
;;==========================================================
SD_id = HDF_SD_START(file)
index1 = HDF_SD_NAMETOINDEX(SD_id, 'Longitude')
sds_id1=HDF_SD_SELECT(SD_id,index1)
HDF_SD_GETDATA,sds_id1,Lon
HDF_SD_ENDACCESS,sds_id1
index2 = HDF_SD_NAMETOINDEX(SD_id, 'Latitude')
sds_id2=HDF_SD_SELECT(SD_id,index2)
HDF_SD_GETDATA,sds_id2,Lat
HDF_SD_ENDACCESS,sds_id2
HDF_SD_END,SD_id
write_tiff,'Lat.tif',Lat,/float
;;;=======================================================================
;;; georeference data
;;;======================================================================
openw,lun,'GCP1.dat',/Get_lun
printf,lun,' gcp[0,l] ','gcp[1,l] ','gcp[2,l] ','gcp[2,l] ','n ','m '
; Numl=50
Numl=nl
; Nums=50
Nums=ns
; s_space=float(ns-1)/float(Numl-1)
; l_space=float(nl-1)/float(Nums-1)
; s_space=float(ns)/float(Numl)
; l_space=float(nl)/float(Nums)
l=0L
gcp= DBLARR(4,Nums*Numl)
rstr = ["Input File :" , "Output File :"]
envi_report_init,rstr,title="output the GCP", base=base,/INTERRUPT
for i=0,Numl-1 do begin
; m=LONG64(i*l_space/5)
; m=LONG64(i*l_space)
for j=0,Nums-1 do begin
; n=LONG64(j*s_space/5)
; n=LONG64(j*s_space)
gcp[0,l]=Lon[j,i];Lon[n,m];lon[LONG64(j*s_space)]
; print,gcp[0,l]
gcp[1,l]=Lat[j,i];Lat[n,m];Lat[LONG64(i*l_space)]
; print,gcp[1,l]
gcp[2,l]=j;LONG64(j*s_space)+1
; print,gcp[2,l]
gcp[3,l]=i;LONG64(i*l_space)+1
; print,gcp[3,l]
printf,lun,gcp[0,l],gcp[1,l],gcp[2,l],gcp[3,l]
l=l+1L
endfor
envi_report_stat, base, i*Nums, Nums*Numl
endfor
envi_report_init, base=base, /finish
free_lun,lun
;;===============================================
;; Projection
;;===============================================
openw,lun,'GCP2.dat',/Get_lun
printf,lun,' ogcp[0,l] ','ogcp[1,l] ','ogcp[2,l] ','ogcp[2,l]'
iproj= ENVI_PROJ_CREATE(/geographic)
params = [6378245.0, 6356863.0,0.000000,105.000000,0.0, 0.0,25.0,47.0]
datum = 'WGS 84'
name = 'China'
units = ENVI_TRANSLATE_PROJECTION_UNITS('Meters')
oproj = ENVI_PROJ_CREATE(type=9,name=name, datum=datum, params=params,units=units)
ENVI_CONVERT_PROJECTION_COORDINATES, gcp[0,*], gcp[1,*],iproj,outx, outy, oproj
; print,gcp[0,*]
ogcp= DBLARR(4,Nums*Numl)
l=0L
for i=0,Numl-1 do begin
for j=0,Nums-1 do begin
ogcp[0,l]=outx[i*Nums+j]
ogcp[1,l]=outy[i*Nums+j]
ogcp[2,l]=gcp[2,l]
ogcp[3,l]=gcp[3,l]
printf,lun,ogcp[0,l],ogcp[1,l],ogcp[2,l],ogcp[3,l]
l=l+1L
endfor
endfor
free_lun,lun
pixel_size = [926.62543314, 926.62543314]
out_name=path+prefix+arg1+'.0'+timestr+'_'+suffix+'.tif'
envi_file_query, fid_map, ns=ns, nl=nl, nb=nb
dims = [-1, 0, ns-1, 0, nl-1]
print,ns,nl
envi_doit, 'envi_register_doit', w_fid=fid_map, w_pos=0,w_dims=dims,method=6,$
out_name=out_name,pts=ogcp, proj=oproj,r_fid=r_fid1,pixel_size=pixel_size,X0=min(ogcp[0,*]),$
Y0=max(ogcp[1,*]),XSIZE=max(ogcp[0,*])-min(ogcp[0,*]),YSIZE=max(ogcp[1,*])-min(ogcp[1,*])
; XSIZE=ceil((max(ogcp[0,*])-min(ogcp[0,*]))/pixel_size[0]),$
; YSIZE=ceil((max(ogcp[1,*])-min(ogcp[1,*]))/pixel_size[1])
envi_file_mng, id=r_fid1, /remove
envi_file_mng, id=fid_map, /remove,/delete
close,/all
; envi_batch_exit
end
4、example_resize_doit.pro
pro example_resize_doit,day,path,prefix,suffix,fill_value
;
; First restore all the base save files.
;
envi, /restore_base_save_files
;
; Initialize ENVI and send all errors
; and warnings to the file batch.txt
;
envi_batch_init, log_file='batch.txt'
;
; Open the input file
file=file_search(path+prefix+day+'.*_'+suffix+'.tif')
in_fid=LON64ARR(size(file,/n_elements)+1-1)
dims=LON64ARR(5,size(file,/n_elements)+1-1)
pos=LON64ARR(1,size(file,/n_elements)+1-1)
use_see_through=intarr(size(file,/n_elements)+1-1)
see_through_val=intarr(size(file,/n_elements)+1-1)
for f=0,size(file,/n_elements)+1-2 do begin
envi_open_file, file[f], r_fid=fid
if (fid eq -1) then return
in_fid[f]=fid
envi_file_query, fid, ns=ns, nl=nl, nb=nb, DATA_TYPE=dtype
dims[*,f]=[-1,0, ns-1,0, nl-1]
pos[*,f]=0
use_see_through[f] = 1
see_through_val[f] = fill_value
endfor
;================================================
; get corner coordinates
;------------------------------------------------
ul_x=dblarr(size(file,/n_elements)+1-1)
ul_y=dblarr(size(file,/n_elements)+1-1)
lr_x=dblarr(size(file,/n_elements)+1-1)
lr_y=dblarr(size(file,/n_elements)+1-1)
for f=0,size(file,/n_elements)+1-2 do begin
ENVI_CONVERT_FILE_COORDINATES, in_fid[f], [0,dims[2,f]], $
[0,dims[4,f]], XMap, YMap,/to_map
ENVI_CONVERT_FILE_COORDINATES, in_fid[0], xf, yf, XMap, YMap
ul_x[f]=xf[0]
ul_y[f]=yf[0]
lr_x[f]=xf[1]
lr_y[f]=yf[1]
endfor
ulx=min(ul_x)
uly=min(ul_y)
x0=LON64ARR(size(file,/n_elements)+1-1)
y0=LON64ARR(size(file,/n_elements)+1-1)
for f=0,size(file,/n_elements)+1-2 do begin
x0[f]=LONG64(floor(ul_x[f]-ulx))
y0[f]=LONG64(floor(ul_y[f]-uly))
endfor
ulx=floor(ulx)
uly=floor(uly)
lrx=ceil(max(lr_x))
lry=ceil(max(lr_y))
ns_result=LONG64(lrx-ulx)+1
nl_result=LONG64(lry-uly)+1
ENVI_CONVERT_FILE_COORDINATES, in_fid[0], ulx,uly,$
XMap, YMap,/to_map
params = [6378245.0, 6356863.0,0.000000,105.000000,0.0, 0.0,25.0,47.0]
datum = 'WGS 84'
name = 'China'
units = ENVI_TRANSLATE_PROJECTION_UNITS('Meters')
pixel_size = [926.62543314,926.62543314]
mc = [0, 0, XMap, YMap]
proj = ENVI_GET_PROJECTION(FID = in_fid[0])
map_info = envi_map_info_create(mc=mc, ps=pixel_size,PROJ=proj)
xsize=ns_result*pixel_size[0]
ysize=nl_result*pixel_size[1]
outname=path+prefix+day+'.'+suffix+'.tif'
envi_doit, 'mosaic_doit', fid=in_fid, dims=dims, $
out_name=outname, background=fill_value,$
pixel_size=pixel_size, pos=pos, $
xsize=xsize, ysize=ysize,out_dt=dtype,$
x0=x0, y0=y0,georef=1,map_info=map_info,$
see_through_val=see_through_val,$
use_see_through=use_see_through
for f=0,size(file,/n_elements)+1-2 do begin
envi_file_mng, id=in_fid[f], /remove
endfor
;PTR_FREE,in_fid,dims,pos,use_see_through,see_through_val,ul_x,ul_y,lr_x,lr_y,x0,y0
y0=0
x0=0
lr_y=0
lr_x=0
ul_y=0
ul_x=0
; Exit ENVI
envi_open_file, 'G:\tangronglin\MOD11A1.2008144.LST_Day_1km.tif', r_fid=fid ; directory of mask,to be modified
;------------------------------------------------------
if (fid eq -1) then begin
envi_batch_exit
return
endif
;--------------------------------------------------------------------
envi_open_file, outname, r_fid=fid1; directory and file name to be modified
;---------------------------------------------------------------------
if (fid1 eq -1) then begin
envi_batch_exit
return
endif
envi_file_query, fid, ns=ns, nl=nl, nb=nb
envi_file_query, fid1, ns=ns1, nl=nl1, nb=nb1
ENVI_CONVERT_FILE_COORDINATES, fid, [0,ns-1],[0,nl-1], XMap, YMap,/to_map
ENVI_CONVERT_FILE_COORDINATES, fid1, xf, yf, XMap, YMap
dims = [-1, xf[0], xf[1], yf[0], yf[1]]
pos = lindgen(nb1)
; Perform the resize calculation.
; Make the output image twice as
; large in both X and Y. Use
; bilinear interpolation.
;--------------------------------------------------------
outname=path+'out\'+prefix+day+'.'+suffix+'.tif'
data1=float(ENVI_GET_DATA(fid=fid1,dims=[-1, 0, ns1-1, 0, nl1-1],pos=pos))
result=fltarr(ns,nl)
result=data1(xf[0]:xf[1],yf[0]:yf[1])
data=result
map_info=envi_get_map_info(fid=fid)
;print,map_info
ENVI_WRITE_ENVI_FILE,data,out_name=outname,$
map_info=map_info,/NO_OPEN
envi_file_mng, id=fid, /remove
envi_file_mng, id=fid1, /remove
;PTR_FREE,data1,result,data,dims,pos,map_info,outname
CLOSE,/all
end
已锁定
（作者置顶）