読者です 読者をやめる 読者になる 読者になる

どせいたんさき。

ナスダヨー

IDL 講習会(中級編)のメモ

先日参加した国立天文台の IDL 講習会(中級編)で作成したスクリプトなどをメモ.実習中にばばっと作成したものなのでコメントがひどかったり無駄があったりするかもしれない.改善点など教えていただけると助かります.

irc_read_fits

irc の fits 画像を読み込むためのプロシージャ. list オプションをつければ N?.lst からドドドッと読み込むことができる. fits 構造体の何番目を読み込むかをオプションで指定できるようにすればいろいろと捗るはず.

Pro irc_read_fits, filename, datacube, list=list
;; irc_read_fits.pro
;;
;; read irc imaging files
;;
  if ~keyword_set(list) then begin
     filelist = filename
  endif else begin
     print, 'if KEYWORD:list is set, filename is read as file list'
     print, '## under construction ##'
     readcol, filename, filelist, FORMAT='A'
  endelse

  nimg=n_elements(filelist)
  print, '% loading basic header information...'
  tmpimage = mrdfits(filelist[0],0,tmphdr,/silent)
  axis1=sxpar(tmphdr,'NAXIS1')
  axis2=sxpar(tmphdr,'NAXIS2')
  datacube=intarr(axis1,axis2,nimg)
  print, '% start read fits files.'
  for i=0, nimg-1 do begin
     ;; read 3d image from filaname
     tmpimage = mrdfits(filelist[i],0,tmphdr)
     ;; check header information
     print, 'FILENUM : ', i
     print, 'FILENAME: ', filelist[i]
     print, 'CHANNEL : ', sxpar(tmphdr,'CHANNEL')
     print, 'FILTER  : ', sxpar(tmphdr,'FILTER')
     ;; save 2d image from tmpimage to datacube
     datacube[*,*,i] = tmpimage[*,*,1]
  endfor
end

irc_subfits.pro

3 次元 fits array から特定の fits ファイルを引き算するためのプロシージャ.引数が string の時はファイルを読み込むようにしてある.バイアスやスカイを引くときに使った.

Pro irc_subfits, in_datacube, out_datacube, fitsfile, mask=mask
;; irc_subfits.pro
;;
;; subtract 2d image from 3d datacube
  ;; obtain array size
  ;; 3d_dim should be a vector contains three elements
  inputdim = size(in_datacube,/DIMENSIONS)
  nimg = inputdim[2]
  
  if ~keyword_set(mask) then mask = (intarr(inputdim[0],inputdim[1])+1)

  ;; read 2d fits file
  if size(fitsfile,/TYPE) eq 7 then begin
     ;; load 2d fits from file
     tmpfits = mrdfits(fitsfile)
  endif else begin
     ;; load 2d fits from array
     tmpfits = fitsfile
  endelse
  
  ;; creage datacube for subtract
  subcube = intarr(inputdim[0],inputdim[1],inputdim[2])
  for i=0, nimg-1 do begin
     subcube[*,*,i] = tmpfits
  endfor
  
  ;; subtract subcube from in_datacube
  out_datacube = in_datacube - subcube
end

irc_divfits

同じく 3 次元 fits array を引数にするプロシージャ.今度は割り算.フラットで割るときに使った.

Pro irc_divfits, in_datacube, out_datacube, fitsfile, mask=mask
;; irc_divfits.pro
;;
;; divide 2d image from 3d datacube
  ;; obtain array size
  ;; 3d_dim should be a vector contains three elements
  inputdim = size(in_datacube,/DIMENSIONS)
  nimg = inputdim[2]
  
  if ~keyword_set(mask) then mask=(intarr(inputdim[0],inputdim[1])+1)

  ;; read 2d fits file
  if size(fitsfile,/TYPE) eq 7 then begin
     ;; load 2d fits from file
     tmpfits = mrdfits(fitsfile)
  endif else begin
     ;; load 2d fits from array
     tmpfits = fitsfile
  endelse
  
  ;; creage datacube for divide
  divcube = intarr(inputdim[0],inputdim[1],inputdim[2])
  for i=0, nimg-1 do begin
     divcube[*,*,i] = (tmpfits * mask)
  endfor
  
  ;; divtract divide from in_datacube
  out_datacube = in_datacube / divcube
end

irc_imshift

shift and add をするためのプロシージャ. starpos は基準となる星のピクセル位置(整数).

Pro irc_imshift, in_datacube, out_datacube, starpos
;; irc_imshift.pro
;;
;; shift image using position of star
  ;; obtain array size
  ;; 3d_dim should be a vector contains three elements
  inputdim = size(in_datacube,/DIMENSIONS)
  nimg = inputdim[2]

  ;; calculate shift value

  ;; shift image
  out_datacube = in_datacube
  for i=0, nimg-1 do begin
     shift_d = starpos[*,i] - round(total(starpos,2)/nimg)
     out_datacube[*,*,i] = shift(in_datacube[*,*,i],$
                                -shift_d[0],-shift_d[1])
  endfor
end

irc_findstar

指定された領域を切り出して 2d-Gaussian でフィットして中心位置を検出するプロシージャ.結局サブピクセルでの shift and add をやっていないのでこの時点ではあまり意味がないスクリプトになっているかもしれない.

