;+ ; NAME: ; sphere_smooth ; PURPOSE: ; Smooths data on a sphere ; CATEGORY: ; sat/idl/toolbox/math ; CALLING SEQUENCE: FUNCTION sphere_smooth, FF , $ irange = irange , $ igrid = igrid , $ ; Input grid odim = odim_ , $ ; Output grid orange = orange , $ ogrid = ogrid , $ width = width , $ gauss = gauss , $ degrees = degrees , $ weightfnc = weightfnc , $ fillall = fillall , $ fillbad = fillbad , $ fillgood = fillgood , $ min_weight = min_weight ; INPUTS: ; FF array[N,M,L]; type: any numerical type ; N: longitudinal dimension ; M: latitudinal dimension ; L: number of arrays to be smoothed ; OPTIONAL INPUT PARAMETERS: ; irange=irange array[2,2]; type: float; default: [[0,360],[-90,90]] ; Specifies the range of the input grid on which ; F is defined ; irange[*,0]: longitude range ; irange[*,1]: latitude range ; The grid specifies a skymap in ; equirectangular projection ; igrid=igrid array[2,N,M]; type: float ; Instead of giving irange, the grid can be specified ; explicitly: ; igrid[0,*,*]: longitudes ; igrid[1,*,*]; latitudes ; Note that the grid must still be equirectangular ; ; odim=odim array[2]; type: integer; default: [N,M] ; dimensions of smoothed output skymap ; The dimensions must be smaller then/equal to ; the input dimensions and must be a factor of ; the input dimensions. ; orange=orange array[2,2]; type: float; default: irange ; Specifies the range of the output grid on which ; R is defined. ; orange[*,0]; longitude range ; orange[*,1]: latitude range ; Again, this is a equirectangular grid. ; ogrid=ogrid array[2,n,m]; type: float ; Instead of giving odim and orange, the output ; grid can be specified explicitly: ; grid[0,*,*]: longitudes ; grid[1,*,*]; latitudes ; width=width scalare; type: float; default: 1 degree ; width used to select points near each grid point ; to be used for a weighted mean. ; /gauss forces use of a Gaussian weighting function ; with width specified in keyword 'width'. ; /degrees if set, all angles are in degrees; default: radians ; ; weight_fnc=weight_fnc ; name of user-defined function to be used to calculate ; weight (see PROCEDURE). ; ; min_weight threshold for the weights. Smoothed values are ignored ; if the total weight of the contributing input bins ; is less then min_weight ; ; /fillbad only fill-in 'bad' values (NaN in input array FF) ; with smoothed values. This is only allowed if the ; input dimension (N,M) are the same as the output ; dimension odim. ; /fillgood only replace good values with smoothed values ; /fillall fill all values (both good and bad) by smoothed ; values, i.e., same as setting both /fillbad, /fillgood. ; This is also the default if neither /fillbad nor ; /fillgood is set. ; OUTPUTS: ; R array[n,m]; type: same as FF ; smoothed output skymap ; OPTIONAL OUTPUT PARAMETERS: ; igrid=igrid if not defined on input, they are returned as calculated ; ogrid=ogrid .. from the *dim and *range keywords. ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; InitVar, ToRadians, ToDegrees, sphere_distance, AngleRange ; ArrayLocation, gridgen ; PROCEDURE: ; > If neither /gauss, nor weight_fnc is set than all points within ; 'width' (angular separation) of a grid point will be averaged with ; the same weight. ; > If /gauss is set a Gaussian weight is used in the calculation ; of the average (with weight exp(-(elo/width)^2). Contributing points ; will extend to about 3*width away in angular separation. ; > A user-defined function weight_fnc can be used to calculate the weight. ; The function should have the form ; W = weight_fnc(elo, width=width, degrees=degrees) ; with 'elo' (in units indicated by keyword 'degrees') the elongation, ; and 'width' the input argument 'width' to this function. ; MODIFICATION HISTORY: ; JAN-2008, Paul Hick (UCSD/CASS) ; APR-2008, Paul Hick (UCSD/CASS) ; Bugfix. Weighting was done wrong. ; APR-2012, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ; Fixed problem with /fillbad keyword. Procedure will now abort ; if /fillbad is set and (some) of the points from the output grid ; are not present in the input grid (this is necessary to be able ; to retain good values from input array in the output array. ;- InitVar, fillbad , /key InitVar, fillgood, /key InitVar, fillall , /key fillall OR= 1-fillbad AND 1-fillgood InitVar, degrees, /key rpm = float(ToRadians(degrees=degrees)) dpm = float(ToDegrees(degrees=degrees)) InitVar, min_weight, 0.0 ; Pick up input dimensions. ; Output dimensions are same as input, if not specified idim = size(FF,/dim) InitVar, odim_, idim, set=odim more_than_one = n_elements(idim) EQ 3 CASE more_than_one OF 0: n_layer = 1 1: BEGIN n_layer = idim[2] IF n_elements(odim) EQ 2 THEN odim = [odim,n_layer] END ENDCASE ; If no function to calculate weights was defined explicitly, ; then do either a "box" smoothing (weight 1 inside, weight 0 ; outside the box), or us a Gaussian smoothing have_fnc = IsType(weightfnc,/defined) IF NOT have_fnc THEN BEGIn InitVar, width, 1.0/dpm InitVar, gauss, /key ENDIF ; The input grid (longitude, latitude) CASE IsType(igrid,/defined) OF 0: BEGIN InitVar, irange, [ [0.0,360.0], [-90.0,90.0] ]/dpm igrid = gridgen(idim[0:1], range=irange, /open, /multid) END 1: irange = [ [igrid[0,[0,idim[0]-1],0]],[igrid[1,0,[0,idim[1]-1]]] ] ENDCASE lng_igrid = reform(igrid[0,*,*]) ; Input longitudes as 2D array for one layer lat_igrid = reform(igrid[1,*,*]) ; Input latitudes as 2D array for one layer ilng = reform(igrid[0,*,0]) ; Input longitudes as 1D array ilat = reform(igrid[1,0,*]) ; Input latitudes as 1D array ; Grid resolutions (deg/bin) for input grid dilng = (igrid[0,idim[0]-1,0]-igrid[0,0,0])/(idim[0]-1) dilat = (igrid[1,0,idim[1]-1]-igrid[1,0,0])/(idim[1]-1) ; Index for latitude 0 in input array (not an integer in general) ilat_zero = -0.5-irange[0,1]/(irange[1,1]-irange[0,1])*idim[1] ; Number of skybins covering 360 degrees in ; longitude in input grid ilng_360 = round(360.0/dpm/dilng) ilng_full360 = idim[0] EQ ilng_360 IF ilng_full360 THEN message, /info, 'input map covers 360 deg in longitude' CASE IsType(ogrid,/defined) OF 0: BEGIN InitVar, orange, irange ogrid = gridgen(odim[0:1], range=orange, /open, /multid) END 1: orange = [ [ogrid[0,[0,odim[0]-1],0]],[ogrid[1,0,[0,odim[1]-1]]] ] ENDCASE olng = reform(ogrid[0,*,0]) ; Output longitudes as 1D array olat = reform(ogrid[1,0,*]) ; Output latitudes as 1D array ; Grid resolutions (deg/bin) for output grid dolng = (ogrid[0,odim[0]-1,0]-ogrid[0,0,0])/(odim[0]-1) dolat = (ogrid[1,0,odim[1]-1]-ogrid[1,0,0])/(odim[1]-1) ; Only allow rebinning to lower resolution ; (smaller array size) IF dolng LT dilng OR dolat LT dilat THEN BEGIN message, /info, 'resolution of results must be coarser than input resolution' RETURN, -1 ENDIF tmp = [dolng,dolat]/[dilng,dilat] nbin_factor = round(tmp) ; Only allow rebinning to dimensions that are a factor ; of the input dimension. IF abs(tmp[0]-nbin_factor[0]) GT 0.0001 OR abs(tmp[1]-nbin_factor[1]) GT 0.001 THEN BEGIN message, /info, 'dimension of result must be factor of the input dimension' RETURN, -1 ENDIF ; For Gaussian weights or if an explicit weight function ; is to be used, make sure to pick a big enough box. box_width = width*(1+2*max([gauss,have_fnc])) ; Set a box size on the equirectangular map that include ; all the bins that are potentially to be included. ; In the longitudinal direction the width scales as cos(lat)^-1 ; with a maximum of 180 degrees. lat_max = acos(box_width*dpm/180.0)/rpm lng_width = box_width/cos(((abs(olat)+0.5*box_width) < lat_max)*rpm) lat_width = box_width ; Find range in output grid that overlaps with the ; input grid and determine which part of the output ; skymap needs to be filled. ilng_index_near_olng = (olng-ilng[0])/dilng ; Array [onx] ilat_index_near_olat = (olat-ilat[0])/dilat ; Array [ony] ; fillbad can only be used if the output grid points exactly ; match point in the input grid so the good values can be ; copied from input to output array. IF fillbad THEN BEGIN IF (where(abs(ilng_index_near_olng-round(ilng_index_near_olng)) GT 0.0001))[0] NE -1 OR $ (where(abs(ilat_index_near_olat-round(ilat_index_near_olat)) GT 0.0001))[0] NE -1 THEN BEGIN message, /info, '/fillbad is set, but (some) points of output grid do not match input grid' RETURN, -1 ENDIF ENDIF ; Nearest points in input grid for points in output grid ilng_index_near_olng = round(ilng_index_near_olng) ; Array [onx] ilat_index_near_olat = round(ilat_index_near_olat) ; Array [ony] CASE ilng_full360 OF 0: BEGIN tmp = where( 0 LE ilng_index_near_olng AND ilng_index_near_olng LT idim[0] ) IF tmp[0] EQ -1 THEN BEGIN message, /info, 'output longitudes do not overlap input longitudes' RETURN, -1 ENDIF END 1: tmp = indgen(odim[0]) ENDCASE ; Min (>= 0) and max (< odim[0]) index for longitude in output map min_olng_index = min(tmp,max=max_olng_index) ; "Reference" index for longitude in center of output map ref_olng_index = (min_olng_index+max_olng_index)/2 ; Number of output longitudes that need to be filled n_olng = max_olng_index-min_olng_index+1 unit = replicate(1,n_olng) ostep = gridgen(n_olng, range=[min_olng_index,max_olng_index]-ref_olng_index) tmp = where( 0 LE ilat_index_near_olat AND ilat_index_near_olat LT idim[1]) IF tmp[0] EQ -1 THEN BEGIN message, /info, 'output latitudes do not overlap input latitudes' RETURN, -1 ENDIF ; Min (>= 0) and max (< odim[1]) index for latitude in output map min_olat_index = min(tmp,max=max_olat_index) ; Number of output latitudes that need to be filled n_olat = max_olat_index-min_olat_index+1 bad = BadValue(FF) ; Initialize the output array with bad values. R = make_array(dim=odim, value=bad) igrid = reform(igrid,2,idim[0]*idim[1],/overwrite); Input grid as [2,inx*iny] array ; In principle, the do-loop can be reduced further by using North/South symmetry. olat_done = bytarr(n_olat) epsilon = 1.0e-6 FOR j=min_olat_index,max_olat_index DO BEGIN ; Loop over output latitudes ; Check whether latitude has been done already IF olat_done[j] THEN continue ; Check whether opposite latitude is present. Since the weights are symmetric ; across the equator this can get done here too. ; (This potentially reduces the loop over output latitudes by ; a factor of two. There is some speed gain, but not much, ; and this doubles the memory requirements, so I'm not so sure ; this is a good idea; disable for now). j_other = -1 ;CASE abs(olat[j]) GT epsilon OF ;0: j_other = -1 ;1: j_other = (where(abs(olat+olat[j]) LT epsilon))[0] ;ENDCASE ; Mark latitude (and the opposite latitude, if present) as 'done' olat_done[j] = 1 IF j_other NE -1 THEN olat_done[j_other] = 1 ; Check the reference longitude at each output latitude: ; Reference point is P = [olng[ref_olng_index],olat[j]] ; lng_igrid and lat_igrid are [inx,iny] arrays. ipos = AngleRange(lng_igrid-olng[ref_olng_index],/pi,degrees=degrees) ipos = where(abs(ipos) LE lng_width[j] AND abs(lat_igrid-olat[j]) LE lat_width, n_ipos) ; Array [n_ipos] IF n_ipos EQ 0 THEN continue ; No points in box. Shouldn't happen ; ipos now are linear indices into the input map of points ; near the ref point in the center of the output map. ; Elongations for all n_ipos points in box. igrid is [2,inx*iny] array W = sphere_distance( igrid[*,ipos], [olng[ref_olng_index],olat[j]], degrees=degrees) W = float(W) ; Array [n_ipos] ; Weights for all n_ipos points are functions of elongation ; These are weights relative to the ref point in the center of the output map CASE 1 OF have_fnc: $ W = call_function(weight_fnc, W, width=width, degrees=degrees) gauss: BEGIN ; Gaussian weights W /= width W = exp(-W*W) END ELSE: $ ; "Box" weighting W LE= width ENDCASE ; This is primarily intended to reduce memory requirements ; for subsequent calculations: ditch points with very low weights tmp = where(W GE min_weight) IF tmp[0] NE -1 THEN BEGIN n_ipos = n_elements(tmp) ipos = ipos[tmp] W = W [tmp] ENDIF ; Array indices into input array FF for n_ipos points used to fill ; output map at the ref point P = [olng[ref_olng_index],olat[j]] ipos = ArrayLocation(ipos,dim=idim[0:1]) ; Array [2,n_ipos] into input map opos = [ref_olng_index,j] ; Array [2] into output map for ref point ; If the opposite latitude is to be included, set up an 2nd group. ; The olng indices remain the same; the olat index differences ; with olat[j_other] have the opposite sign as for olat[j]. n_opos = 1+(j_other NE -1) IF n_opos NE 1 THEN BEGIN ipos = SuperArray(ipos,n_opos,/trail) ; Array [2,n_ipos,n_opos] ipos[1,*,1] = round(2*ilat_zero-ipos[1,*,1]); Adjust lat index for opposite latitude W = SuperArray(W,n_opos,/trail) ; Array [n_ipos,n_opos] opos = SuperArray(opos,n_opos,/trail) ; Add output indices for opposite latitude opos[1,1] = j_other ; Array [2,n_opos] ENDIF ; Replicate the weights as n_olng identical rows. Each row represents ; one of the longitudes between min_olng_index and max_olng_index. W = SuperArray(W,n_olng,/trail) ; Array [n_ipos, n_opos, n_olng] ; There are n_olng longitudes in the output array at latitute olat[j]. ; We just calculated the points to be included in calculating the ; smoothed value at lng[ref_olng_index]. ; The other longitudes will include the same grouping of points, but ; shifted in longitude by a multiple of nbin_factor[0] input skybins. ; Set up arrays for n_olng groups of n_ipos skybins with longitudinal ; and latitudinal array indices into the input array FF. ; Longitudes: n_olng rows, increasing by nbin_factor[0] from row to row. ; Latitudes : n_olng identical rows. ipos = SuperArray(ipos,n_olng,/trail) ; Array [2,n_ipos,n_opos,n_olng] ipos = SubArray (ipos,element=0,add=replicate(1,n_ipos*n_opos)#(nbin_factor[0]*ostep)) opos = SuperArray(opos,n_olng,/trail) ; Array [2,n_opos,n_olng] opos = SubArray (opos,element=0,add=replicate(1,n_opos)#ostep) ; The longitudinal indices may be outside the valid range [0,idim[0]-1]. ; If irange[1]-irange[0] = 360.0 this can be fixed by ; adding/subtracting the equivalent of 360 to the array index. IF ilng_full360 THEN BEGIN tmp = SubArray(ipos,element=0) tmp mod= ilng_360 tmp += ilng_360 tmp mod= ilng_360 ipos = SubArray(ipos,element=0,replace=tmp) ENDIF x = SubArray(ipos,element=0) y = SubArray(ipos,element=1) tmp = where(0 LE x AND x LT idim[0] AND 0 LE y AND y LT idim[1]) IF tmp[0] EQ -1 THEN continue ; Shouldn't happen ; Convert to linear index into input array ipos = reform(ipos,[2,n_ipos*n_opos*n_olng],/overwrite) ipos = ArrayLocation(ipos[*,tmp], dim=idim, /onedim) ; Pick up function values from input array F = make_array(dim=[n_ipos,n_opos,n_olng,n_layer],value=bad) n_off = n_ipos*n_opos*n_olng m_off = idim[0]*idim[1] FOR i=0,n_layer-1 DO F[tmp+n_off*i] = FF[ipos+m_off*i] W = SuperArray(W,n_layer,/trail) ; Bad values are set to zero, and are given a zero ; weight. This effectively excludes them from ; contributing to the weighted mean. tmp = where(1-finite(F)) IF tmp[0] NE -1 THEN BEGIN F[tmp] = 0 W[tmp] = 0 ENDIF ; Calculate weighted mean for all points at olat[j] ; F,W = Array[n_ipos,n_opos,n_olng,n_layer] ; Sum over n_ipos dimension F = total(W*F,1)/total(W,1) ; Array [n_opos, n_olng, n_layer] ; Index array for lng and lat into the output array R x = SubArray(opos,element=0) ; Array [n_opos, n_olng] y = SubArray(opos,element=1) ; Array [n_opos, n_olng] ; Index arrays for lng and lat into the input array FF x_near = ilng_index_near_olng[x] ; Array [n_opos, n_olng] y_near = ilat_index_near_olat[y] ; Array [n_opos, n_olng] ; Index array for 3rd dimension (nr of layers) ; This is filled with zeroes. If the index of a layer is added ; it can be used to index 3D arrays R and FF z = lonarr(n_opos,n_olng) ; Array [n_opos, n_olng] F = reform(F,n_opos*n_olng,n_layer,/overwrite); Array [n_opos*n_olng, n_layer] FOR i=0,n_layer-1 DO BEGIN CASE 1 OF fillall : R[x,y,i+z] = F[*,i] ; Copy all smoothed values fillgood: BEGIN ; Overwrite only good values with smoothed values good = where( finite(FF[x_near,y_near,i]) ) IF good[0] NE -1 THEN R[x[good],y[good],i+z[good]] = F[good,i] END fillbad : BEGIN ; Overwrite only bad values with smoothed values ; and retain good values from the input array. good = where( finite(FF[x_near,y_near,i]), complement=tmp ) IF tmp [0] NE -1 THEN R[x[tmp ],y[tmp ],i+z[tmp ]] = F [tmp ,i] IF good[0] NE -1 THEN R[x[good],y[good],i+z[good]] = FF[x_near[good],y_near[good],i+z[good]] END ENDCASE ENDFOR ENDFOR ;IF fillbad OR fillall THEN BEGIN ; IF (where(1-finite(R)),n)[0] NE -1 THEN BEGIN ; message, /info, '/fillbad or /fillall set, but there still are NaNs left; fill with gridfill' ; FOR i=0,n_layer-1 DO R[*,*,i] = gridfill(R[*,*,i]) ; IF (where(1-finite(R)))[0] NE -1 THEN message, 'still NaNs left; should not have happened' ; message, /info, 'filled '+strcompress(n,/rem)+' bad values' ; ENDIF ;ENDIF RETURN, R & END