PRO jmap_combine, files, savefile, START_TIME=start_time, END_TIME=end_time, $ START_HT=start_ht, END_HT=end_ht, TYPE=type, REGIONS=regions, DEBUG=debug ;+ ; $Id: jmap_combine.pro,v 1.4 2009/07/29 16:16:45 cooper Exp $ ; NAME: jmap_combine ; ; PURPOSE: combine previously made jmaps into one file ; ; USAGE: jmap_combine, files, savefile ; ; KEYWORDS: ; START_TIME: desired start time of resultant map (TAI) ; END_TIME: desired end time of resultant map (TAI) ; START_HT: desired start height of resultant map (degrees) ; END_HT: desired end height of resultant map (degrees) ; TYPE: 1 for time-wise concatenation, 2 for height extension ; ; WRITTEN: Adam Herbst, 2007/07/26 ; ; $Log: jmap_combine.pro,v $ ; Revision 1.4 2009/07/29 16:16:45 cooper ; modified debug ; ; Revision 1.3 2009/07/20 18:38:41 nathan ; added /DEBUG ; ; Revision 1.2 2009/06/15 16:16:24 cooper ; modified widget, added satellite info ; ; Revision 1.1 2007/11/08 18:56:08 nathan ; moved from adam/jmaps2a ; ;- IF keyword_set(DEBUG) THEN debugon=1 ELSE debugon=0 IF(datatype(files) EQ 'UND') THEN BEGIN message, 'Usage: jmap_combine, files, savefile', /inf RETURN ENDIF nfiles = n_elements(files) hdr = create_jmap_strct2a() hdrs = replicate(hdr,nfiles) FOR i=0,nfiles-1 DO BEGIN junk = read_jmap2a(files[i],header=h, /NODATA) hdrs[i] = h ENDFOR time0 = min(hdrs.start_time) time1 = max(hdrs.end_time) IF(datatype(start_time) EQ 'UND') THEN start_time = time0 IF(datatype(end_time) EQ 'UND') THEN end_time = time1 IF(start_time LT 1.0) THEN $ ; passed as fraction of the total interval start_time = (time0 + start_time * (time1 - time0)) > time0 < time1 IF(end_time LE 2.0) THEN end_time = (time0 + end_time * (time1 - time0)) > time0 < time1 pix_loc = 0 all_ht = dblarr(1,2) all_time = dblarr(1,2) all_dim = lonarr(1,2) all_res = dblarr(1,2) newmap = 0 reg_num = 0 nreg = total(hdrs.nreg) IF debugon THEN BEGIN help,start_time,end_time,time0,time1,nfiles,nreg print,'START_TIME= ',utc2str(tai2utc(start_time),/ecs) print,' END_TIME= ',utc2str(tai2utc(end_time),/ecs) ENDIF reg_count = 0 FOR i=0,nfiles-1 DO BEGIN map = read_jmap2a(files[i], HEADER=h, PIXELS=pix) curr_pix = 0 IF(end_time LE h.start_time OR start_time GE h.end_time) THEN CONTINUE FOR j=0,h.nreg-1 DO BEGIN n_pix = h.reg_dim[j,0]*h.reg_dim[j,1] region = reform(pix[curr_pix:curr_pix+n_pix-1], h.reg_dim[j,0], h.reg_dim[j,1]) curr_time = h.start_time + (h.end_time-h.start_time) * float(h.border[j,0:1]) / h.nx curr_ht = h.min_r + (h.max_r-h.min_r) * float(h.border[j,2:3]) / h.ny IF(end_time LE curr_time[0] OR start_time GE curr_time[1]) THEN CONTINUE new_time = curr_time > start_time < end_time curr_scl = float(h.reg_dim[j,0]-1) / (curr_time[1]-curr_time[0]) curr_x1 = floor((new_time[0]-curr_time[0]) * curr_scl) > 0 curr_x2 = floor((new_time[1]-curr_time[0]) * curr_scl) < (h.reg_dim[j,0]-1) region = (temporary(region))[curr_x1:curr_x2, *] reg_num = [reg_num, (keyword_set(REGIONS) ? regions[i] : reg_count)] reg_count = reg_count + 1 newmap = [temporary(newmap), reform(region, n_elements(region))] all_time = [temporary(all_time), new_time] all_ht = [temporary(all_ht), curr_ht] all_dim = [temporary(all_dim), transpose([curr_x2-curr_x1+1, h.reg_dim[j,1]])] time_res = float(all_dim[hdr.nreg+1,0]) / (all_time[hdr.nreg+1,1] - all_time[hdr.nreg+1,0]) ht_res = float(all_dim[hdr.nreg+1,1]) / (all_ht[hdr.nreg+1,1] - all_ht[hdr.nreg+1,0]) all_res = [temporary(all_res), transpose([time_res,ht_res])] hdr.width[hdr.nreg] = h.width[j] hdr.nreg = hdr.nreg + 1 curr_pix = curr_pix + n_pix pix_loc = [pix_loc, pix_loc[n_elements(pix_loc)-1] + n_elements(region)] IF debugon THEN help, /str,h IF debugon THEN help, j, n_pix, curr_scl,curr_x1,curr_x2,reg_count,time_res,ht_res,curr_pix ENDFOR ENDFOR IF(hdr.nreg EQ 0) THEN BEGIN message, 'No data found in desired interval', /inf RETURN ENDIF newmap = (temporary(newmap))[1:*] all_time = all_time[1:*,*] all_ht = all_ht[1:*,*] all_dim = all_dim[1:*,*] all_res = all_res[1:*,*] reg_num = reg_num[1:*] hdr.min_r = min(all_ht) hdr.max_r = max(all_ht) hdr.start_time = min(all_time) hdr.end_time = max(all_time) get_utc, curr_utc, /ecs hdr.date_made = anytim2tai(curr_utc) ;rad_sort = reverse(sort(all_ht[*,1])) ;hdr.width[0:hdr.nreg-1] = (hdr.width[0:hdr.nreg-1])[rad_sort] map = newmap hdr.nx = ceil((hdr.end_time-hdr.start_time) * max(all_res[*,0])) hdr.ny = ceil((hdr.max_r-hdr.min_r) * max(all_res[*,1])) hdr.proj_type = 1 hdr.pa = h.pa regs = reg_num[uniq(reg_num)] regs = regs[sort(regs)] nreg = n_elements(regs) newhdr = hdr newhdr.width[*] = 0.0 count=0 FOR i=1,n_elements(hdrs)-1 DO BEGIN IF(hdrs[0].satellite EQ hdrs[i].satellite) THEN count+=1 ENDFOR IF(count EQ n_elements(hdrs)-1) THEN newhdr.satellite=hdrs[0].satellite & newhdr.hassat='yes~' reg_map = 0B reg_time = dblarr(nreg,2) reg_ht = dblarr(nreg,2) reg_dim = lonarr(nreg,2) reg_width = dblarr(nreg) reg_loc = lonarr(nreg+1) curr_pix = 0L FOR j=0,nreg-1 DO BEGIN inds = where(reg_num EQ regs[j]) reg_time[j,0] = min(all_time[inds,0]) reg_time[j,1] = max(all_time[inds,1]) reg_ht[j,0] = min(all_ht[inds,0]) reg_ht[j,1] = max(all_ht[inds,1]) reg_dim[j,0] = ceil((reg_time[j,1]-reg_time[j,0]) * max(all_res[inds,0])) reg_dim[j,1] = ceil((reg_ht[j,1]-reg_ht[j,0]) * max(all_res[inds,1])) reg_width[j] = max(hdr.width[inds]) region = bytarr(reg_dim[j,*]) FOR i=0,n_elements(inds)-1 DO BEGIN k = inds[i] section = reform(newmap[pix_loc[k]:pix_loc[k+1]-1], all_dim[k,0], all_dim[k,1]) x_ext = floor(reg_dim[j,0] * (all_time[k,*]-reg_time[j,0]) / (reg_time[j,1]-reg_time[j,0])) y_ext = floor(reg_dim[j,1] * (all_ht[k,*]-reg_ht[j,0]) / (reg_ht[j,1]-reg_ht[j,0])) region[x_ext[0]:x_ext[1]-1, y_ext[0]:y_ext[1]-1] = congrid(section, x_ext[1]-x_ext[0], y_ext[1]-y_ext[0]) if debugon then print,'k=',k,' x_ext=',x_ext,' y_ext=',y_ext ENDFOR curr_pix = curr_pix + n_elements(region) reg_loc[j+1] = curr_pix reg_map = [temporary(reg_map), reform(region, n_elements(region))] if debugon then help,curr_pix,reg_map if debugon then wait,5 ENDFOR reg_sort = reverse(sort(reg_ht[*,1])) newhdr.width[0:nreg-1] = reg_width[reg_sort] curr_pix = 0L map = bytarr(total(reg_dim[*,0]*reg_dim[*,1])) newhdr.nreg = nreg newhdr.version = 2 FOR j=0,nreg-1 DO BEGIN i = reg_sort[j] newhdr.reg_dim[j,*] = reg_dim[i,*] newhdr.border[j,0:1] = long(floor(newhdr.nx * (reg_time[i,*]-newhdr.start_time) / (newhdr.end_time-newhdr.start_time))) newhdr.border[j,2:3] = long(floor(newhdr.ny * (reg_ht[i,*]-newhdr.min_r) / (newhdr.max_r-newhdr.min_r))) n_pix = reg_loc[i+1]-reg_loc[i] map[curr_pix:curr_pix+n_pix-1] = reg_map[reg_loc[i]:reg_loc[i+1]-1] curr_pix = curr_pix + n_pix ENDFOR ;curr_pix = 0 ;FOR j=0,hdr.nreg-1 DO BEGIN ; i = rad_sort[j] ; hdr.reg_dim[j,*] = all_dim[i,*] ; region = newmap[pix_loc[i]:pix_loc[i+1]-1] ; map[curr_pix:curr_pix+n_elements(region)-1] = region ; hdr.border[j,0:1] = floor(hdr.nx * (all_time[i,*]-hdr.start_time) / (hdr.end_time-hdr.start_time)) ; hdr.border[j,2:3] = floor(hdr.ny * (all_ht[i,*]-hdr.min_r) / (hdr.max_r-hdr.min_r)) ; curr_pix = curr_pix + n_elements(region) ;ENDFOR ;stop write_jmap2a,map,savefile,newhdr ;hdr message, 'File saved as ' + savefile, /inf END