Function sp_quick_cube, data, $ scan_info, $ wmap=wmap, $ interp=interp, $ nbad=nbad, $ loud=loud, $ filldata=filldata, $ xytinterp=xytinterp, $ num_interps=num_interps, $ single_line_median_interp=single_line_median_interp, $ debug=debug ; ; INPUT PARAMETERS ; data : data arrays are pulled from the stks=sotsp_stks2struct(l1images) structure externally, see: sp_make_cube_index.pro ; eg, data == utimed = stks.utimed ; ; data == blapp = stks.blapp ; data == btapp = stks.btapp ; ... ; scan_info : from a call like: read_sotsp,sp2,index2,data,scan_info=scan_info,/quiet,/silent (where index2 is a level 2 index). See sp_make_cube_index.pro ; ; ; KEYWORD ARGUMENTS ; wmap : wmap is generated here, indexing where missing (and interpolated) columns are. ; interp : ; nbad : This is a vector of externally determined bad column numbers that are passed in and used to mark ; : (as opposed to auto-detect) a bad column (and a candidate for interpolation) in wmap ; loud : verbose output in some cases ; filldata : If filldata=0, then fill missing data with 0, ; otherwise fill with mean of all good data ; When calling, filldata=1 for conti, dopp2, dopp1 only. ; if /xytinterp, then interpolate using average delta between slit positions when calling, ; xytinterp = 1 should be called for slong, slat, utimed only ; xytinterp : if /xytinterp, then interpolate using average delta between slit positions ; num_interps : The maximum number of missing columns to fill with interpolated values. ; single_line_median_interp: ; debug : Causes slightly more verbose software execution. ; ; TSTART STRING '2014-02-20T23:33:11.005' ; TEND STRING '2014-02-21T01:31:58.923' ; IDL> vcube2 = sp_quick_cube(vint,scan_info,wmap=wmap,/interp,nbad=[304,1390,2463]) ;NOTE: I attempted to incorporate the comments from the following March 6, 2018 letter from Ted Tarbell, and am leaving ; the text in place here in the code for reference in case it is determined that this issue needs more work. ; ;Hi, ; ... ;However, I'm puzzled about how the interpolation is actually being done. Looking at (a small number of) movies, ;it looks like it is in effect copying the last good data over long stretches, more than just 5 columns. ;I thought for more than 5 columns (adjustable, num_interps) it would fill with 0 or mean of all the good data. ; ;Also, I thought I understood what it was doing down to about line 187, where wmap = -3 means temporal interpolation, ;and presumably all the spatial interpolated points after that have wmap = -2, and all filled points wmap = -1. ;Is that right? I printed numbers of -3, -2 and -1 cases and they don't really look like that. ; ;Down at the bottom (line 449), it seems to be ignoring the keyword, xytinterp, which was to be used for utimed, slat and slong; ;those are just the times/coordinates of the images and so have (very, very nearly) constant deltas between slit positions. ;Does it fall into that code for all the cubes? It should only do this for utimed, slat & slong and not go through all the ;previous interpolation code for those 3 variables. ; ;I have a new version of sp_make_cube_index too, but I'll hang on to it a little longer while you look at this. Then you can make gobs and gobs of these (including videos) and I'll look again. ; ;Thanks, ; ;Ted DEFAULT_MAX_NUM_INTERPS = 5 if ~keyword_set(num_interps) then num_interps = DEFAULT_MAX_NUM_INTERPS DEFAULT_FILLER_VALUE = 0.0 ;what datatype should this be? WMAP_FILLED_POINTS_FLAG = -1 SPATIAL_INTERPOLATION_FLAG = -2 TEMPORAL_INTERPOLATION_FLAG = -3 dsize = size(data) ny = dsize(2) ntot = dsize(1) slits = all_vals(scan_info.slitpos) nsummed = scan_info.scn_sum nx = scan_info.nslitpos/nsummed if (n_elements(slits) ne nx) then slits = scan_info.spmapctr-(scan_info.nslitpos/2)*scan_info.scn_step + $ indgen(nx)*scan_info.scn_step*nsummed hh = histogram(scan_info.slitpos) ;I am not sure of the role of hh. pgs nt = max(hh) cube = fltarr(nx,ny,nt) it=0 wlast = -1 wmap = -1+0*lonarr(nx,nt) ;nx is the number of columns in one image/data product, nt is the number of them ;This is the initial population of the data cube and the wmap. These will be updated via interpolations. ;stop, 'make sure data is non-zero. and step through the for loop' for i=0,ntot-1 do begin if (it le nt-1) then begin w = where(slits eq scan_info.slitpos(i)) if (w lt wlast+1) then it = it+1 if (it le nt-1) then begin cube(w,*,it) = data(i,*) ;print, 'cube being assigned a data value for i = ', i wmap(w,it) = i endif wlast = w endif endfor ;stop, 'check2a' if (it lt nt-1) then begin nt = it+1 cube = cube(*,*,0:it) wmap = wmap(*,0:it) endif ;stop, 'check2b' if (keyword_set(loud)) then print,'ntot,ny,nx,nt ',ntot,ny,nx,nt ;adjust wmap bad columns indices by adding externally passed in 'bad' column values. ..I guess it is assumed that whoever created this vector of bad column values did so from ;a set of columns where the duplicate columns had already been removed -- if not, we might have mismatching indices for intended-bad and marked-bad. Is there any way to check for ;or warn about this in a reasonable way? if keyword_set(nbad) then begin nb = n_elements(nbad) for i=0,nb-1 do begin w = where(wmap(*,*) eq nbad(i),count) if (count gt 0) then wmap(w) = WMAP_FILLED_POINTS_FLAG ;currently, this sets all endfor if (keyword_set(loud)) then print, 'nbad ', nbad endif ;This is the first of two interpolation sections in sp_quick_cube.pro. The interpolation in this section is spatial. (The section after this one performs temporal interpolation.) ;When invoked from sp_make_cube_index.pro the following SP data products are spatially interpolated: ;yes, spatially interpolated: blapp, btapp, conti, dopp1, dopp2, magincl, magazim, rufaz, ptot, bzcomp, fillfac (not all of these are written out by sp_make_cube_index.pro as output products). ;no, not spatially interpolated: magfield, slat, slong, utimed if keyword_set(interp) then begin if ((~keyword_set(xytinterp)) OR (nt eq 1)) then begin ;New version from Ted (slightly amended with parentheses for algorithmic clarity); what is the role of nt? num_pixels = n_elements(cube(0,*,0)) if (keyword_set(loud)) then print, 'NUM_PIXELS = ', num_pixels for i=0,nt-1 do begin w = where(wmap[*,i] eq WMAP_FILLED_POINTS_FLAG, count) ;where naively 'bad', and this should skip already been interpolated cols, because WMAP_FILLED_POINTS_FLAG == -1 values and they'll will be replaced by -2. ;w = where(wmap[*,i] eq WMAP_FILLED_POINTS_FLAG, count) ; where naively 'bad' -- no interpolations, yet. (should be identical to the prior line, which had "()" instead of '[]' -- remove.) if (keyword_set(loud)) then print,'Missing columns interpolated ',count ; ;note, the line below returns the indices of good values (as well as non -1 negative "bad" values), should it be "gt" rather than "ne" to make sure no -2 or -3 values are references?? w2 = where(wmap[*,i] ne WMAP_FILLED_POINTS_FLAG, ncount_w2a) ; replaced -1 by WMAP_FILLED_POINTS_FLAG to match construction above 20190425 ; ;filler is calculated here. Would it be better if the filler calculated below were the median, instead of the mean? --maybe try that. Right now the different ;data products end up looking different when they are joined as one multi-image product. if keyword_set(filldata) then begin filler=mean(cube[w2,*,i]) ;check this -- is this how to do it? Looking at this now, I think this only works because the -2 and -3 values haven't been set, yet. endif else begin filler=DEFAULT_FILLER_VALUE ;currently DEFAULT_FLLER_VALUE == 0 endelse ;then, filler is used to fill all 'bad' values. if count gt 0 then begin cube[w,*,i] = filler ;if any bad columns are found, set all of the 'naively bad' columns to the filler value. endif if (count ge 1) then begin ;count == the columns in the ith frame that is filled with WMAP_FILLED_POINTS_FLAG -- so these are bad, and this is the first time through the interpolation loop. for j=0,count-1 do begin ;count eq number of columns marked bad. jj = w(j) if keyword_set(debug) then print, 'in interpolation loop, processing column jj = '+ strtrim(jj,2) + ' for image i = '+ strtrim(i,2) interpolate_complete = 0; This is reset for each bad column, then when an interpolation scheme is applied it is set to 1/complete and we exit interpolation candidates for this column. ;EDGE CASES -- START ;is this the LEFT-most column (of the ith frame?) and the column has bad values? Can't interpolate, instead copy the if ((jj eq 0) && (interpolate_complete eq 0) && wmap(0,i) eq WMAP_FILLED_POINTS_FLAG) && ((size(wmap,/dimensions))[0] gt 1) then begin ;new 20200513: && ((size(wmap,/dimensions))[0] gt 1) if (wmap(1,i) ne WMAP_FILLED_POINTS_FLAG) then cube(0,*,i) = cube(1,*,i) else cube(0,*,i) = filler ;(is wmap(1,i) really the next column over that I am testing to make sure it is good, which is what I should be doing, or is this the next *good* column over, which maybe I shouldn't use. wmap(0,i) = SPATIAL_INTERPOLATION_FLAG ;mark the column as spatially interpolated. (I don't know if it really qualifies as spatial interpolation.) interpolate_complete = 1; if (keyword_set(loud)) then print, 'in edge case 1' endif ;Is this the RIGHT-most column (of the ith frame?) and the column has bad values? if ((jj eq nx-1) && (interpolate_complete eq 0) && (wmap(nx-1,i) eq WMAP_FILLED_POINTS_FLAG)) then begin if(wmap(nx-2,i) ne WMAP_FILLED_POINTS_FLAG) then cube(nx-1,*,i) = cube(nx-2,*,i) else cube(nx-1,*,i) = filler wmap(nx-1,i) = SPATIAL_INTERPOLATION_FLAG ;It looks like I might be marking both of these as spatially interpolated, even if I am assigning median filler values. interpolate_complete = 1 if (keyword_set(loud)) then print, 'in edge case 2' endif ;EDGE CASES -- END ; ;INTERNAL CASES -- START if ((jj gt 0) && (jj lt nx-1) && (interpolate_complete eq 0)) then begin ;not an edge case and interpolation is incomplete. ;SINGLE LINE INTERPOLATION (TWO CHOICES) -- START if ~keyword_set(single_line_median_interp) then single_line_median_interp = 0 ;so, the default approach is the box-pixel approach. if (((wmap(jj-1,i) ne WMAP_FILLED_POINTS_FLAG) && (wmap(jj-1,i) ne SPATIAL_INTERPOLATION_FLAG)) && $ ((wmap(jj+1,i) ne WMAP_FILLED_POINTS_FLAG) && (wmap(jj+1,i) ne SPATIAL_INTERPOLATION_FLAG))) then begin if (single_line_median_interp eq 1) then begin if (keyword_set(loud)) then print, 'Using OLD linear interpolation approach for one column only.' cube(jj,*,i) = 0.5*(cube(jj-1,*,i)+cube(jj+1,*,i)) endif else begin ;use the box-interpolation scheme. if (keyword_set(loud)) then print, 'Using NEW linear interpolation approach for one column only.' num_rows = (size(cube[0,*,0]))(2) ;is this robust for all cases, e.g., movie and non movie cases? TODO: check this out. ;for non edge interpolation ;maybe these (1/6.) weights should be different. I should probably test. cube(jj,1:num_rows-2,i) = (1/6.)*(cube(jj-1,1:num_rows-2,i)+cube(jj+1,1:num_rows-2,i)) ;straight across, top and bottom pixel cube(jj,1:num_rows-2,i) = cube(jj,1:num_rows-2,i)+ (1/6.)*(cube(jj-1,0:num_rows-3,i)+cube(jj+1,2:num_rows-1,i)) ;low to high across the pixel cube(jj,1:num_rows-2,i) = cube(jj,1:num_rows-2,i)+ (1/6.)*(cube(jj-1,2:num_rows-1,i)+cube(jj+1,0:num_rows-3,i)) ;high to low across the pixel ;still haven't filled the top and bottom pixel cube(jj,0,i) = (1/4.)*(cube(jj-1,0,i)+cube(jj+1,0,i)+cube(jj-1,1,i)+cube(jj+1,1,i)) ; bottom pixel, using only 4 values to interpolate cube(jj,num_rows-1,i) = (1/4.)*(cube(jj-1,num_rows-1,i)+cube(jj+1,num_rows-1,i)+cube(jj-1,num_rows-2,i)+cube(jj+1,num_rows-2,i)) ;top pixel, using only 4 values to interpolate endelse ; wmap(jj,i) = SPATIAL_INTERPOLATION_FLAG; this is a true spatial interpolation, unlike the edge cases above. interpolate_complete = 1; endif ;SINGLE LINE INTERPOLATION (TWO CHOICES) -- END ; a spline interpolation here would be fine to cover gaps of ; up to DEFAULT_MAX_NUM_INTERPS columns (or is over-ridden by num_interps), as long as it ; doesn't use any columns with wmap = -1 = WMAP_FILLED_POINTS_FLAG ;===================================================================== ;experiment: generalized number of columns: right hand spline -- START num_filler_cols_init=1 ; 2 ==> anything after nearest neighbor. ;===================================================================== ;this is what we are trying to pattern match, below for two cases, left and right.: ; (A.) if (0 eq interpolate_complete) && $ g-bbbb-gg ;complete flag not set ; (B.) (jj gt 0) && (jj lt nx-4) && $ ;within allowed range of jj values ; (C.) (wmap(jj-1,i) gt WMAP_FILLED_POINTS_FLAG) && (wmap(jj+4,i) gt WMAP_FILLED_POINTS_FLAG) && (wmap(jj+5) ne WMAP_FILLED_POINTS_FLAG) && $ ;valid values to calculate spline ; (D.) (wmap(jj,i) eq WMAP_FILLED_POINTS_FLAG) && (wmap(jj+1,i) eq WMAP_FILLED_POINTS_FLAG) && (wmap(jj+2,i) eq WMAP_FILLED_POINTS_FLAG) && (wmap(jj+3,i) eq WMAP_FILLED_POINTS_FLAG) num_filler_cols = num_filler_cols_init for_loop_range = num_filler_cols_init+num_interps-1 for num_filler_cols = num_filler_cols_init, for_loop_range do begin ;num_interps is currently set to 5 if (0 eq interpolate_complete) && ((jj+num_filler_cols) le nx-2) then begin ;Determine whether we should enter ;test condition (A for both LEFT and RIGHT cases. ;initialization of logic values -- "nothing up my sleeve, or left over from a previous calculation" contig_count_test = 0 contig_inner_bad_test = 0 boundary_vals_cols_valid = 0 LEFT_SPLINE_CONDITIONS_VALID = 0 ;initialize to failure RIGHT_SPLINE_CONDITIONS_VALID = 0 ;initialize to failure ;move these above and make a constant -- END ;This code sets a flag to determine if we can proceed with a right or left spline or neither. boundary_vals_cols_valid = 0 contig_count_test = 0 count_bad = 0 contig_inner_bad_test_temp = 0 ;initialize logical flags to failure ;TWO_AND_TWO_PIXEL_SPLINE_CONDITIONS_VALID boundary_vals_cols_valid = (jj gt 1) && $ (jj lt (nx - 3 - num_filler_cols)); ;test condition B. contig_count_test = ((wmap(jj-2,i) gt WMAP_FILLED_POINTS_FLAG) && (wmap(jj-1,i) gt WMAP_FILLED_POINTS_FLAG)) && $ ((wmap(jj+num_filler_cols,i) gt WMAP_FILLED_POINTS_FLAG) && (wmap(jj+num_filler_cols+1, i) gt WMAP_FILLED_POINTS_FLAG)) ;test condition C. contig_inner_bad_test_TEMP = where(wmap(jj:jj+num_filler_cols-1, i) eq WMAP_FILLED_POINTS_FLAG, count_bad) ;preperatory step (for test condition D, below). ; if (count_bad eq (num_filler_cols)) then contig_inner_bad_test = 1 else contig_inner_bad_test = 0 ;test condition D. ; TWO_AND_TWO_SPLINE_CONDITIONS_VALID = ((boundary_vals_cols_valid ne 0) && (contig_count_test ne 0) && (contig_inner_bad_test eq 1)) if (TWO_AND_TWO_SPLINE_CONDITIONS_VALID && (interpolate_complete eq 0)) then begin if (keyword_set(loud)) then box_message, ['In new two sided 2x2 Spline code','i = '+strtrim(string(i),2), $ 'jj-2 = '+strtrim(string(jj-2),2), 'wmap(jj-2,i) = '+strtrim(string(wmap(jj-2,i)),2),$ 'jj-1 = '+strtrim(string(jj-1),2), 'wmap(jj-1,i) = '+strtrim(string(wmap(jj-1,i)),2),$ 'jj = '+strtrim(string(jj),2) , 'wmap(jj,i) = '+strtrim(string(wmap(jj,i)),2) ,$ 'jj+1 = '+strtrim(string(jj+1),2), 'wmap(jj+1,i) = '+strtrim(string(wmap(jj+1,i)),2),$ 'jj+2 = '+strtrim(string(jj+2),2), 'wmap(jj+2,i) = '+strtrim(string(wmap(jj+2,i)),2) ] ;x_rspline definition and initialization x_rlspline = fltarr(4); These are the xvalues. The are consecutive on 'the left', skip two values (to be interpolated), and then have the last actual value. x_rlspline(0) = jj-2 x_rlspline(1) = jj-1 x_rlspline(2) = jj+num_filler_cols x_rlspline(3) = jj+num_filler_cols+1 ; t_rlspline = findgen(num_filler_cols + 2 + 2)+jj-2 sp_out_temp = fltarr( num_filler_cols+2+2)*0. y_rlspline = fltarr(2+2)*0.;these are the actual values in the neighboring cube entries. ; for index9=0, num_pixels-1 do begin ; for each pixel y_rlspline = fltarr(4)*0.;these are the actual values in the neighboring cube entries. y_rlspline(0) = cube(jj-2, index9,i); y_rlspline(1) = cube(jj-1, index9,i); y_rlspline(2) = cube(jj+num_filler_cols, index9,i); y_rlspline(3) = cube(jj+num_filler_cols+1,index9,i); ; if (keyword_set(loud)) then print, 'num_filler_cols = ', num_filler_cols if (keyword_set(loud)) then print, 'x_rlspline = ', x_rlspline if (keyword_set(loud)) then print, 'y_rlspline = ', y_rlspline ;t_rspline ; sp_out_temp = spline(x_rlspline, y_rlspline, t_rlspline, 100,/double) ;perform a spline calculation 'across' the bad, consecutive jj values if (keyword_set(loud)) then print, 'sp_out_temp = ', sp_out_temp if (keyword_set(loud)) then print, '' if (keyword_set(loud)) then print, 'Before correction: cube(jj-2:jj+num_filler_cols+1,index9,i) = ', cube(jj-2:jj+num_filler_cols+1,index9,i) cube(jj:jj+num_filler_cols-1,index9,i) = sp_out_temp[2:2+(num_filler_cols-1)] if (keyword_set(loud)) then print, 'After correction: cube(jj-2:jj+num_filler_cols+1,index9,i) = ',cube(jj-2:jj+num_filler_cols+1,index9,i) if (keyword_set(loud)) then print, '' wmap(jj:jj+num_filler_cols-1,i) = SPATIAL_INTERPOLATION_FLAG endfor ; jj = jj+num_filler_cols-1; no need to re-run code over already-splined values delvar, index9, sp_out_temp, x_rlspline, y_rlspline, t_rlspline interpolate_complete = 1 TWO_AND_TWO_SPLINE_CONDITIONS_VALID = 0; not needed because will be reassigned in code above in next loop-- but included for the sake of more robust code. endif ;"RIGHT_SPLINE" test for availability: g-(jj)(b*num_filler_cols-1)-gg boundary_vals_cols_valid = (jj gt 0) && $ (jj lt (nx - 2 - num_filler_cols)); ;test condition B. contig_count_test = (wmap(jj-1,i) gt WMAP_FILLED_POINTS_FLAG) && $ ((wmap(jj+num_filler_cols,i) gt WMAP_FILLED_POINTS_FLAG) && (wmap(jj+num_filler_cols+1, i) gt WMAP_FILLED_POINTS_FLAG)) ;test condition C. contig_inner_bad_test_TEMP = where(wmap(jj:jj+num_filler_cols-1, i) eq WMAP_FILLED_POINTS_FLAG, count_bad) ;preperatory step (for test condition D, below). ; if (count_bad eq (num_filler_cols)) then contig_inner_bad_test = 1 else contig_inner_bad_test = 0 ;test condition D. ; RIGHT_SPLINE_CONDITIONS_VALID = ((boundary_vals_cols_valid ne 0) && (contig_count_test ne 0) && (contig_inner_bad_test eq 1)) if (RIGHT_SPLINE_CONDITIONS_VALID && (interpolate_complete eq 0)) then begin if (keyword_set(loud)) then box_message, ['In new LARGE RANGE Right Spline code','i = '+strtrim(string(i),2), $ 'j-1 = '+strtrim(string(jj-1),2), 'wmap(jj-1,i) = '+strtrim(string(wmap(jj-1,i)),2),$ 'j = '+strtrim(string(jj),2) , 'wmap(jj,i) = '+strtrim(string(wmap(jj, i)),2) ,$ 'j+1 = '+strtrim(string(jj+1),2), 'wmap(jj+1,i) = '+strtrim(string(wmap(jj+1,i)),2),$ 'j+2 = '+strtrim(string(jj+2),2), 'wmap(jj+2,i) = '+strtrim(string(wmap(jj+2,i)),2) ] ;x_rspline definition and initialization x_rspline = fltarr(3); These are the xvalues. The are consecutive on 'the left', skip two values (to be interpolated), and then have the last actual value. x_rspline(0) = jj-1 x_rspline(1) = jj+num_filler_cols x_rspline(2) = jj+num_filler_cols+1 ; t_rspline = findgen(num_filler_cols+1+2) + (jj-1) sp_out_temp = fltarr(num_filler_cols+1+2)*0. y_rspline = fltarr(1+2)*0.;these are the actual values in the neighboring cube entries. for index10=0, (num_pixels-1) do begin ;for each pixel ;y_rspline definition and initialization y_rspline[0] = cube(jj-1, index10,i); y_rspline[1] = cube(jj+num_filler_cols ,index10,i); y_rspline[2] = cube(jj+num_filler_cols+1,index10,i); ; if (keyword_set(loud)) then print, 'num_filler_cols = ', num_filler_cols if (keyword_set(loud)) then print, 'x_rspline = ', x_rspline if (keyword_set(loud)) then print, 'y_rspline = ', y_rspline ; ;t_rspline sp_out_temp = spline(x_rspline, y_rspline, t_rspline, 100,/double) ;perform a spline calculation 'across' the bad, consecutive jj values if (keyword_set(loud)) then print, 'Before correction: cube(jj-1:jj+num_filler_cols+1,index10,i) = ', cube(jj-1:jj+num_filler_cols+1,index10,i) cube(jj:jj+num_filler_cols-1,index10,i) = sp_out_temp[1:1+(num_filler_cols)-1] if (keyword_set(loud)) then print, 'After correction: cube(jj-1:jj+num_filler_cols+1,index10,i) = ', cube(jj-1:jj+num_filler_cols+1,index10,i) if (keyword_set(loud)) then print, '' wmap(jj:jj+num_filler_cols-1,i) = SPATIAL_INTERPOLATION_FLAG endfor jj = jj+num_filler_cols-1; no need to re-run code over already-splined values ;MAKE SURE IT IS jj AND NOT j delvar, index10, sp_out_temp, x_rspline, y_rspline, t_rspline interpolate_complete = 1 RIGHT_SPLINE_CONDITIONS_VALID = 0; not needed because will be reassigned in code above in next loop-- but included for the sake of more robust code. endif ;END RIGHT SPLINE CODE ;"LEFT_SPLINE" test for availability: gg-(jj)(b*num_filler_cols-1)-g boundary_vals_cols_valid = (jj gt 1) && $ (jj lt (nx - 1 - num_filler_cols)) ; ;test condition B. contig_count_test = ((wmap(jj-2,i) gt -1) && (wmap(jj-1,i) gt WMAP_FILLED_POINTS_FLAG)) && $ (wmap(jj+num_filler_cols, i) gt WMAP_FILLED_POINTS_FLAG) ;test condition C. if (jj+num_filler_cols-1) ge (size(wmap))[1] then continue; return, -1 ;this filter may need to be put in more places ; will this one-d return work? contig_inner_bad_test_TEMP = where(wmap(jj:jj+num_filler_cols-1, i) eq WMAP_FILLED_POINTS_FLAG, count_bad) ;preperatory step (for test condition D, below). ; if (count_bad eq (num_filler_cols)) then contig_inner_bad_test = 1 else contig_inner_bad_test = 0 ;test condition D. LEFT_SPLINE_CONDITIONS_VALID = ((boundary_vals_cols_valid ne 0) && (contig_count_test ne 0) && (contig_inner_bad_test eq 1)) ; if (LEFT_SPLINE_CONDITIONS_VALID && (interpolate_complete eq 0)) then begin if (keyword_set(loud)) then box_message, ['In new LARGE RANGE LEFT Spline code'] if (keyword_set(loud)) then print, 'num_filler_cols = ', num_filler_cols ;x_lspline definition and initialization x_lspline = fltarr(2+1); These are the xvalues. The are consecutive on 'the left', skip two values (to be interpolated), and then have the last actual value. x_lspline(0) = jj-2 x_lspline(1) = jj-1 x_lspline(2) = jj+num_filler_cols ; t_lspline = findgen(num_filler_cols+2+1) + (jj-2) sp_out_temp = fltarr(num_filler_cols+2+1)*0. y_lspline = fltarr(2+1);these are the actual values in the neighboring cube entries. for index8=0, (num_pixels-1) do begin ;for each pixel y_lspline[0] = cube(jj-2, index8,i); y_lspline[1] = cube(jj-1, index8,i); y_lspline[2] = cube(jj+num_filler_cols,index8,i); if (keyword_set(loud)) then print, 'num_filler_cols = ', num_filler_cols if (keyword_set(loud)) then print, 'x_lspline = ', x_lspline if (keyword_set(loud)) then print, 'y_lspline = ', y_lspline ; sp_out_temp = spline(x_lspline, y_lspline, t_lspline, 100,/double) ;perform a spline calculation 'across' the bad, consecutive jj values if (keyword_set(loud)) then print, 'Before correction: cube(jj-2:jj+num_filler_cols,index8,i) = ', cube(jj-2:jj+num_filler_cols,index8,i) cube(jj:jj+num_filler_cols-1,index8,i) = sp_out_temp[2:2+num_filler_cols-1] ;WHY WOULDN't THESE JUST GO TO -1 == WMAP_FILLED_POINTS_FLAG if (keyword_set(loud)) then print, 'After correction: cube(jj-2:jj+num_filler_cols,index8,i) = ', cube(jj-2:jj+num_filler_cols,index8,i) if (keyword_set(loud)) then print, '' wmap(jj:jj+num_filler_cols-1,i) = SPATIAL_INTERPOLATION_FLAG endfor j = j+num_filler_cols-1; no need to re-run code over already-splined values delvar, index8, sp_out_temp, x_lspline, y_lspline, t_lspline interpolate_complete = 1 boundary_vals_cols_valid = 0 & contig_count_test = 0 & count_bad = 0 & contig_inner_bad_test_temp = 0 ;initialize logical flags to failure ;NEW LEFT_SPLINE_CONDITIONS_VALID = 0; not needed because will be reassigned in code above in next loop-- but included for the sake of more robust code. endif ;END LEFT SPLINE CODE endif;END OF NOT AN EDGE CASE and interpolation is incomplete. endfor;END of FOR LOOP through all possible j values. ;INTERNAL CASES -- END endif; END of SPATIAL INTERPOLATION CODE ;I need to verify that these edge cases are what ted wanted, if he wanted them. And, also, that the order of the temporal vs. spatial interpolations is correct. ;EDGE CASES -- START ;new test interpolation accross time images for start and stop ends START if (jj eq 0) && (i eq 0) && (wmap(jj,i) eq WMAP_FILLED_POINTS_FLAG) then begin if (keyword_set(loud)) then print, 'in edge case 5' if (keyword_set(loud)) then print, 'First pack missing in first data: copy next image data.' ; ;overwrite if certain special sub-cases: if ((i+1) le (nt-1)) && ((wmap(jj,i+1) ge 0)) then begin ;0 is the starting marker of good data cube(jj,*,i) = cube(jj,*,i+1) wmap(jj,i) = TEMPORAL_INTERPOLATION_FLAG; if (keyword_set(loud)) then print, 'First packet missing in first data: copy next image data.' interpolate_complete = 1 endif endif if (jj eq 0) && (i gt 0) && (i lt (nt-1)) && (wmap(jj,i) eq WMAP_FILLED_POINTS_FLAG) then begin if (keyword_set(loud)) then print, 'in edge case 6' if (keyword_set(loud)) then print, 'first packet missing, interior image: linear interpolation in time.' cube(jj,*,i) = 0.5*(cube(jj,*,i-1) + cube(jj,*,i+1)) wmap(jj,i) = TEMPORAL_INTERPOLATION_FLAG interpolate_complete = 1 endif if (jj eq (count-1)) && (i gt 0) && (i lt (nt-1)) && (wmap(jj,i) eq WMAP_FILLED_POINTS_FLAG) then begin if (keyword_set(loud)) then print, 'in edge case 7' if (keyword_set(loud)) then print, 'last pack, interior image: linear interpolation in time' cube(jj,*,i) = 0.5*(cube(jj,*,i-1) + cube(jj,*,i+1)) wmap(jj,i) = TEMPORAL_INTERPOLATION_FLAG interpolate_complete = 1 endif if (jj eq (count-1)) && (i gt 0) && (i eq (nt-1)) && (wmap(jj,i) eq WMAP_FILLED_POINTS_FLAG) then begin if (keyword_set(loud)) then print, 'in edge case 8' if (keyword_set(loud)) then print, 'last packet, last image: repeat the previous frame where the missing data is.' cube(jj,*,i) = cube(jj,*,i-1) wmap(jj,i) = TEMPORAL_INTERPOLATION_FLAG interpolate_complete = 1 endif ;EDGE CASES END endfor endif endfor endif; if ~keyword_set(xytinterp) endif ;if keyword_set(interp) ; ;; xytinterp was set so do this interpolation instead ;xytinterp, is to be used for utimed, slat and slong --see above note from Ted. At one point it wasn't being handled ;properly. Is that still the case? ;; Find the average delta between columns (at center) within the same image if keyword_set(xytinterp) then begin tlast = -1 ;should tlast be set to WMAP_FILLED_POINTS_FLAG here? ;original if n_elements(w) gt 1 then begin ;maybe it isn't supposed to be n_elementsf(w) gt than 1, but, rather, w gt 1? if (n_elements(w) gt 0) and (w[0] ge 0) then begin for jj = 0,nt-1 do begin count = 0; initialize in case we are going to skip because 1-d and therefore we cannot interpolate if (size(wmap))[1] gt 1 then begin; wrapper to catch 1-d case w = where(wmap(0:nx-2,jj) gt WMAP_FILLED_POINTS_FLAG and wmap(1:nx-1,jj) gt WMAP_FILLED_POINTS_FLAG, count) ;changed to gt to capture only non-interpolated values. I hope this was the intent. if (count gt 0) then begin if (keyword_set(loud)) then print, 'in edge case 9' delta = cube(w+1,ny/2,jj)-cube(w,ny/2,jj) if (tlast eq -1) then deltas = delta else deltas = [deltas,delta] ;should tlast be equal to WMAP_FILLED_POINTS_FLAG here? tlast = jj endif endif ;end of 1-d catch wrapper case endfor ;was "end" if (tlast gt -1) then delta = mean(deltas) else delta = 0.; ;Interpolate using delta from nearest good data for i=0,nt-1 do begin w = where(wmap(*,i) eq WMAP_FILLED_POINTS_FLAG, count) if (count ge 1) then begin wgood = where(wmap(*,i) ne WMAP_FILLED_POINTS_FLAG,countgood) for k=0,count-1 do begin ;THIS DOESN'T LOOK RIGHT TO ME -- the use of countgood, below. if (countgood eq 1) then begin cube(w(k),*,i) = cube(wgood(0),*,i)+delta*(w(k)-wgood(0)) endif else begin wdist = abs(w(k)-wgood) mindist = min(wdist,minw) cube(w(k),*,i) = cube(wgood(minw),*,i)+delta*(w(k)-wgood(minw)) endelse endfor endif endfor endif endif ;end of /xytinterp path -- it should only be used for utimed, slat and slong. In sp_make_cube_index.pro it is only and always called for utimed, slat and slong. return,cube end