;+
; PROCEDURE get_data_thmasi
;
; :DESCRIPTION:
; Extract the ASI data at the designated time from the tplot variable
; and make the structure.
;
; :PARAMS:
; asi_vn: tplot variable name (as string)
;
; :KEYWORDS:
; set_time: set the time (UNIX time) to get 2D data (can be an array)
; altitude: set the altitude on which the image data will be mapped.
; The default value is 110 (km).
; aacgm: set to obtain the position data in the AACGM coordinates
; data: output data as the structure which have the tags as follows;
; name: tplot variable name
; set_time: the designated time (array)
; dat_time: the actual time (array)
; data: image data
; alti: altitude (scalar)
; azim: azimuth angle for each pixel of the images
; elev: elevation angle for each pixel of the images
; center_glat: geographic latitude for the center of each pixel
; center_glon: geographic longitude
; corner_glat: geographic latitude for the corner of each pixel
; corner_glon: geographic longitude
; center_mlat: aacgm latitude for the center of each pixel
; center_mlon: aacgm longitude
; center_mlt: magnetic local time
; corner_mlat: aacgm latitude for the corner of each pixel
; corner_mlon: aacgm longitude
; corner_mlt: magnetic local time
;
; :AUTHOR:
; Yoshimasa Tanaka (E-mail: ytanaka@nipr.ac.jp)
;
; :HISTORY:
; 2014/07/08: Created
;
;-
pro get_data_thmasi, asi_vn, set_time=set_time, $
cal=cal, altitude=altitude, aacgm=aacgm, data=data
;----- initialize the map2d environment -----;
map2d_init
;----- output data -----;
data=''
;===== check parameters =====;
npar=n_params()
if npar lt 1 then return
;----- if asi_vn is the index number for tplot var -----;
vn = tnames(asi_vn)
print, vn
if total(vn eq '') gt 0 then begin
print, 'given tplot var(s) does not exist?'
return
endif
prefix = strmid( vn, 0, 8 )
dtype = strmid( prefix, 4,3 ) ; ast or asf
stn = strmid( vn, 8,4 ) ;3-letter station code
if strpos(dtype,'ast') eq 0 then is_thumb=1 else is_thumb=0
if is_thumb then begin
nx=32 & ny=32
endif else begin
nx=256 & ny=256
endelse
;----- set_time -----;
if ~keyword_set(set_time) then begin
t0 = !map2d.time
get_timespan, tr
if t0 ge tr[0] and t0 le tr[1] then begin
set_time = t0
endif else begin
set_time = (tr[0]+tr[1])/2. ; take the center of the designated time range
endelse
endif
;----- altitude -----;
if ~keyword_set(altitude) then altitude=110. ; 90, 110, 150, 250km
;----- aacgm -----;
if ~keyword_set(aacgm) then aacgm=!map2d.coord
;----- Load the cal file -----;
if keyword_set(cal) then calstr = cal[i] else $
thm_load_asi_cal, stn, calstr
;----- check altitude -----;
case is_thumb of
0: begin
idx = where( strpos( calstr.vars[*].name, vn+'_alti' ) eq 0 )
if idx[0] ne -1 then altvec = *(calstr.vars[idx[0]].dataptr)
ialt=where(fix(altvec/1000.) eq fix(altitude), cnt)
if cnt eq 0 then begin
print, 'no position data for the designated altitude!!!'
return
endif
alti=altitude
end
1: alti=110.
endcase
;----- obtain azimuth and elevation angle -----;
idx = where( strpos( calstr.vars[*].name, vn+'_elev' ) eq 0 )
if idx[0] ne -1 then elev = *(calstr.vars[idx[0]].dataptr)
idx = where( strpos( calstr.vars[*].name, vn+'_azim' ) eq 0 )
if idx[0] ne -1 then azim = *(calstr.vars[idx[0]].dataptr)
;----- obtain corner position data -----;
if not is_thumb then begin ;For asf
idx = where( strpos( calstr.vars[*].name, vn+'_glat' ) eq 0 )
if idx[0] ne -1 then corner_glat = reform( (*(calstr.vars[idx[0]].dataptr))[ialt, *, *] ) ;[257, 257]
idx = where( strpos( calstr.vars[*].name, vn+'_glon' ) eq 0 )
if idx[0] ne -1 then corner_glon = reform( (*(calstr.vars[idx[0]].dataptr))[ialt, *, *] ) ;[257, 257]
;----- obtain center position data -----;
center_glat=fltarr(nx, ny)
center_glon=fltarr(nx, ny)
for ix=0, nx-2 do begin
for iy=0, ny-2 do begin
center_glat[ix, iy]=mean(corner_glat[ix:ix+1, iy:iy+1])
center_glon[ix, iy]=mean(corner_glon[ix:ix+1, iy:iy+1])
endfor
endfor
endif else begin ;For ast
idx = where( strpos( calstr.vars[*].name, vn+'_glat' ) eq 0 )
if idx[0] ne -1 then corner_glat = *(calstr.vars[idx[0]].dataptr) ;[4, 1024]
idx = where( strpos( calstr.vars[*].name, vn+'_glon' ) eq 0 )
if idx[0] ne -1 then corner_glon = *(calstr.vars[idx[0]].dataptr) ;[4, 1024]
;----- obtain center position data -----;
center_glat=fltarr(nx*ny)
center_glon=fltarr(nx*ny)
for ipxl=0L, nx*ny-1 do begin
center_glat[ipxl]=mean(corner_glat[*, ipxl])
center_glon[ipxl]=mean(corner_glon[*, ipxl])
endfor
endelse
;----- obtain mlat mlon -----;
corner_mlat='' & corner_mlon=''
center_mlat='' & center_mlon=''
if keyword_set(aacgm) then begin
;----- Load the S-H coefficients -----;
ts = time_struct(set_time)
aacgmloadcoef, ts.year
;----- corner mlat mlon -----;
altmat = corner_glat & altmat[*] = alti ;***** altitude in km *****;
aacgmconvcoord, corner_glat, corner_glon, altmat, $
corner_mlat, corner_mlon, err, /to_aacgm
;----- center mlat mlon -----;
altmat = center_glat & altmat[*] = alti ;***** altitude in km *****;
aacgmconvcoord, center_glat, center_glon, altmat, $
center_mlat, center_mlon, err, /to_aacgm
endif
;----- initialize output arrays -----;
ntime=n_elements(set_time)
set_time_all = ''
dat_time_all = ''
image_all = fltarr(ntime, nx, ny)
if keyword_set(aacgm) then begin
center_mlt_all = fltarr(ntime, nx, ny)
corner_mlt_all = fltarr(ntime, nx+1, ny+1)
endif else begin
center_mlt_all = ''
corner_mlt_all = ''
endelse
for itime=0L, ntime-1 do begin
stime=time_double(set_time[itime])
;----- obtain image data for the designated time -----;
get_data, vn, data=d
tidx = nn(d.x, stime)
if tidx lt 0 then continue
image = reform(d.y[tidx, *, *])
dtime=d.x[tidx]
if is_thumb then begin ;array rotation mimicing Line 173 in thm_mosaic_array.pro
image = rotate( image, 8 )
endif else begin
bkgd = mean( image[0:10,0:10] ) ;Define the background count by averaging counts near the bottom-left corner
image_sbtrctd = image - bkgd ;Subtraction of background count
image = image_sbtrctd
endelse
;----- check if data for the designated time is obtained or not. -----;
crt_dt=2.
dt = abs(stime - d.x[tidx])
if dt lt crt_dt then note = ' (ok)' else note = ' !!! not within '+string(fix(crt_dt))+' sec !!!'
print, '========== '+vn+' =========='
print, 'designated time: '+time_string(stime)
print, ' ASI data time: '+time_string(dtime), tidx, note
d = 0L ;initialize the variable to save the memory
;----- append array -----;
append_array, set_time_all, stime
append_array, dat_time_all, dtime
image_all[itime, *, *] = image
;----- ontain mlt -----;
if keyword_set(aacgm) then begin
ts = time_struct(stime)
;----- center mlt -----;
yrs = fix(center_mlat) & yrs[*] = ts.year
yrsec = long(center_mlat) & yrsec[*] = long( (ts.doy-1)*86400. + ts.sod )
mlt = aacgmmlt( yrs, yrsec, center_mlon )
mlt = ( ( mlt + 24. ) mod 24. ) / 24.*360. ; [deg]
igt=where(mlt gt 180., cnt)
if cnt gt 0 then mlt[igt] -= 360.
center_mlt_all[itime, *, *] =mlt
;----- corner mlt -----;
yrs = fix(corner_mlat) & yrs[*] = ts.year
yrsec = long(corner_mlat) & yrsec[*] = long( (ts.doy-1)*86400. + ts.sod )
mlt = aacgmmlt( yrs, yrsec, corner_mlon )
mlt = ( ( mlt + 24. ) mod 24. ) / 24.*360. ; [deg]
igt=where(mlt gt 180., cnt)
if cnt gt 0 then mlt[igt] -= 360.
corner_mlt_all[itime, *, *] =mlt
endif
endfor
data={name:vn, set_time:set_time_all, dat_time:dat_time_all, $
data:image_all, alti:alti, azim:azim, elev:elev, $
center_glat:center_glat, center_glon:center_glon, $
corner_glat:corner_glat, corner_glon:corner_glon, $
center_mlat:center_mlat, center_mlon:center_mlon, center_mlt:center_mlt_all, $
corner_mlat:corner_mlat, corner_mlon:corner_mlon, corner_mlt:corner_mlt_all}
end