www.pudn.com > batch_reproject_resize_v03.rar > batch_reproject_resize_v03.pro


; batch_reproject_resize v 0.3 and batch_reproject_resize_reverse v 0.3 
 
; Author: Jonathan A. Greenberg (greenberg@ucdavis.edu) 
 
; batch_reproject_resize is designed to batch reproject and crop (or expand) a set 
;	of images to a single base image's projection and size.  This is 
;	helpful with preparing large sets of data for the purposes of change detection. 
 
; batch_reproject_resize_reverse is designed to reproject and crop (or expand) a single 
;	image based on a set of input base images.  This can be helpful if, for instance, 
;	you want a set of DEMs associated with a large number of flightline images, all pulled 
;	from the same source DEM. batch_reproject_resize_reverse requires batch_reproject_resize. 
 
; Dependencies: COMPARE_STUCT (included), 
;					from http://idlastro.gsfc.nasa.gov/ftp/pro/structure/compare_struct.pro 
 
; Known issues: I haven't done a lot of error checking yet. 
;	- if an image has an arbitrary projection it will throw an error. 
;	- hemispheres may cause problems (can someone check this?) 
;	- no output directory parameter yet, all output files go where the input files are. 
;	- possibly do a header check and set background value to the data ignore value. 
 
; Usage: 
;	Open in ENVI/IDL, go to run-> Compile batch_reproject_resize_vXX.pro or 
;		place this in your save_add directory.  At the ENVI command line: 
; 
;	batch_reproject_resize[,INFILES=string array,BASEFILE=string,INSUFFIX=string,OUTSUFFIX=string,$ 
;		WARP_METHOD={0 | 1 | 2 | 3},RESAMPLING={0 | 1 | 2},background=double,outbasename=string array 
; 
;	INFILES (optional): a string array containing the full path to each of the input files.  If 
;		no INFILES are given, a GUI window to select the files will be opened.  All infiles must 
;		have a proper projection set up in the header. 
;	BASEFILE (optional): a string containing the full path to the base file that all of the INFILES 
;		will be warped and cropped (or expanded) to.  If no BASEFILE is given, a GUI window to select 
;		the file will be opened. The basefile must have a proper projection set up in the header. 
;	INSUFFIX (optional): a string containing the input files' suffix (default is '.img').  The code 
;		uses this to set the output filename by replacing the INSUFFIX portion of the filename with OUTSUFFIX. 
;	OUTSUFFIX (optional): a string containing the output files' suffix (default is '_reproject.img').  The code 
;		uses this to set the output filename by replacing the INSUFFIX portion of the filename with OUTSUFFIX. 
;	WARP_METHOD (optional): keyword to set warp method as described in ENVI_CONVERT_FILE_MAP_PROJECTION.  We 
;		HIGHLY recommend leaving this set to its default (3, Rigorous (pixel-by-pixel)), as other settings 
;		can cause poor results. 
;	RESAMPLING (optional): keyword to set resampling method as described in ENVI_CONVERT_FILE_MAP_PROJECTION. 
;		Default is 0 (nearest neighbor).  Other settings can be 1: bilinear or 2: cubic convolution. 
;	BACKGROUND (optional): value to set the background values to during the reproject and any expansion 
;		of the input files.  The default is 0. 
;	OUTBASENAME (optional): don't use this, its there for compatibility with batch_reproject_resize_reverse. 
 
; batch_reproject_resize_reverse,[,INFILES=string,BASEFILE=string,INSUFFIX=string,OUTSUFFIX=string,$ 
;		WARP_METHOD={0 | 1 | 2 | 3},RESAMPLING={0 | 1 | 2},background=double 
;	INFILES (optional): a string array containing the full path to each of the input files.  If 
;		no INFILES are given, a GUI window to select the files will be opened.  All infiles must 
;		have a proper projection set up in the header. 
;	BASEFILE (optional): a string containing the full path to the base file that will be warped and 
;		cropped/expanded to each of the INFILES.  If no BASEFILE is given, a GUI window to select the 
;		 file will be opened.  The basefile must have a proper projection set up in the header. 
;	INSUFFIX (optional): a string containing the input files' suffix (default is '.img').  The code 
;		uses this to set the output filename by replacing the INSUFFIX portion of the filename with OUTSUFFIX. 
;	OUTSUFFIX (optional): a string containing the output files' suffix (default is '_reproject.img').  The code 
;		uses this to set the output filename by replacing the INSUFFIX portion of the filename with OUTSUFFIX. 
;	WARP_METHOD (optional): keyword to set warp method as described in ENVI_CONVERT_FILE_MAP_PROJECTION.  We 
;		HIGHLY recommend leaving this set to its default (3, Rigorous (pixel-by-pixel)), as other settings 
;		can cause poor results. 
;	RESAMPLING (optional): keyword to set resampling method as described in ENVI_CONVERT_FILE_MAP_PROJECTION. 
;		Default is 0 (nearest neighbor).  Other settings can be 1: bilinear or 2: cubic convolution. 
;	BACKGROUND (optional): value to set the background values to during the reproject and any expansion 
;		of the input files.  The default is 0. 
 
;   batch_reproject_resize and batch_reproject_resize_reverse is Copyright (C) 2008 Jonathan A. Greenberg 
 
;   This program is free software: you can redistribute it and/or modify 
;    it under the terms of the GNU General Public License as published by 
;    the Free Software Foundation, either version 3 of the License, or 
;    (at your option) any later version. 
 
;    This program is distributed in the hope that it will be useful, 
;    but WITHOUT ANY WARRANTY; without even the implied warranty of 
;    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
;    GNU General Public License for more details. 
 
;    You should have received a copy of the GNU General Public License 
;    along with this program.  If not, see . 
 
 
;+ 
; NAME: 
;       COMPARE_STRUCT 
; PURPOSE: 
;       Compare all matching tag names and return differences 
; 
; EXPLANATION: 
;       Compare all matching Tags names (except for "except_Tags") 
;       between two structure arrays (may have different struct.definitions), 
;       and return a structured List of fields found different. 
; 
; CALLING SEQUENCE: 
;       diff_List = compare_struct( struct_A, struct_B [ EXCEPT=, /BRIEF, 
;                                    /FULL, /NaN, /RECUR_A, /RECUR_B ) 
; INPUTS: 
;       struct_A, struct_B : the two structure arrays to compare. 
;       Struct_Name : for internal recursion use only. 
; OPTIONAL INPUT KEYWORDS: 
;               EXCEPT = string array of Tag names to ignore (NOT to compare). 
;               /BRIEF = number of differences found for each matching field 
;                                               of two structures is printed. 
;               /FULL = option to print even if zero differences found. 
;               /NaN = if set, then tag values are considered equal if they 
;                      are both set to NaN 
;               /RECUR_A = option to search for Tag names 
;                               in sub-structures of struct_A, 
;                               and then call compare_struct recursively 
;                               for those nested sub-structures. 
;               /RECUR_B = search for sub-structures of struct_B, 
;                               and then call compare_struct recursively 
;                               for those nested sub-structures. 
;       Note: 
;               compare_struct is automatically called recursively 
;               for those nested sub-structures in both struct_A and struct_B 
;               (otherwise cannot take difference) 
; OUTPUT: 
;       Returns a structure array describing differences found, 
;       which can be examined using print,diff_List or help,/st,diff_List. 
; PROCEDURE: 
;       Match Tag names and then use where function on tags. 
; EXAMPLE: 
;       Find the tags in the !X system variable which are changed after a 
;       simple plot. 
;       IDL> x = !X              ;Save original values 
;       IDL> plot, indgen(25)    ;Make a simple plot 
;       IDL> help,/str,compare_struct(x,!X)    ;See how structure has changed 
; 
;            and one will see that the tags  !X.crange and !X.S are changed 
;            by the plot. 
; MODIFICATION HISTORY: 
;       written 1990 Frank Varosi STX @ NASA/GSFC (using copy_struct) 
;       modif Aug.90 by F.V. to check and compare same # of elements only. 
;       Converted to IDL V5.0   W. Landsman   September 1997 
;       Added /NaN keyword W. Landsman  March 2004 
;- 
 
function compare_struct, struct_A, struct_B, EXCEPT=except_Tags, Struct_Name, $ 
                                        FULL=full, BRIEF=brief, NaN = NaN, $ 
                                        RECUR_A = recur_A, RECUR_B = recur_B 
 
   common compare_struct, defined 
   if N_params() LT 2 then begin 
       print,'Syntax - diff_List = compare_struct(struct_A, struct_B ' 
       print,'         [EXCEPT=, /BRIEF, /FULL, /NaN, /RECUR_A, /RECUR_B ]' 
       if N_elements(diff_List) GT 0 then return, diff_List else return, -1 
   endif 
 
        if N_elements( defined ) NE 1 then begin 
 
                diff_List = { DIFF_LIST, Tag_Num_A:0, Tag_Num_B:0, $ 
                                                Field:"",  Ndiff:0L } 
                defined = N_tags( diff_List ) 
          endif else diff_List = replicate( {DIFF_LIST}, 1 ) 
 
        Ntag_A = N_tags( struct_A ) 
        if (Ntag_A LE 0) then begin 
                message," 1st argument must be a structure variable",/CONTIN 
                return,diff_List 
           endif 
        Ntag_B = N_tags( struct_B ) 
        if (Ntag_B LE 0) then begin 
                message," 2nd argument must be a structure variable",/CONTIN 
                return,diff_List 
           endif 
 
        N_A = N_elements( struct_A ) 
        N_B = N_elements( struct_B ) 
 
        if (N_A LT N_B) then begin 
 
                message,"comparing "+strtrim(N_A,2)+" of first structure",/CON 
                message,"to first "+strtrim(N_A,2)+" of "+strtrim(N_B,2)+ $ 
                        " in second structure",/CONTIN 
 
                diff_List = compare_struct( struct_A, struct_B[0:N_A-1], $ 
                                                EXCEPT=except_Tags, $ 
                                                RECUR_A = recur_A, $ 
                                                RECUR_B = recur_B, $ 
                                                FULL=full, BRIEF=brief ) 
                return,diff_List 
 
          endif else if (N_A GT N_B) then begin 
 
                message,"comparing first "+strtrim(N_B,2)+" of "+ $ 
                        strtrim(N_A,2)+" in first structure",/CON 
                message,"to "+strtrim(N_B,2)+" in second structure",/CONTIN 
 
                diff_List = compare_struct( struct_A[0:N_B-1], struct_B, $ 
                                                EXCEPT=except_Tags, $ 
                                                RECUR_A = recur_A, $ 
                                                RECUR_B = recur_B, $ 
                                                FULL=full, BRIEF=brief ) 
                return,diff_List 
           endif 
 
        Tags_A = tag_names( struct_A ) 
        Tags_B = tag_names( struct_B ) 
        wB = indgen( N_elements( Tags_B ) ) 
        Nextag = N_elements( except_Tags ) 
 
        if (Nextag GT 0) then begin 
 
                except_Tags = [strupcase( except_Tags )] 
 
                for t=0,Nextag-1 do begin 
 
                        w = where( Tags_B NE except_Tags[t], Ntag_B ) 
                        Tags_B = Tags_B[w] 
                        wB = wB[w] 
                  endfor 
           endif 
 
        if N_elements( struct_name ) NE 1 then sname = "." $ 
                                          else sname = struct_name + "." 
 
        for t = 0, Ntag_B-1 do begin 
 
                wA = where( Tags_A EQ Tags_B[t] , nf ) 
 
                if (nf GT 0) then begin 
 
                     tA = wA[0] 
                     tB = wB[t] 
 
                     NtA = N_tags( struct_A.(tA) ) 
                     NtB = N_tags( struct_B.(tB) ) 
 
                     if (NtA GT 0 ) AND (NtB GT 0) then begin 
 
                        if keyword_set( full ) OR keyword_set( brief ) then $ 
                                                print, sname + Tags_A[tA], " :" 
 
                        diffs = compare_struct( struct_A.(tA), struct_B.(tB), $ 
                                                sname + Tags_A[tA], $ 
                                                EXCEPT=except_Tags, $ 
                                                FULL=full, BRIEF=brief ) 
                        diff_List = [ diff_List, diffs ] 
 
                      endif else if (NtA LE 0) AND (NtB LE 0) then begin 
 
                           if keyword_set(NaN) then begin 
                                  x1 = struct_b.(tB) 
                                  x2 = struct_a.(tA) 
                                  g = where( finite(x1) or finite(x2), Ndiff ) 
                                  if Ndiff GT 0 then $ 
                                    w = where( x1[g] NE x2[g], Ndiff ) 
                           endif else $ 
                            w = where( struct_B.(tB) NE struct_A.(tA) , Ndiff ) 
 
                                if (Ndiff GT 0) then begin 
                                        diff = replicate( {DIFF_LIST}, 1 ) 
                                        diff.Tag_Num_A = tA 
                                        diff.Tag_Num_B = tB 
                                        diff.Field = sname + Tags_A[tA] 
                                        diff.Ndiff = Ndiff 
                                        diff_List = [ diff_List, diff ] 
                                   endif 
 
                                if keyword_set( full ) OR $ 
                                  (keyword_set( brief ) AND (Ndiff GT 0)) then $ 
                                   print, Tags_A[tA], Ndiff, FORM="(15X,A15,I9)" 
 
                        endif else print, Tags_A[ta], " not compared" 
 
                 endif 
          endfor 
 
        if keyword_set( recur_A ) then begin 
 
                for tA = 0, Ntag_A-1 do begin 
 
                   if N_tags( struct_A.(tA) ) GT 0 then begin 
 
                        diffs = compare_struct( struct_A.(tA), struct_B, $ 
                                                sname + Tags_A[tA], $ 
                                                EXCEPT=except_Tags, $ 
                                                RECUR_A = recur_A, $ 
                                                RECUR_B = recur_B, $ 
                                                FULL=full, BRIEF=brief ) 
                        diff_List = [ diff_List, diffs ] 
                     endif 
                  endfor 
          endif 
 
        if keyword_set( recur_B ) then begin 
 
                for tB = 0, Ntag_B-1 do begin 
 
                   if N_tags( struct_B.(tB) ) GT 0 then begin 
 
                        diffs = compare_struct( struct_A, struct_B.(tB), $ 
                                                sname + Tags_B[tB], $ 
                                                EXCEPT=except_Tags, $ 
                                                RECUR_A = recur_A, $ 
                                                RECUR_B = recur_B, $ 
                                                FULL=full, BRIEF=brief ) 
                        diff_List = [ diff_List, diffs ] 
                     endif 
                  endfor 
          endif 
 
        w = where( [diff_List.Ndiff] GT 0, np ) 
        if (np LE 0) then w = [0] 
 
return, diff_List[w] 
end 
 
 
pro batch_reproject_resize,infiles=raster,basefile=basefilename,insuffix=insuffix,outsuffix=outsuffix,$ 
	warp_method=warp_method,resampling=resampling,background=background,outbasename=outbasename 
 
	if n_elements(raster) eq 0 then begin 
		rasterfilenames=envi_pickfile(title='Please select input raster(s):',/multiple_files) 
	    if rasterfilenames[0] eq '' then return 
	endif else rasterfilenames=raster 
	numfiles=n_elements(rasterfilenames) 
	rasterfids=lonarr(numfiles) 
	for i=0,numfiles-1 do begin 
		envi_open_file, rasterfilenames[i], r_fid=tempfid 
	   	rasterfids[i]=tempfid 
	endfor 
 
	if n_elements(basefilename) eq 0 then begin 
		basefilename=envi_pickfile(title='Please select base image to crop to...') 
	    if basefilename eq '' then return 
	endif 
 
	if n_elements(insuffix) eq 0 then insuffix='.img' 
	if n_elements(outsuffix) eq 0 then outsuffix='_reproject.img' 
 
	if insuffix eq outsuffix then begin 
		print,'insuffix must be different from outsuffix, exiting.' 
		return 
	end 
	tempsuffix='_temp.img' 
 
	if n_elements(warp_method) eq 0 then warp_method=3 
	if n_elements(resampling) eq 0 then resampling=0 
	if n_elements(background) eq 0 then background=0 
 
 
	rasterfilebasenames=strarr(numfiles) 
	rasterfiletempnames=strarr(numfiles) 
 
if n_elements(outbasename) lt numfiles then begin 
	for i=0,numfiles-1 do begin 
		rasterfilebasenames[i]=$ 
			file_dirname(rasterfilenames[i],/mark_directory)+file_basename(rasterfilenames[i],insuffix)+outsuffix 
		rasterfiletempnames[i]=$ 
			file_dirname(rasterfilenames[i],/mark_directory)+file_basename(rasterfilenames[i],insuffix)+tempsuffix 
	endfor 
endif else begin 
	for i=0,numfiles - 1 do begin 
	rasterfilebasenames[i]=$ 
			outbasename+outsuffix 
	rasterfiletempnames[i]=$ 
			outbasename+tempsuffix 
	endfor 
endelse 
 
	envi_open_file,basefilename,r_fid=basefid 
	envi_file_query,basefid,ns=basens,nl=basenl,nb=basenb,xstart=xstart,ystart=ystart,dims=basedims 
	o_proj = envi_get_projection(fid=basefid) 
	map_info = envi_get_map_info(fid=basefid) 
	o_pixel_size=map_info.ps[0:1] 
	inherit = envi_set_inheritance(basefid, basedims,[0],/full) 
 
	hemispherical_correction=[1,-1] 
	base_corners=fltarr(2,4) 
; Right now we assume N. hemisphere -- this should be changed 
	base_corners[*,0]=map_info.mc[2:3]-(map_info.mc[0:1]*o_pixel_size)*hemispherical_correction 
	base_corners[*,1]=base_corners[*,0]+(o_pixel_size*[basens,0])*hemispherical_correction 
	base_corners[*,2]=base_corners[*,0]+(o_pixel_size*[0,basenl])*hemispherical_correction 
	base_corners[*,3]=base_corners[*,0]+(o_pixel_size*[basens,basenl])*hemispherical_correction 
 
for i=0,numfiles-1 do begin 
; First, we reproject the entire image. 
	; If the image is already in the right projection, skip it... 
 
	tempfid=rasterfids[i] 
	tempoutname=rasterfiletempnames[i] 
	envi_file_query,tempfid,ns=ns,nl=nl,nb=nb,data_ignore_value=data_ignore_value,data_type=data_type 
	temp_proj=envi_get_projection(fid=tempfid) 
	if n_elements(compare_struct(temp_proj,o_proj)) gt 1 then begin 
	dims = [-1l, 0, ns-1, 0, nl-1] 
	pos=lindgen(nb) 
	if n_elements(data_ignore_value) eq 0 then data_ignore_value=background 
 
	envi_convert_file_map_projection, fid=tempfid, $ 
		pos=pos, dims=dims, o_proj=o_proj, $ 
	    o_pixel_size=o_pixel_size, $ 
	    out_name=tempoutname, warp_method=warp_method, $ 
	    resampling=resampling, background=data_ignore_value, $ 
	    r_fid=tempoutfid 
 
endif else begin 
	tempoutfid=tempfid 
endelse 
 
; Now we crop the image to the base image. 
	envi_file_query,tempoutfid,ns=outns,nl=outnl,nb=outnb 
	tempout_map_info=envi_get_map_info(fid=tempoutfid) 
	temp_o_proj = envi_get_projection(fid=tempoutfid) 
	tempout_corners=fltarr(2,4) 
 
; Right now we assume N. hemisphere -- this should be changed 
	tempout_corners[*,0]=tempout_map_info.mc[2:3]-(tempout_map_info.mc[0:1]*o_pixel_size)*hemispherical_correction 
	tempout_corners[*,1]=tempout_corners[*,0]+(o_pixel_size*[outns,0])*hemispherical_correction 
	tempout_corners[*,2]=tempout_corners[*,0]+(o_pixel_size*[0,outnl])*hemispherical_correction 
	tempout_corners[*,3]=tempout_corners[*,0]+(o_pixel_size*[outns,outnl])*hemispherical_correction 
 
	if tempout_corners[0,0] ge base_corners[0,0] then begin 
		tempout_xs = 0 
		tempout_xspos = long(abs(tempout_corners[0,0] - base_corners[0,0])/o_pixel_size[0])+1 
	endif else begin 
		tempout_xs = long(abs(tempout_corners[0,0] - base_corners[0,0])/o_pixel_size[0]) 
		tempout_xspos = 0 
	endelse 
 
	if tempout_corners[0,1] lt base_corners[0,1] then begin 
		tempout_xe = outns-1 
	endif else begin 
		tempout_xe = outns-long(abs(tempout_corners[0,1] - base_corners[0,1])/o_pixel_size[0])-2 
	endelse 
 
	tempout_xepos=tempout_xspos+(tempout_xe-tempout_xs) 
 
	if tempout_corners[1,0] lt base_corners[1,0] then begin 
		tempout_ys = 0 
		tempout_yspos = long(abs(tempout_corners[1,0] - base_corners[1,0])/o_pixel_size[1]) 
	endif else begin 
		tempout_ys = long(abs(tempout_corners[1,0] - base_corners[1,0])/o_pixel_size[1]) 
		tempout_yspos = 0 
	endelse 
 
	if tempout_corners[1,3] ge base_corners[1,3] then begin 
		tempout_ye = outnl-1 
		tempout_yepos = long(abs(tempout_corners[1,3] - base_corners[1,3])/o_pixel_size[1])-1 
	endif else begin 
		tempout_ye = outnl-long(abs(tempout_corners[1,3] - base_corners[1,3])/o_pixel_size[1])-1 
		tempout_yepos = 0 
	endelse 
 
	; Write out file 
	temppos=lindgen(outnb) 
	tempoutname=rasterfilebasenames[i] 
	openw,tempoutunit,tempoutname,/get_lun 
 
	tempoutarr=make_array([basens,outnb],type=data_type,value=background) 
 
;	fltarr(basens,nb)+float(background) 
 
	for j=0,tempout_yspos-1 do writeu,tempoutunit,tempoutarr 
 
	for j=tempout_ys,tempout_ye do begin 
		tempoutarr[tempout_xspos:tempout_xepos,*]=$ 
			envi_get_slice(fid=tempoutfid, line=j, pos=temppos, xs=tempout_xs, xe=tempout_xe) 
		writeu,tempoutunit,tempoutarr 
	endfor 
 
	tempoutarr=fltarr(basens,nb)+float(background) 
	for j=0,tempout_yepos-1 do writeu,tempoutunit,tempoutarr 
 
	close,tempoutunit 
	free_lun,tempoutunit 
 
	envi_setup_head,fname=tempoutname,nb=nb,/open,/write,ns=basens,nl=basenl,$ 
		data_type=data_type,interleave=1,inherit=inherit,data_ignore_value=data_ignore_value 
 
 
endfor 
 
end 
 
 
pro batch_reproject_resize_reverse,infiles=raster,basefile=basefile,insuffix=insuffix,outsuffix=outsuffix,$ 
			warp_method=warp_method,resampling=resampling,background=background 
 
	if n_elements(raster) eq 0 then begin 
		rasterfilenames=envi_pickfile(title='Please select base raster(s) to crop to...:',/multiple_files) 
	    if rasterfilenames[0] eq '' then return 
	endif else rasterfilenames=raster 
	numfiles=n_elements(rasterfilenames) 
	rasterfids=lonarr(numfiles) 
;	for i=0,numfiles-1 do begin 
;		envi_open_file, rasterfilenames[i], r_fid=tempfid 
;	   	rasterfids[i]=tempfid 
;	endfor 
 
	if n_elements(basefile) eq 0 then begin 
		basefilename=envi_pickfile(title='Please select input raster...') 
	    if basefilename eq '' then return 
	endif 
 
	if n_elements(outdirectory) eq 0 then begin 
		outdirectory=dialog_pickfile(title='Please select output directory...',/directory) 
	    if outdirectory eq '' then return 
	endif 
 
	if n_elements(insuffix) eq 0 then insuffix='.img' 
	if n_elements(outsuffix) eq 0 then outsuffix='_reproject.img' 
	tempsuffix='_temp.img' 
 
	if n_elements(warp_method) eq 0 then warp_method=3 
	if n_elements(resampling) eq 0 then resampling=0 
	if n_elements(background) eq 0 then background=0 
 
	for i=0,numfiles-1 do begin 
 
	tempbasename=outdirectory+file_basename(rasterfilenames[i],insuffix) 
 
	batch_reproject_resize,infiles=basefilename,basefile=rasterfilenames[i],insuffix=insuffix,outsuffix='_crop.img',$ 
			warp_method=warp_method,resampling=resampling,background=background,outbasename=tempbasename 
 
	endfor 
end