|
Hello,
I want to put another question regarding mosaicking of georefrenced images when they have different map projection i.e. they are in different UTM zone. My code is as under:
;Setup parameters for mosaic_doit
pro georef_mosaic_setup, fids=fids, dims=dims, out_ps=out_ps, $
xsize=xsize, ysize=ysize, x0=x0, y0=y0, map_info=map_infos
compile_opt strictarr, hidden
; some basic error checking
;
if keyword_set(dims) then $
if n_elements(fids) ne n_elements(dims[0,*]) then dims=0
;
if n_elements(fids) lt 2 then begin
xsize = -1
ysize = -1
x0 = -1
y0 = -1
return
endif
; if no DIMS passed in
;
nfiles = n_elements(fids)
if (keyword_set(dims) eq 0) then begin
dims = fltarr(5, nfiles)
for i=0, nfiles-1 do begin
envi_file_query, fids[i], ns=ns, nl=nl
dims[*,i] = [-1L, 0, ns-1, 0, nl-1]
endfor
endif
; - compute the size of the output mosaic (xsize and ysize)
; - store the map coords of the UL corner of each image since you'll need it later
;
UL_corners_X = dblarr(nfiles)
UL_corners_Y = dblarr(nfiles)
east = -1e34
west = 1e34
north = -1e34
south = 1e34
for i=0,nfiles-1 do begin
pts = [ [dims[1,i], dims[3,i]], $ ; UL
[dims[2,i], dims[3,i]], $ ; UR
[dims[1,i], dims[4,i]], $ ; LL
[dims[2,i], dims[4,i]] ] ; LR
envi_convert_file_coordinates, fids[i], pts[0,*], pts[1,*], xmap, ymap, /to_map
UL_corners_X[i] = xmap[0]
UL_corners_Y[i] = ymap[0]
east = east > max(xmap)
west = west max(ymap)
south = south 1])
map_info = envi_map_info_create(proj=proj, mc=[0,0,west,north], ps=out_ps)
temp = bytarr(10,10)
envi_enter_data, temp, map_info=map_info, /no_realize, r_fid=tmp_fid
PRINT, proj
; find the x and y offsets for the images
;
x0 = lonarr(nfiles)
y0 = lonarr(nfiles)
for i=0,nfiles-1 do begin
envi_convert_file_coordinates, tmp_fid, xpix, ypix, UL_corners_X[i], UL_corners_Y[i]
x0[i] = xpix
y0[i] = ypix
endfor
envi_file_mng, id=tmp_fid, /remove, /no_warning
end
;============
pro ex_georef_mosaic_full
envi, /restore_base_save_files
envi_batch_init, log_file='batch.txt'
;
;Open the input files
dir = 'G:/landsat_band5/2005/'
CD, dir
pr1='161028_028'
pr2='160028_028'
pr3='161029_029'
pr4='162029_029'
y='2005'
ex='_B50.TIF'
m ='09'
fname1 = FILE_SEARCH('*'+pr1+y+m+'*'+ex)
fname2 = FILE_SEARCH('*'+pr2+y+m+'*'+ex)
fname3 = FILE_SEARCH('*'+pr3+y+m+'*'+ex)
fname4 = FILE_SEARCH('*'+pr4+y+m+'*'+ex)
; fname = [fname1, fname2, fname3, fname4]
envi_open_data_file, fname1, r_fid=fid1
if (fid1 eq -1) then begin
return
endif
;
envi_open_data_file, fname2, r_fid=fid2
if (fid2 eq -1) then begin
return
endif
envi_open_data_file, fname3, r_fid=fid3
if (fid3 eq -1) then begin
return
endif
envi_open_data_file, fname4, r_fid=fid4
if (fid4 eq -1) then begin
return
endif
;
;Build the necessary keywords by querying the input images
;
envi_file_query, fid1, ns=file1_ns, nl=file1_nl, nb=file1_nb
envi_file_query, fid2, ns=file2_ns, nl=file2_nl, nb=file2_nb
envi_file_query, fid3, ns=file3_ns, nl=file3_nl, nb=file3_nb
envi_file_query, fid4, ns=file4_ns, nl=file4_nl, nb=file4_nb
;
;create fid array
fids = [fid1,fid2, fid3,fid4]
map_info1 = envi_get_map_info(fid=fid1)
map_info2 = envi_get_map_info(fid=fid2)
map_info3 = envi_get_map_info(fid=fid3)
map_info4 = envi_get_map_info(fid=fid4)
map_infos= [map_info1,map_info2,map_info3,map_info4]
print, map_infos
;
;create band arrays for POS keyword
pos1 = lindgen (file1_nb)
pos2 = lindgen (file2_nb)
pos3 = lindgen (file3_nb)
pos4 = lindgen (file4_nb)
pos = [[pos1],[pos2],[pos3],[pos4]]
;
;set output pixel size. This could also be queried from the input image.
out_ps = [30., 30.]
;
;call georef_mosaic_setup to calculate the dims, xsize, ysize, x0, y0,
;and map_info structure. Pass the FIDs of the two input files and output
;pixel size you want. If you want to use a spatial subset,
;pass in the dims already filled out
;
georef_mosaic_setup, fids=fids, out_ps=out_ps, dims=dims, xsize=xsize, ysize=ysize,$
x0=x0, y0=y0, map_info=map_infos
;
;setup the see through values and output bands names
use_see_through = [[1],[1],[1],[1]]
see_through_val = [[0],[0],[0],[0]]
dir = 'G:/landsat_band5/newMosac/'
CD, dir
out_name = y+'_'+m+'mosaic'
;
;Call the doit. Use a background value of 0 and set the output data type to byte.
;
envi_doit, 'mosaic_doit', fid=fids, pos=pos, dims=dims, out_name=out_name, $
r_fid=out_fid, xsize=xsize, ysize=ysize, x0=x0, y0=y0, georef=1, $
out_dt=1, pixel_size=out_ps, background= 255, see_through_val=see_through_val, $
use_see_through=use_see_through
;
;exit ENVI
;
envi_batch_exit
end
The program works well for images of same UTM zone but as the resultant mosaiced image takes one of the either zone then the different zone images dont fit together.
|