;+
; Name:
;
; DROT_RASTER
;
; Purpose:
;
; Differentially rotate solar raster maps
;
; Category:
;
; Maps
;
; Calling sequence:
;
; map_drot=DROT_RASTER(map [ , t_ref ] [ , REF_MAP=map_ref ]
; [ , ROLL=roll_angle , RCENTER=roll_center ]
; [ , /ADD_RDUR , /VERBOSE , SPHERE=0 , /KEEP ] )
;
; Inputs:
;
; map : input map (size of data array: [Nx,Ny])
; t_ref : reference time for rotation (def.: current time)
;
; Keyword parameters:
;
; REF_MAP = reference map for rotation and gridding. Its TIME, ROLL
; ROLL_CENTER values supersede those specified by t_ref,
; ROLL and RCENTER. Values of FULL_SIZE, SAME_GRID and MEMSAVE
; (see below) are also ignored.
; ROLL = final roll angle of the image
; RCENTER = center of rotation (if rolling)
; ADD_RDUR = add to the output map a tag (RDUR) containing an array that
; estimates the duration of the rotation interval on a
; pixel-by-pixel basis.
; VERBOSE = issue messages warning about off-limb pixels
; SPHERE = keyword passed to ARCMIN2HEL throuth ROT_XY (IN, default: 1;
; note that the default in ROT_XY is 0)
; KEEP = passed to DROT_COORD: keeps same same P,B0,R values when
; rotating: use TIME tag as reference epoch for computing these
; values
; FULL_SIZE = expand image size to fit the entire rotated FOV.
; SAME_GRID = retains old grid (ignored when dealing with an array of maps)
; MEMSAVE = for arrays of maps, do not store the rotated coordinates
; in a sigle array, but rather regrid images based on estimates
; of the new positions of the FOV corners.
; all keywords accepted by INTERP2D, except /REGULAR
;
; Outputs:
;
; map_drot : new map (differentially rotated)
;
; Common blocks:
;
; None
;
; Calls:
;
; VALID_MAP,GET_UTC,ANYTIM2TAI,UNPACK_MAP,SOHO_VIEW,USE_SOHO_VIEW,
; USE_EARTH_VIEW,GET_MAP_XRANGE,GET_MAP_YRANGE,GET_ARR_CENTER,ROT_XY,
; REP_TAG_VALUE,ADD_PROP,ROLL_XY,DROT_COORD,PB0R
;
; Description:
;
; A map is considered to represent a raster if it contains a tag,
; named START, giving the start times of each exposure (column) of the
; raster.
; The coordinates of each column (exposure) are differentially rotated
; to the reference time (t_ref). The data array is then interpolated on
; a regular grid. If an array of maps is given, the new grid is the
; same for all the maps, and covers the largest rectangular FOV included
; in (or including, if FULL_SIZE is set) the rotated FOV's of all maps.
; For a single map, the original grid can be also used (keyword SAME_GRID,
; useful for small time intervals).
; When ROLL and/or RCENTER, are specified, the resulting maps are
; also rotated to the new roll angle.
; When a reference map is provided, its TIME, ROLL, ROLL_CENTER, and
; grid parameters override all other input values.
;
; Side effects:
;
; None
;
; Restrictions:
;
; - Assumes that rasters are built in the E-W direction or viceversa,
; not N-S (see DROT_COORD).
; - Mixing maps (including MAP_REF) with both SOHO and EARTH coordinates
; is not allowed.
;
; Notes:
;
; The algorithm employed here (for zero roll angles) is the following:
;
; 1) Differentially rotate coordinates (-> ROT_XY)
; 2) Interpolate data on the new (irregular) grid (-> INTERP2D,REGULAR=0)
;
; On the other hand, the algorithm employed by DROT_MAP (current version
; as of January 2001) is the following:
;
; 1) Differentially rotate coordinates (-> ROT_XY)
; 2) Regrid on a regular grid
; 3) Rotate back to initial time (-> ROT_XY)
; 4) Interpolate over old grid (-> INTERP2D,REGULAR=1)
;
; In the latter case, the advantage is the ability of interpolating over
; a regular grid (faster algorithm); the disadvantage is that two
; differential rotations are required. Moreover, here Pass 3 cannot be
; coded in one step, since for rasters there is no unique reference time
; to rotate back to, and, more important, Pass 2 mixes together pixels
; taken at different times.
;
; Modification history:
;
; V. Andretta, 26/Jan/1999 - Written
; V. Andretta, 28/Feb/1999 - Added ROLL, RCENTER, REF_MAP, SPHERE
; keywords
; V. Andretta, 11/Apr/1999 - Added KEEP keyword
; V. Andretta, 20/Apr/1999 - Specified ambiguous keyword "R" in ROT_XY;
; fixed erroneous definition of ROLL/ROLL_CENTER kwds for arrays of maps
; V. Andretta, 19/Jan/2001 - Informs INTERP2D that after rotation, grid
; is no longer uniformly spaced (REGULAR=0 was the default in an
; earlier version of that routine).
; - Fixed a problem with center of grid in the presence of a roll angle.
; - Added SAME_GRID, FULL_SIZE, and MEMSAVE keywords.
; - Set KEEP to 1 by default.
; - Attempts to take round-off errors into account when computing
; grid parameters.
; V. Andretta, 5/Dec/2001 - Taken into account changes in UNPACK_MAP
; regarding tag and keyword names specifying roll properties
; (ROLL -> ROLL_ANGLE, RCENTER -> ROLL_CENTER).
;
; Contact:
;
; andretta@na.astro.it
;-
;===============================================================================
;
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;
;
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function DROT_RASTER,map,t_ref,REF_MAP=map_ref $
,ROLL=roll_angle,RCENTER=roll_center,ADD_RDUR=add_rdur $
,VERBOSE=verbose,_EXTRA=_extra,SPHERE=sphere,KEEP=keep $
,FULL_SIZE=full_size,SAME_GRID=same_grid,MEMSAVE=memsave
ON_ERROR,2
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;%%% Check input and some definitions
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
map_drot=0
if N_PARAMS() lt 1 then begin
PRINT,'%E> DROT_RASTER: Usage: map_drot=DROT_RASTER,map,t_ref'
PRINT,' or: map_drot=DROT_RASTER,map,REF_MAP=ref_map'
RETURN,map_drot
endif
if VALID_MAP(map) eq 0 then begin
PRINT,'%E> DROT_RASTER: Input structure not a map'
RETURN,map_drot
endif
Nmaps=N_ELEMENTS(map)
verb=KEYWORD_SET(verbose)
if N_ELEMENTS(sphere) eq 0 then sphere=1
do_add_rdur=KEYWORD_SET(add_rdur)
using_soho_view=SOHO_VIEW()
;% Set keyword KEEP to 1 by default
if N_ELEMENTS(keep) eq 0 then keep=1
keep_angles=KEYWORD_SET(keep)
;% Get reference time and roll properties from input parameters/keywords
;% Get roll values for all the maps
map_roll=FLTARR(Nmaps)
map_rcenter=FLTARR(2,Nmaps)
for i=0,Nmaps-1 do begin
UNPACK_MAP,map(i),ROLL_ANGLE=curr_roll,ROLL_CENTER=curr_rcenter
curr_roll=FLOAT(curr_roll) mod 360
map_roll(i)=curr_roll
map_rcenter(*,i)=curr_rcenter
endfor
;% Initial estimate of reference time
GET_UTC,t0 & t0=ANYTIM2TAI(t0)
if not VALID_MAP(map_ref) then begin
;% Estimate (or use input) reference time for rotation...
if N_ELEMENTS(t_ref) ne 0 then begin
errmsg=''
t0=ANYTIM2TAI(t_ref(0),ERRMSG=errmsg)
if errmsg ne '' then PRINT,'%W> DROT_RASTER: '+errmsg
endif
;% ...and final roll angle and center of rotated map(s)
if N_ELEMENTS(roll_angle) ne 0 then $
roll=REPLICATE(FLOAT(roll_angle(0)) mod 360,Nmaps) $
else $
roll=map_roll
if N_ELEMENTS(roll_center) eq 2 then $
rcenter=TRANSPOSE([[REPLICATE(FLOAT(roll_center(0)),Nmaps)] $
,[REPLICATE(FLOAT(roll_center(1)),Nmaps)]]) $
else $
rcenter=map_rcenter
endif
;% Check if map represents a raster
sz=SIZE(map(0).data)
Nx=sz(1)
Ny=sz(2)
multi_time=0
if TAG_EXIST(map(0),'START') then begin
if N_ELEMENTS(map(0).start) ne Nx then $
PRINT,'%W> DROT_RASTER: ' $
+'Number of START times inconsistent with DATA 1st dimension' $
else $
multi_time=1
endif
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;%%% Create new grid
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;%%%%%
if VALID_MAP(map_ref) then begin
;% Get grid parameters from reference map; also get reference time and
;% roll parameters
UNPACK_MAP,map_ref(0),t0,TAG='TIME' $
,NX=Nxg,NY=Nyg,XC=xcg0,YC=ycg0,DX=dxg,DY=dyg $
,ROLL_ANGLE=roll,ROLL_CENTER=rcenter
roll=FLOAT(roll) mod 360
xpg0=GET_MAP_XP(map_ref(0))
ypg0=GET_MAP_YP(map_ref(0))
;% Final roll angle for all maps
roll=REPLICATE(roll,Nmaps)
rcenter=TRANSPOSE([[REPLICATE(FLOAT(rcenter(0)),Nmaps)] $
,[REPLICATE(FLOAT(rcenter(1)),Nmaps)]])
;%%%%%
endif else begin
if Nmaps gt 1 or KEYWORD_SET(same_grid) eq 0 then begin
;% Define a common grid for all the maps to be rotated
if KEYWORD_SET(memsave) then begin
;% Define new field of view
;%%% Define current field of view
xr=FLTARR(2,2,Nmaps) ; => first coordinate: E-W, second coordinate: N-S
yr=FLTARR(2,2,Nmaps)
tr=DBLARR(2,Nmaps)
;% Times of the first and last column
if multi_time eq 1 then begin
tr(0,*)=map.start(0)
tr(1,*)=map.start(Nx-1)
endif else begin
tr(0,*)=ANYTIM2TAI(map.time)
tr(1,*)=ANYTIM2TAI(map.time)
endelse
;% Times and coordinates of the FOV corners
for i=0,Nmaps-1 do begin
xr(*,0,i)=GET_MAP_XRANGE(map(i))
xr(*,1,i)=xr(*,0,i)
yr(0,*,i)=GET_MAP_YRANGE(map(i))
yr(1,*,i)=yr(0,*,i)
endfor
;%%% Differentially rotate current field of view
xr_drot=xr
yr_drot=yr
view_type=GET_MAP_PROP(map(0),/SOHO)>0
if view_type then USE_SOHO_VIEW else USE_EARTH_VIEW
for i=0,Nmaps-1 do begin
;%%% Roll back boundaries of field of view
curr_roll=map_roll(i)
curr_rcenter=map_rcenter(*,i)
xr_i=xr(*,*,i)
yr_i=yr(*,*,i)
if curr_roll ne 0 then ROLL_XY,xr(*,*,i),yr(*,*,i) $
,-curr_roll,CENTER=curr_rcenter,xr_i,yr_i
xr_drot_i=xr_i
yr_drot_i=yr_i
;%%% If KEEP keyword is set, compute here P,B0,R values; these
;%%% values will be passed to ROT_XY
if keep_angles then begin
angles=PB0R(map(i).time,SOHO=SOHO_VIEW())
P_sun=angles(0)
B0_sun=angles(1)
;% Convert solar radius to observer's distance, in units of
;% solar radii
D_sun=1./ATAN(angles(2)/(60.*!RADEG))
endif
;%%% Differentially rotate field of view
for iix=0,1 do begin
crd=ROT_XY(xr_i(iix,*),yr_i(iix,*),TSTART=tr(iix,i),TEND=t0 $
,OFFLIMB=offlimb,INDEX=disk_index $
,BACK_INDEX=back_index $
,P=P_sun,B0=B0_sun,R0=D_sun,KEEP=keep_angles $
,SPHERE=KEYWORD_SET(sphere))
;% Points still on disk:
if disk_index(0) ge 0 then begin
xr_drot_i(iix,disk_index)=crd(disk_index,0)
yr_drot_i(iix,disk_index)=crd(disk_index,1)
endif
;% Points which have rolled over behind the limb: leave them where
;% they crossed the limb
if back_index(0) ge 0 then begin
;% Compute latitude of points behind the visible disk
lat=ARCMIN2HEL(xr_i(iix,back_index)/60,yr_i(iix,back_index)/60 $
,DATE=tr(iix,i),SOHO=SOHO_VIEW() $
,P=P_sun,B0=B0_sun,R0=D_sun $
,SPHERE=KEYWORD_SET(sphere))
lat=lat(0,*)
;% If not done already, compute solar B0
if keep_angles eq 0 then begin
angles=PB0R(tr(iix,i),SOHO=SOHO_VIEW())
B0_curr=angles(1)
endif else $
B0_curr=B0_sun
;% Longitude of the limb at each latitude (would be 180 degrees
;% if B0=0)
lon_limb=-TAN(B0_curr/180.*!DPI)*TAN(lat/180.*!DPI)
good_val=WHERE(ABS(lon_limb) le 1,Ngood_val)
if Ngood_val gt 0 then begin
lon_limb(good_val)=ACOS(lon_limb(good_val))/!DPI*180.
;% Transform new coordinates back to cartesian
arc=HEL2ARCMIN(lat(good_val),lon_limb(good_val) $
,DATE=tr(iix,i),SOHO=SOHO_VIEW() $
,P=P_sun,B0=B0_sun,R0=D_sun)*60
xr_drot_i(iix,back_index(good_val))=arc(0,*)
yr_drot_i(iix,back_index(good_val))=arc(1,*)
endif
endif
endfor
;%%% Return FOV coordinates to original roll angle
if curr_roll ne 0 then ROLL_XY,xr_drot_i,yr_drot_i $
,+curr_roll,CENTER=curr_rcenter,xr_drot_i,yr_drot_i
xr_drot(*,*,i)=xr_drot_i
yr_drot(*,*,i)=yr_drot_i
endfor
if using_soho_view then USE_SOHO_VIEW else USE_EARTH_VIEW
;%%% New field of view
if KEYWORD_SET(full_size) then begin
x1=MIN(xr_drot(0,*,*))
x2=MAX(xr_drot(1,*,*))
y1=MIN(yr_drot(*,0,*))
y2=MAX(yr_drot(*,1,*))
endif else begin
x1=MAX(xr_drot(0,*,*))
x2=MIN(xr_drot(1,*,*))
y1=MAX(yr_drot(*,0,*))
y2=MIN(yr_drot(*,1,*))
endelse
endif else begin
;%%% Define array of rotated coordinates
xp_drot=MAKE_ARRAY(Nx,Ny,Nmaps,TYPE=SIZE(GET_MAP_XP(map(0)),/TYPE))
yp_drot=xp_drot
;%%% Rotate coordinates for all the rasters in the array
for i=0,Nmaps-1 do begin
xyr=DROT_COORD(map(i),t0,VERBOSE=verb,SPHERE=KEYWORD_SET(sphere) $
,KEEP=keep_angles)
xp_drot(*,*,i)=xyr(*,*,0)
yp_drot(*,*,i)=xyr(*,*,1)
xyr=0
endfor
;%%% New field of view
if KEYWORD_SET(full_size) then begin
x1=MIN(xp_drot( 0 , * ,*))
x2=MAX(xp_drot(Nx-1, * ,*))
y1=MIN(yp_drot( * , 0 ,*))
y2=MAX(yp_drot( * ,Ny-1,*))
endif else begin
x1=MAX(xp_drot( 0 , * ,*))
x2=MIN(xp_drot(Nx-1, * ,*))
y1=MAX(yp_drot( * , 0 ,*))
y2=MIN(yp_drot( * ,Ny-1,*))
endelse
endelse
;% Define new grid
;% Machine-specific "Floating-point precision"
eps=(MACHAR(DOUBLE=SIZE(x1,/TYPE) eq 5)).eps
;% New pixel size
dxg=MIN(map.dx)
dyg=MIN(map.dy)
;% New image size
if KEYWORD_SET(full_size) then begin
Nxg=CEIL(((x2-x1)/dxg+1)*(1-2*eps))>1
Nyg=CEIL(((y2-y1)/dyg+1)*(1-2*eps))>1
endif else begin
Nxg=FLOOR(((x2-x1)/dxg+1)*(1+2*eps))>1
Nyg=FLOOR(((y2-y1)/dyg+1)*(1+2*eps))>1
endelse
;% New center
xc=0.5*(x1+x2)
x0=xc-0.5*(Nxg-1)*dxg
yc=0.5*(y1+y2)
y0=yc-0.5*(Nyg-1)*dyg
;% New grid
xpg0=(x0+dxg*FINDGEN(Nxg))#REPLICATE(1,Nyg)
ypg0=REPLICATE(1,Nxg)#(y0+dyg*FINDGEN(Nyg))
xcg0=GET_ARR_CENTER(xpg0)
ycg0=GET_ARR_CENTER(ypg0)
endif else begin
;% Just use current grid
UNPACK_MAP,map(0),NX=Nxg,NY=Nyg,XC=xcg0,YC=ycg0,DX=dxg,DY=dyg
xpg0=GET_MAP_XP(map(0))
ypg0=GET_MAP_YP(map(0))
endelse
endelse
;%%% New START times
tpg=REPLICATE(t0,Nxg)
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;%%% Define new map(s)
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;% Create template map
map_drot0=map(0)
data_drot=MAKE_ARRAY(DIM=[Nxg,Nyg],VALUE=map(0).data(0))
map_drot0=REP_TAG_VALUE(map_drot0,data_drot,'DATA') ; DATA
map_drot0.xc=xcg0 ; XC
map_drot0.yc=ycg0 ; YC
map_drot0.dx=dxg ; DX
map_drot0.dy=dyg ; DY
map_drot0=REP_TAG_VALUE(map_drot0,tpg,'START') ; START
rtime=ANYTIM2CAL(t0,FORM=9) ; RTIME
if TAG_EXIST(map_drot0,'RTIME') then $
map_drot0=REP_TAG_VALUE(map_drot0,rtime,'RTIME') $
else $
ADD_PROP,map_drot0,RTIME=rtime
if do_add_rdur then begin ; RDUR
rdur=DOUBLE(data_drot)
if TAG_EXIST(map_drot0,'RDUR') then $
map_drot0=REP_TAG_VALUE(map_drot0,TEMPORARY(rdur),'RDUR') $
else $
ADD_PROP,map_drot0,RDUR=TEMPORARY(rdur)
endif
if TAG_EXIST(map_drot0,'ROLL') then $ ; ROLL
map_drot0=REP_TAG_VALUE(map_drot0,roll(0),'ROLL') $
else $
ADD_PROP,map_drot0,ROLL=roll(0)
if TAG_EXIST(map_drot0,'ROLL_CENTER') then $ ; ROLL_CENTER
map_drot0=REP_TAG_VALUE(map_drot0,rcenter(*,0),'ROLL_CENTER') $
else $
ADD_PROP,map_drot0,ROLL_CENTER=rcenter(*,0)
;% Replicate template map to create the output map(s)
map_drot=REPLICATE(map_drot0,Nmaps)
map_drot0=0
map_drot.roll=roll
map_drot.roll_center=rcenter
;% Copy all other tags into the new map(s)
name=TAG_NAMES(map)
if do_add_rdur then begin
others=WHERE(name ne 'DATA' $
and name ne 'XC' $
and name ne 'YC' $
and name ne 'DX' $
and name ne 'DY' $
and name ne 'START' $
and name ne 'RTIME' $
and name ne 'RDUR' $
and name ne 'ROLL' $
and name ne 'ROLL_CENTER' $
,N_others)
endif else begin
others=WHERE(name ne 'DATA' $
and name ne 'XC' $
and name ne 'YC' $
and name ne 'DX' $
and name ne 'DY' $
and name ne 'START' $
and name ne 'RTIME' $
and name ne 'ROLL' $
and name ne 'ROLL_CENTER' $
,N_others)
endelse
for itag=0,N_others-1 do map_drot.(others(itag))=map.(others(itag))
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;%%% Differentially rotating map data
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=0,Nmaps-1 do begin
;% Rotate coordinates
if VALID_MAP(map_ref) eq 0 $
and (Nmaps gt 1 or KEYWORD_SET(same_grid) eq 0) $
and KEYWORD_SET(memsave) eq 0 then begin
xpr=xp_drot(*,*,i)
ypr=yp_drot(*,*,i)
endif else begin
xyr=DROT_COORD(map(i),t0,VERBOSE=verb,SPHERE=KEYWORD_SET(sphere) $
,KEEP=keep_angles)
xpr=xyr(*,*,0)
ypr=xyr(*,*,1)
xyr=0
endelse
;% Roll back coordinates to zero roll angle
if map_roll(i) ne 0 then $
ROLL_XY,xpr,ypr,-map_roll(i),CENTER=map_rcenter(*,i),xpr,ypr
;% Roll new coordinates to the final angle
xpg=xpg0
ypg=ypg0
if roll(i) ne 0 then $
ROLL_XY,xpg,ypg,-roll(i),CENTER=rcenter(*,i),xpg,ypg
;% Resample data
if verb then PRINT,'%I> DROT_RASTER: Resampling data array...'
map_drot(i).data=INTERP2D(map(i).data,xpr,ypr,xpg,ypg,[Nx,Ny] $
,REGULAR=0,_EXTRA=_extra)
;% Estimate rotation interval on a pixel-by-pixel basis;
;% N.B.: no fancy interpolation keywords used here (_extra keywords); only
;% using here keyword EXTRAPOLATE.
if do_add_rdur then begin
if multi_time then $
dt=t0-map(i).start $
else $
dt=REPLICATE(t0-ANYTIM2TAI(map(i).time),Nx)
dt=dt#REPLICATE(1,Ny)
map_drot(i).rdur=INTERP2D(dt,xpr,ypr,xpg,ypg,[Nx,Ny] $
,REGULAR=0,/EXTRAP)
endif
endfor
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;%%% End
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RETURN,map_drot
END