Pro irc_findstar, in_datacube, starpos, center, sigma=sigma, area=area
;; irc_findstar.pro
;;
;; find stars around center
  ;; obtain array size
  ;; 3d_dim should be a vector contains three elements
  inputdim = size(in_datacube,/DIMENSIONS)
  xmax = inputdim[0]-1
  ymax = inputdim[1]-1
  nimg = inputdim[2]

  if ~keyword_set(area) then area=15
  if ~keyword_set(sigma) then sigma=5s

  starpos = intarr(2,nimg)
  for i=0, nimg-1 do begin
     tmpcutout = in_datacube[((center[0,i]-area)>0):((center[0,i]+area)<xmax), $
                             ((center[1,i]-area)>0):((center[1,i]+area)<ymax),i]
     fitinfo = gauss2dfit(tmpcutout, fitres)
     starpos[*,i] = [fitres[4]+center[0,i]-area, $
                     fitres[5]+center[1,i]-area ]
  endfor

end

irc_combine

3 次元 fits array の足し合わせを実行するためのプロシージャ.一応 sigma-clipping 対応. median オプションをつけると通常使用されているアルゴリズムでの sigma-clipping になる.本来ならそれをデフォルトにするべきだけども……. x,y 方向にループを回していないので実行速度はまずまず.

Pro irc_combine, in_datacube, out_image, sigclip=sigma, median=median
;; irc_combine.pro
;;
;; combine 3d datacube to 2d image
  ;; obtain array size
  ;; 3d_dim should be a vector contains three elements
  inputdim = size(in_datacube,/DIMENSIONS)
  xsize = inputdim[0]
  ysize = inputdim[1]
  nimg  = inputdim[2]
  
  if (~keyword_set(sigma)) then sigma=!values.f_infinity

  ;; create datamask for nan from datacube
  mask = fltarr(xsize,ysize,nimg) + 1
  idx = where(finite(in_datacube,/nan),wflg)
  tmp_datacube = in_datacube
  if (wflg gt 0) then begin
     tmp_datacube[idx] = 0 & mask[idx] = 0
  endif

  ;; calculate mean value and create data cube
  mean_datacube = fltarr(xsize,ysize,nimg)
  ndata = total(mask,3)
  idx = where(ndata eq 0, wflg)
  if (wflg gt 0) then ndata[idx] = !values.f_nan
  mean_image = (total(tmp_datacube,3)/ndata)
  for i=0, nimg-1 do begin
     mean_datacube[*,*,i] = mean_image
  endfor

  ;; calculate deviation datacube
  ;; calculate variance datacube
  dev_datacube = in_datacube - mean_datacube
  var_datacube = fltarr(xsize,ysize,nimg)
  dof = total(mask,3) - 1.      ;degree of freedom
  idx = where(dof le 0, wflg)
  if (wflg gt 0) then dof[idx] = !values.f_nan
  for i=0, nimg-1 do begin
     var_datacube[*,*,i] = (dev_datacube[*,*,i])^2/dof
  endfor

  ;; if keyword_median is on
  if (keyword_set(median)) then begin
     median_image = median(in_datacube,DIMENSION=3)
     for i=0, nimg-1 do begin
        mean_datacube[*,*,i] = median_image
     endfor
     dev_datacube = in_datacube - mean_datacube
  endif

  ;; calculate sigma datacube
  sig_image    = sqrt(total(var_datacube,3))
  sig_datacube = fltarr(xsize,ysize,nimg)
  for i=0, nimg-1 do begin
     sig_datacube[*,*,i] =  sig_image
  endfor

  ;; create datamask for sigma-clipping
  mask = fltarr(xsize,ysize,nimg) + 1.
  idx = where(dev_datacube gt (sigma * sig_datacube), wflg)
  if (wflg gt 0) then mask[idx] = 0.

  ;; if masked, an element of tmp_datacube will be 0
  tmp_datacube = in_datacube
  idx = where(mask eq 0, wflg)
  if (wflg gt 0) then tmp_datacube[idx] = 0.

  ;; subtract subcube from in_datacube
  ndata = total(mask,3)
  idx = where(ndata eq 0, wflg)
  if (wflg gt 0) then npar[idx] = !value.f_nan
  out_image = (total(tmp_datacube,3)/ndata)
end

im_meanclip

meanclip という idlastro に含まれるプロシージャを利用して sigma-clipping を実行するためのプロシージャ. for loop で回した場合どれくらい遅くなるかをチェックするために作成.実際に計測したわけではないが 2048x2048 array x 10 でだいたい 10 倍位遅くなるみたい.

Pro im_meanclip, image, outimage, sigimage, sigclip=sigma, maxiter=maxiter
  s = size(image,/DIMENSIONS)
  xs = s[0]
  ys = s[1]
  im = s[2]
  
  if (~keyword_set(sigma)) then sigma = 3.
  if (~keyword_set(maxiter)) then maxter = 5

  outimage = fltarr(xs,ys)
  sigimage = fltarr(xs,ys)
  for i=0, xs-1 do begin
     for j=0, ys-1 do begin
        tmp_m = !values.f_nan
        tmp_s = !values.f_nan
        idx = where(finite(image[i,j,*]),wflg)
        if (wflg gt 0) then begin
           meanclip, image[i,j,*], $
                     tmp_m, tmp_s, $
                     clipsig=sigma,$
                     maxiter=maxiter
        endif
        outimage[i,j] = tmp_m
        sigimage[i,j] = tmp_s
     endfor
  endfor
end