• kyoto_load_dst
• kyoto_load_ae
• ace_swe_load
• ace_mfi_load
• omni_hro_load
• time_double(), time_string()
• time_struct()
• CDF_TT2000形式の時刻データの取り扱いについて (メモ)
• tplot_save, tplot_restore
• tplot_ascii
• print_tinfo
• add_data (ex, add_data, 'data1', 'data2', newname='data3' )
• dif_data
• mult_data
• div_data
• avg_data
• deriv_data
• timespan, get_timespan
• tlimit
• tplot_options
• options
• ctime
• ylim, zlim
• popen, pclose
• makepng
• cdf2tplot
• print_cdf_info
• split_vec, join_vec
• clean_spikes
• tinterpol_mxn
• time_average
• tsmooth_in_time
• thigh_pass_filter (e.g., thigh_pass_filter, ‘tplot_var’, 30. )
• tsub_average
• tclip, tdeflag, tdegap
• nn()
• append_array
• bin1d
• bin2d
• plotxyz
• scat_plot
• tdpwrspc, wav_data
• tkm2re
• xyz_to_polar
• geopack_sphcar
• cotrans
• geopack_conv_coord
• thm_cotrans
• aacgm_conv_coord, aacgm_mlt(), aacgm_mlong(), aacgm_load_coef
• thm_mom_fx
• thm_mom_plot
• thm_local
• thm_ucb
• istp_local
• istp_cdaweb
• kyoto_load_symasy

Quick-look DstデータをWDC/Kyotoからダウンロード、kyoto_dst というtplot変数として読み込む。

Provisional AU/AL/AO/AE をダウンロード、それぞれtplot変数に格納。デフォルトではAE indexのみ。kyoto_load_ae, datatype=’au al ae’ のようにdatatypeキーワードで欲しい指数を指定すればそれらを読み込むことができる。datatype=’all’ を指定すると全て(au, al, ae, ao) のデータを読み込む。

ACE SWE data(太陽風)をCDAWeb siteからダウンロード、tplot変数に格納。Keyword無しだとkey parameterをもってくる。ace_swe_load,datatype=’h0’ と実行すればL2 dataをもってくる。

ace_swe_loadのMFI (IMF) 版

OMNI2 data (地球前面のbow shockまでtime shiftさせた太陽風、IMF データ、及び各種地磁気指数、太陽黒点数など)をダウンロードして、tplot変数に格納。/res5min というkeywordを付けると5 min value (Tsyganenko modelの計算によく使われる)、デフォルトは1 min valueをもってくる。

TDASを使っていると、tplot変数の中の時刻表記に使う倍精度小数のUnix timeの計算や、通常の時刻形式 (例えば '2010-04-03/13:30:00' などの文字列)との間の変換を頻繁に行うことになる。そのためには time_double(), time_string(), time_struct() を使いこなす必要がある。

time_double() はTDASでの標準的な日時の表現である、倍精度小数で表したUnix time、つまり”1970年1月1日00:00:00 UT (1970-01-01/00:00:00 TAI)からの積算秒数を倍精度小数であらわしたもの” に変換する関数。普通は ’YYYY-MM-DD/hh:mm:ss’ のような文字列を引数とするが、time_struct()で生成されるTDASのtime構造体も引数にとることができる。ベクトル化されているので、時刻文字列を配列で与えると、それぞれについて計算を行い、結果を同じ数の配列で返す。またEPOCHキーワードをセットすると、入力値としてCDF_EPOCH形式を取ることもできる。

例) 
unix_time = time_double( ‘2008-02-02/16:30:00’ ) 
    unix_time = time_double( ['2008-02-02/00:00','2008-02-02/12:00'] ) 
    unix_time = time_double( cdf_epoch_time, /epoch )

time_string() は上記のtime_double()が吐く倍精度小数のUnix time、またはtime_struct()関数が吐くtime構造体を入力として、実際の日時を表す文字列(例えばYYYY-MM-DD/hh:mm:ss) を返す関数。文字列のフォーマットはTFORMAT keywordで明示的に指定することもできる。

例) 
dblt = time_double( ‘2007-12-14/15:00’ )
print, time_string( dblt ) 
->  2007-12-14/15:00:00 が出力される
    print, time_string( dblt, tformat='YYYYMMDD' ) 
      -> 20071214 が出力される

time_double()などで計算される倍精度小数Unix timeやTDAS標準の文字列による時刻表記(’yyyy-mm-dd/hh:mm:dd’)を引数にとり、TDAS特有のtime構造体を返す関数。Time構造体(ここではこう呼ぶことにする、正式名は不明)の中身は以下。

IDL> ts=time_struct( time_double('2007-03-23'))
IDL> help,ts,/str

  YEAR            INT           2007            ; 4-digit year
  MONTH           INT              3            ; month
  DATE            INT             23            ; day
  HOUR            INT              0            ; hour
  MIN             INT              0            ; minute
  SEC             INT              0            ; second
  FSEC            DOUBLE           0.0000000  ; fractal seconds
  DAYNUM          LONG            732757      ; day since 0 AD
  DOY             INT             82            ; day of year
  DOW             INT              4            ; day of week
  SOD             DOUBLE           0.0000000  ; second of day
  DST             INT              0          ;if 1, daylight saving time
  TZONE           INT              0         ; timezone (e.g. PST: -8)
  TDIFF           INT              0         ; hours from UTC

変換処理はベクトル化されており、引数に配列を与えることもでき、高速に動く。TDASの時刻形式からyear, …, secなどの日時を求めたり、データ解析でよく使われる day of year (上記のtime構造体の DOY) や second of day (同 SOD) を計算するのにも使われる。

Van Allen Probes衛星より、CDFデータフォーマットの新しい時刻形式であるCDF_TT2000形式が採用され始めた。ERG衛星データのCDFファイルでもこの時刻形式を採用している。この実体は、64ビット整数(LONG64)でJ2000.0 (2000-01-01/12:00:00 TT) からのnanosecond換算での経過時間によって時刻を表したものである。またIDL用の最新バージョンのCDF dynamic link moduleをインストールすることで、leap secondの補正にも対応する。

CDF_TT2000形式の時刻をUnix time に変換するには、time_double()にTT2000キーワードをセットすればよい。また変換を行う前にcdf_leap_second_init を実行することで、最新のleap second table をNASA/CDFサイトからダウンロードした上で必要な環境変数の設定が自動的に行われる。

IDL> cdf_leap_second_init 
IDL> unix_time = time_double( tt2000_time, /TT2000 ) 

なおtime_string(), time_struct() は直接CDF_TT2000形式からの変換に対応していないので、まずtime_double() でUnix timeに変換した後で、これらの関数に渡すようにする。

また、独自にCDFファイルを生成する時など、逆に通常の時刻形式からCDF_TT2000形式に変換したい場合は、"cdf_tt2000" コマンドを用いる。

2014-03-17/12:00:00 UTC を変換するには
IDL> cdf_tt2000, tt2000_time, 2014, 3, 17, 12, 0, 0, /compute_epoch 

逆変換するにはbreakdown_epochキーワードを使う
IDL> cdf_tt2000, tt2000_time, yyyy, mm, dd, hr, min, sec, msec, nsec, /breakdown_epoch 

tplot_saveはtplot変数の中身をファイルにセーブする。Tplot_restoreはそのファイルを読み込んで、tplot変数をIDLメモリー上に回復する。

tplot_save, ‘thd_esa_en_eflux’, filename=’thd_esa’
    - >	thd_esa.tplot というファイルにセーブされる

tplot_restore, filename=’thd_esa.tplot’  

また同じtplot変数名だが異なる時刻範囲が保存されたtplot saveファイルを一度にリストアすることができる(appendキーワード)。

tplot_restore, file=[ 'thd_esa_20141101.tplot', 'thd_esa_20141102.tplot' ], /append

tplot変数を指定しないでtplot_saveを実行すると、IDLメモリー上にある全てのtplot変数をセーブしようとするので注意。

tplot変数の中身をファイルにascii dumpする。

Ex) tplot_ascii, ‘tha_state_pos’ 
->	tha_state_pos.txt　というファイルができる

tplot変数の中のデータ変数やメタデータを簡易表示する。timeキーワードをセットすると、データの先頭・最後の時刻ラベルも表示される。

IDL> print_tinfo, 'OMNI_HRO_1min_BZ_GSM', /time
*** Variable: OMNI_HRO_1min_BZ_GSM
Start time: 2017-06-01/00:00:00
End time: 2017-06-30/23:59:00
** Structure <16c9978>, 2 tags, length=518400, data length=518400, refs=1:
   X               DOUBLE    Array[43200]
   Y               FLOAT     Array[43200]

tplot変数同士(data1, data2) の足し算をして、結果を新しいtplot変数 (data3)に格納する。注意する点は、data1とdata2が同じ次元のデータで、かつデータ数(つまり時刻ラベルの数)が同じである必要がある。

tplot変数同士の引き算。使い方および注意点はadd_data と同様。

tplot変数の掛け算。

tplot変数の割り算。

時間平均を計算して新しいtplot変数に格納する。

avg_data,'thc_fgs_gsm', 120   ;;Th-Cの磁場データから120秒平均を計算する
  -> 結果はthc_fgs_gsm_avg というtplot変数に格納。もちろん格納するtplot変数名を指定することもできる。

ただしTDASに最初から入っているavg_data.proは非常に処理が遅いため、数万個レベルの配列(例えば1秒値で1日分くらいのデータ)が入ったtplot変数に対して使うのは避けた方がよい。もっと高速で動く改造版がmember_contribリポジトリの misc/avg_data.pro にある。

tplot変数の時間微分をとって、新しいtplot変数に格納する。THEMIS gmag dataからDC成分を除いて変動成分のみ(Pi2など)を見る時に便利。

Timespanは時間幅を指定する。ここで指定した時間幅はthm_load_xxx等のデータをロードするprocedure, またtplotなどのプロットprocedureの時間幅となる。Get_timespanは、timespanで指定した時間幅(表からは直接見えない構造体に格納、保持される)を知るために使う。

THEMIS> timespan, ‘2004-03-12/01:00:00’, 2, /day    
;;2004年3月12日1時0分 UTから2日分 (3/14 01:00 UTまで)を指定

THEMIS> get_timespan, t   ; t は2成分の倍精度小数の配列になる 
THEMIS> Print, time_string(t)
2004-03-12/01:00:00 2004-03-14/01:00:00 

tplotで表示されたプロットの任意の点を開始時間、終了時間の順で左クリックを行うことで、その2点で挟まれた時間幅で再プロットする。インタラクティブにプロットの時間幅を変えたい場合に非常に便利。

Ex) 	tlimit      ;; マウスクリックで指定した時間幅のプロットにする
	tlimit, /last  ;; １つ前の時間幅に戻す
	tlimit, /full   ;; timespanで指定されている全時間幅に戻す

開始時間を左クリックした後でそれ以前の時間を左クリックするとtlimit, /lastと同じ動作をする。また、開始時間を左側の余白(現在のプロットの開始時刻より左)、終了時間を右側の余白(現在のプロットの開始時刻より右)で順にそれぞれ左クリックすると、tlimit, /fullと同じ動作をする。これら機能を利用するとわざわざtlimit, /lastや/fullを使わずにtlimitとマウスクリックのみでプロットの時間幅の拡大縮小を繰り返すことが可能になる。

MacOSではこのtlimitでのマウスクリックがうまく動かないことがある。workaround がTDASのクリックリファレンスに載っている。 http://themis.ssl.berkeley.edu/downloads/THEMIS_Software_Users_Quick_Reference_Guide_v8_00.pdf

Tplotでプロットする際の、プロット全体に関する様々な描画オプションを設定する。よく使うのはtitleとvar_label。

Ex) tplot_options, ‘title’, ‘ESA-SST-FGM summary plot’  ;;プロットの1番上にタイトル
   ;;プロットの1番下に衛星軌道のラベルを付加
   Tplot_options, ‘var_label’, ‘tha_state_pos_x’ 

Options, 'tplot変数名', 'option名', optionの内容 という使い方で、tplotでプロットする際の、各tplot変数のプロットに関する様々な描画オプションを設定する。Tplot変数名には*, ?などのワイルドカードも使用可能。また変数名の代わりに、tplot変数に割り当てられている番号を使ってもよい。

Ex) options, ‘tha_fgm_gsm’, ‘ytitle’, ‘B [nT]’   
;; tha_fgm_gsmをプロットすると縦軸に B [nT] というラベルが付く

tplotで表示されたプロットの任意の点をクリックすることで、その点の時刻と縦軸の値等を変数に格納する。

ctime, time, y   ;;time: UNIX timeでの時刻、y: 縦軸の値

またctime.proに類似したものとして，crosshairs.proがある。これはプロットの任意の点をクリックすることで，その点のnormal座標が変数に格納される。

tplot変数を引数に取り、それらがtplotでプロットされる際のy軸、z軸のプロット範囲を指定する。[yz]lim, ‘tplot変数名’, min, max, [0 or 1]　というふうにして使う。tplot変数名はワイルドカードが使えるし(例: ‘th?_fgs’)、変数名を文字列配列として複数同時に与えることもできる(例: [‘tha_fgs’,’thb_fgs’])。また４つ目の引数は0にすると linearスケール、1にするとlogスケールでのプロットになる。minとmaxに両方ゼロを与えると全データ値を含むような範囲に自動設定される。このylim, zlimを実行した後、再びtplotを実行した時にプロット範囲の変更が反映される。

ylim, ‘thg_mag_fsim’, -400,400, 0   
;;[-400,+400]の範囲をlinearスケールでプロットするように設定

Tplot等で描画されるプロットをポストスクリプトファイルに出力するためのprocedure。使い方はまずpopenの後ろに出力PSファイル名を付けて実行し、tplot等の描画ルーチンを実行。最後にpcloseで出力ファイルを閉じる。

Ex) popen, ‘plot1.ps’    ;;plot1.ps というファイルを開ける
    tplot              ;;直前にプロットした内容で再描画
    pclose            ;;plot1.psを閉じる

landscape キーワードをオンにする ( /landscape のように書く)とデフォルトのPortrait (縦長)ではなく横長のPSファイルになる。

現在指定されているプロットウィンドウをPNGファイルとして出力する。

tplot, [‘tha_fgs_x’, ‘tha_fgs_y’,’tha_fgs_z’]
makepng, ‘tha_fgs_plot’
-> tha_fgs_plot.png というファイルに出力される。pngという拡張子が不要であることに注意。

Common Data Format (CDF)ファイルを読み込んで、時間変動する変数をtplot変数に格納する。

Ex) cdf2tplot, ‘CDFファイル名’, varnames=’Vx’  
;; Vxという名前の変数をtplot変数として読み込む

cdf2tplot, ‘CDFファイル名’,/all ;CDFファイル内の全ての変数をtplot変数にする

CDFファイルに格納されているメタデータ(g/v-属性)と、データ変数の配列構造などの一覧を画面に表示する。CDFファイルの中身のチェックをする際に非常に便利。

 IDL> print_cdf_info, 'omni_hro_1min_20170601_v01.cdf'
 File: /Users/horit/work/data/istp/omni/omni_cdaweb/hro_1min/2017/omni_hro_1min_20170601_v01.cdf
 INFO = STRUCT   =  --(6 Tags/5408 Bytes)-->
    FILENAME     = STRING   =       
     'omni_hro_1min_20170601_v01.cdf'
 INQ          = STRUCT   =  --(9 Tags/80 Bytes)-->
 NDIMS    = LONG     = 0
 DECODING = STRING   = 'HOST_DECODING'
 ENCODING = STRING   = 'NETWORK_ENCODING'
 MAJORITY = STRING   = 'ROW_MAJOR'
 MAXREC   = LONG     = -1
 NVARS    = LONG     = 0
 NZVARS   = LONG     = 47
 NATTS    = LONG     = 46
 DIM      = LONG[1]  = [0]
 G_ATTRIBUTES = STRUCT   =  --(27 Tags/736 Bytes)-->
  PROJECT                    = STRING   = 'NSSDC'
  DISCIPLINE                 = STRING   = 'Space Physics>Interplanetary Studies'
  SOURCE_NAME                = STRING   = 'OMNI (1AU IP Data)>Merged 1 minute Interplantary OMNI data'
  DATA_TYPE                  = STRING   = 'HRO>Definitive 1minute'
  DESCRIPTOR                 = STRING   = 'IMF and Plasma data'
  .... ....
  .... ....
 VARIABLES:
       NAME NUM IS_ZVAR  DATATYPE TYPE NUMATTR NUMELEM RECVARY NUMREC NDIMEN D[0] D[1] D[2] D[3] D[4] D[5]
      Epoch   0       1 CDF_EPOCH    5      -1       1       1  43200      0    0    0    0    0    0    0
         YR   1       1  CDF_INT4    3      -1       1       1  43200      0    0    0    0    0    0    0
        Day   2       1  CDF_INT4    3      -1       1       1  43200      0    0    0    0    0    0    0
         HR   3       1  CDF_INT4    3      -1       1       1  43200      0    0    0    0    0    0    0
     Minute   4       1  CDF_INT4    3      -1       1       1  43200      0    0    0    0    0    0    0
   .... ....
   .... ....

G-属性名＋中身(ただし長い文字列は途中でtruncateされる)に引き続いて、各データ変数の名前、変数型(DATATYPE),時刻方向の配列数(NUMREC)、配列の次元(NDIMEN)の一覧が表示される。端末を横長にしておくと後半のデータ変数一覧表示が変な位置で改行されず見易くなる。

なお/print_att を付けて実行すると、データ変数一覧のあとに、1) g-属性の完全リスト、2) v-属性のリストも表示される。

ベクトル量がデータとして入っているtplot var. に対して、成分ごとに別々のtplot var.を作って代入する。

Ex) split_vec, ‘tha_fgm_gsm’   ;;磁場データを成分分解
->	tha_fgm_gsm_x, tha_fgm_gsm_y, tha_fgm_gsm_z が生成される

またsplit_vecの反対で、同じ型のデータを持つtplot変数のデータ配列を１つにまとめて新しいtplot変数にする join_vec というコマンドもある。

 IDL> join_vec, 'tha_fgm_gsm_'+[ 'x', 'y', 'z' ], 'tha_fgm_gsm_new'

tplot変数を引数に取り、値のスパイクを取り除く(NaNで埋める)。実際には、ある幅でスムージングされた値と元データ値との差をとり、それが閾値を超えたかどうかで、スパイクかそうでないかを判断している。 以下の例ではthg_fsmi_x のデータを時刻方向に3配列分(0.5秒値なら1.5秒bin、配列数をns というkeywordで指定)の幅でrunning averageをとることでスムージングし、元データの値をx, スムージングされた値をxsとすると、 |x| > nsft|xs| / (ns-1+ft) という条件を満たす場合に、スパイク(つまりエラー値)と判定する。ftもkeywordで指定できる。またスパイクが除去されたデータ値はnew_name keywordで指定される新しいtplot変数に格納される。

Ex) clean_spikes, ‘thg_fsmi_x’, ns=3, ft=10, new_name=’thg_fsmi_x_cln’ 

IDL procedure “interpol” の多次元版で、tplot var.を入力にとる。時刻ラベルが荒いデータが入ったtplot var.を内挿して、時刻ラベルが細かいデータを生成する。例えば通常1分値の衛星軌道データ(以下の例中のtha_state_pos)を内挿して、磁場データ(tha_fgm_gsm, 3秒値)の時刻ラベルに合わせた軌道データを生成したい場合は次のようにする。

Ex) tinterpol_mxn, ‘tha_state_pos’, ‘tha_fgm_gsm’, newname=’tha_state_pos_intpl’

1つ目の引数のtplot var.が内挿されるデータ。2つ目の引数のtplot var. の時刻ラベルに合わせて内挿値が計算される。つまりnewnameで指定するtplot var.に入るデータは、2つ目の引数のtplot var.と全く同じ時刻ラベルが付く。

Tplot var.ではなく、IDLの普通の配列に対して、time binごとの平均をとって、結果の配列を返す。かなり便利な上に、内部でベクトル化されたルーチンを使っているので処理も速い。他にもoptionがあるが、基本的には以下のようにして使う。d.x, d.yをそれぞれget_data で生成された時系列の時刻ラベル配列、及びデータ配列とすると

newdata = time_average(d.x, d.y, newtime=newtime, trange=[‘2008-02-02/00:00’,’2008-02-02/12:00’], resolution=60.) 

これでtrangeで指定したtime interval 中で60 sec (1 min) bin ごとのデータの平均値を入れた配列が newdata に代入される。また”newtime” keyword で指定した変数に、bin分けした際の各binの時刻(bin中央の時刻)が入った配列が代入される。元のデータより細かいresolution (binの時間幅)を指定しないこと。 なおtplot変数を引数に取り、このprocedureと同じ処理をするのがavg_data である。

引数のtplot var.の時系列データに対してrunning averageを計算し、結果を別のtplot var.として生成する。Running averageの結果は、平均する前のtplot変数と同じ時刻タグを付けて返される。

Ex) tsmooth_in_time, ‘tplot_var’, 60., newname=’new_tplot_var’ 
60 sec windowでrunning averageをとった結果が new_tplot_var に入る。

Tplot var.に対してHigh pass filterを適用して、結果を新たなtplot var.に格納する。2つ目の引数でband filterの下限値(単位は秒)を指定する。↑の例だと30 sec windowでrunning averageをとったものを元データが差し引くことにより、30秒より速い変動成分のみを抜きだす。

データの全時間平均値を差し引いて、結果を新しいtplot var.に代入する。例えば地磁気データのゼロ線をそろえるのに使う。median keyword をオンにすると、平均ではなくmedianを差し引く。

Ex) tsub_average, ‘thg_mag_fsim’, out_name=’thg_mag_fsim_subtracted’

Tplot変数の中のエラー値を除去して、時間的にデータが連続になるように補間する、という処理でよく使われる。これらはセットで使われることが多い。具体的には、tclipはあるmin, maxの範囲外にある値をNaNと置き換える。tdegapはmaxgap keywordで指定される時間幅以内の時刻の抜けについて、抜けている時刻のデータをNaNとしてデータ配列に挿入する。さらにtdegapはNaNの値を周囲の時刻のデータから補間する。例えばthg_gill_xという名前のtplot変数で表されるデータには時刻抜けがあり、また正常値の範囲(以下では-1000 nT ~ +1000 nT)から外れているデータ値を含んでおり、それらを除去する場合、

tclip, ‘thg_gill_x’, -1000., +1000., /overwrite 　
  ;; [-1000,1000]の範囲外の値をNaNで置き換える。/overwriteを指定すると元のtplot変数に上書き

tdegap, ‘thg_gill_x’, maxgap=20, /overwrite 
;; 20時刻配列分、ここでthg_gill_xが0.5秒値なら10秒以下の時刻抜けについて、NaNという値の配列を挿入する。Maxgapを指定しないと、全ての時刻抜けをNaNで埋める。

tdeflag, ‘thg_gill_x’, ‘linear’, maxgap=20, /overwrite 
;;NaNまたはInfとなっている値を補間する。２つ目の引数に’linear’ を指定すると線形補間、’repeat’ を指定すると直前の正常値を繰り返すことで補間する。

また以上の例では/overwrite を指定してもとのtplot変数に上書きしたが、newname というkeywordで新しいtplot変数の名前を指定すれば、その変数名で結果が格納される。

tplot変数、またはget_dataで取り出したデータ構造体(x:Unix time形式の時刻、y:データ)を引数にとり、指定した時刻に1番近い時刻タグの配列番号を返す関数。

idx = nn( data, time ) 
;;data: データ構造体、またはtplot変数名
;;time: Unix time形式の時刻。データ構造体の中で、この時刻に1番近い
;;データの配列番号がidxに返る。

使用例としては、例えばdataがOMNI2データ(time shiftを考慮した太陽風・IMFデータ、1分値か5分値)として、ある時刻(time)でのOMNI2データを知りたい時に、

idx = nn( data, time )
print, data.y[idx] 

とすればいい。 このtimeは時刻を配列で与えることも可能。idxは同じ数の配列が返る。処理がベクトル化されているので、forループで回して1回1回nn()を呼ぶよりはるかに高速。

ある配列(または構造体)を別の配列(別の構造体)の最後に継ぎ足すためのprocedure。ある配列をa0, それに継ぎ足す配列をa1 とすれば、

append_array, a0, a1 

とすればいい。第1次元の方向に継ぎ足しされるので、第2次元かそれより高次の次元数はa0とa1とで一致している必要がある。また構造体同士の場合、メンバーの種類、数、型が一致している必要がある。これはIDLでの a0 = [a0, a1] として行われる処理と同値である。a0 = [a0,a1]と比較して、append_arrayではa0が空配列でもエラーにならない。例えば、

a0=’’ 
append_array, a0, a1 

とすると、a0はa1と等しくなる。配列を次から次へと継ぎ足していくような処理をスクリプトで書いた場合、a0が空かどうかのチェックをしなくていいので、スクリプトをすっきり書ける。

ある時系列のデータ配列を、それと同時刻数の1次元配列の値でbin分けし、各binでの平均、標準偏差、medianなどを計算してくれる、超便利procedure。

具体例としては、 ni: イオンの密度の時系列（1次元配列） ygsm: 同じ時系列での衛星の位置のY座標 (niと同じ数の1次元配列、単位はRe) の時に、

bin1d, ygsm, ni, -15., 15., 2., binnm, ygsmc, ave, stdv, med 

とすると、ygsmのうち[-15, 15]の範囲を2 Reごとのbinにわけて、それぞれのbinに入るデータ数がbinnmに, そのbinの中心でのygsmの値がygsmc, 各binでのniの平均、標準偏差、medianがそれぞれave, stdv, med という配列に返る。Binnm, ygsmc, ave, stdv, medは同じ数の1次元配列となる。

bin1dの2次元版。bin1dと同じで、内部処理がベクトル化されているため非常に高速で処理されるのが特徴。 Xgsm, ygsm, niは共通の時刻タグを持つ時系列で、それぞれ衛星のX座標 [Re]、Y座標 [Re]、イオン密度を表す1次元配列である時、

bin2d, xgsm, ygsm, ni, binsize=[1,2], $
    xrange=[-10,-20],$ yrange=[-15,15], $
    xcenters=xc, ycenters=yc, averages=ave, $
    medians=med, stdevs=stdv, binhistogram=binnm

とすると、xrangeの範囲でxgsmを1 Reごと、yrangeの範囲でygsmを2 Reごとにわけた2次元binを作り、そのbinに入るデータ数をbinnm、各binでのniの平均、median、標準偏差をそれぞれave, med, stdvに格納する。また各binの中心のX,Y座標はそれぞれxc, yc に入る。binnm, ave, med, stdv は同じ数の2次元配列になり、xc, ycはbinnmなどの1次元目、及び2次元目と同じ数の1次元配列となる。 また上記の例で、binsizeの代わりに、binumというkeywordを与えると、bin幅ではなく、全体で幾つのbinに切るかを指定できる。例えばbinum=[20,30] を与えると、xgsmをxrangeの範囲で20分割、ygsmをyrangeの範囲で30分割した、20x30の2次元binにわけて、平均、標準偏差等を計算する。

データから2次元平面上にcolor contourを描画するprocedure。bin2dとセットで使われることが多く、bin2dから返る配列をそのままplotxyzの引数に与えることができる。例えば、

plotxyz, x, y, z 

が基本形で、ここでzはcontourとして表される2次元配列で、これがM×Nの配列とすれば、xはMの次元の座標を表す1次元配列、またyはNの次元の座標を表す1次元配列となる。それぞれ、M個、N個の1次元配列である。つまり上記のbin2dの例の返り値を使って、

plotxyz, xc, yc, ave 

とやれば、X-Y平面でのniのcolor contourを簡単に表示することができる。

なおplotxyzには描画オプションがいろいろある。詳細については (tdas dir)/idl/themis/examples/thm_crib_plotxyz.pro を参照のこと。

scatter plotを描画する。２つのtplot変数を与え、1つ目が横軸、2つ目が縦軸の値として、scatter plotが作られる。なお２つのtplot変数は要素数が同じである必要がある。

IDL> timespan, ‘2008-11-22’
IDL> thm_load_esa, probe=’c’, level=2, $
datatype=’peir_density peir_avgtemp’
IDL> scat_plot, ‘thc_peir_density’, ‘thc_peir_avgtemp’ 

どちらもtplot変数を引数にとり、tplot変数内の時系列データのdynamic 周波数spectrumを計算し、newname keywordで指定される新しい名前のtplot変数に結果を格納する。ただし前者はFFT、後者はWavelet変換を用いる。wav_dataのみ、trange keywordで求めたい時間幅を与えるか、maxpoints keywordで最大の点数を指定する必要がある。そうしないとデータ点数が多い場合、マウスクリックで時間幅を指定させられることになる。maxpoints=2l^22などど非常に大きな値を指定しておくとメモリがいっぱいになる点数までは時間幅を指定せずとも計算を行ってくれる。

Ex) 
split_vec, ‘thg_mag_fsmi’
tdpwrspc, ‘thg_mag_fsmi_x’   ; X成分のみFFTでdynamic spectrum求める 
　　　　　　　　　　　　　　; 結果はthg_mag_fsmi_x_dpwrspc というtplot変数に格納
wav_data, ‘thg_mag_fsmi_x’, $
trange = time_double(‘2004-09-10/00:00’)+[0.,86400.]
;; 2004-09-10 00:00 – 2004-09-11 00:00の範囲をwavelet変換
       ;; 結果のdynamic spectrumはthg_mag_fsmi_x_wv_powというtplot変数に格納

どちらも歯抜けの無い時間間隔が均一なデータを入力として与える必要がある。歯抜けはないものの途中でデータ点の時間間隔が変わってしまっているデータを与える場合、tdpwrspc ではnbox, nshift キーワードを明示的に指定すればその箱幅・シフトでdynamic spectrumを計算してくれる。

IDL> tdpwrspc, 'thg_mag_fsmi_x', nbox=512, nshift=256 

tdpwrspc は実際にはdpwrspc に処理を投げているので、周波数スペクトル計算の詳細についてはdpwrspc の中身を参照する。

Km -> Re の変換。座標の計算でよく使う。そのままだとxxxx_RE という新たなtplot var.を作ってしまうので、通常は”replace” keywordをつけて実行する。また”km” keywordを指定すると、逆にRekmという変換をする。またTDASに入っているtkm2re.pro そのままだと1 Re = 6374 km 換算となっているので、例えばGEOPACKで仮定しているように1 Re = 6371.2 kmで計算したい場合は、tkm2re.pro を自分のディレクトリにコピーして、中身を書き変えて使う必要がある。

Ex)  tkm2re, ‘tplt_var’, /replace

直交座標の値が入ったN×3の配列、またはtplot変数を引数に取り、球座標に変換する。

Ex)  get_data, ‘thc_fgs’, data=d    ;; d.yはN×3 の配列、磁場3成分
     xyz_to_polar, d, magnitude=r, theta=the, phi=phi 
       ;; r, the, phi はdと同様の構造体となる
     xyz_to_polar, d.y, magnitude=r, theta=the, phi=phi
      ;; N×3配列を引数にすることもできる
     xyz_to_polar, ‘thc_fgs’   ;; tplot変数も引数に取れる
->	thc_fgs_mag, thc_fgs_th, thc_fgs_ph という３つのtplot変数に結果が格納される

tplot変数ではなく、スカラーまたは配列を引数に取り、直交座標と極座標とを相互に変換(x,y,z  r, theta, phi )をする。DLMであり、かつベクトル化されているため、高速に処理される。どちらの方向に変換するかはTO_RECT, TO_SPHERE keywordで指定する。theta, phiの単位はラジアンであるが、DEGREE keywordを指定すると孤度(°)になる。

Ex) geopack_sphcar, r, theta, phi, x, y, z, /TO_RECT
geopack_sphcar, x, y, z, r, theta, phi, /TO_SPHERE

tplot変数を引数にとる万能座標変換procedure。 Ex) ace_mfi_load ;;ACE衛星の磁場データをロード cotrans, ’ace_K0_mfi_BGSEc’, ’ace_mag_Bgsm’, /GSE2GSM でGSM変換したIMF dataを生成

地球物理でよく使われる座標系の変換を行う。これもベクトル化されたDLMゆえに、非常に高速に処理される。FROM_XXX, TO_XXX (XXXは座標系名)というkeywordで変換元、変換先を指定する。変換可能な座標系は、GEO (地理座標), MAG (dipole地磁気座標), GEI (慣性座標), SM, GSM, GSE。このDLMを使って座標変換する前に、必ずgeopack_recalcに変換する値の時刻を入れて実行し、各種係数を求めておく必要がある。

Ex) 地理座標系の(x0,y0,z0) をGSM座標系に変換する例 (2000年11月30日13:20 UT)
IDL> geopack_recalc, 2000, 11, 30, 13, 20, 0, /date  ;;時刻を指定
IDL> geopack_conv_coord, x0, y0, z0, x1, y1, z1, /FROM_GEO, /TO_GSM 
	;;GSM座標系に変換された結果は(x1,y1,z1)となる 

使用できるkeywordのリスト: /FROM_GEO, /FROM_MAG, /FROM_GEI, /FROM_SM, /FROM_GSM, /FROM_GSE, /TO_GEO, /TO_MAG, /TO_GEI, /TO_SM, /TO_GSM, /TO_GSE

THEMISデータを格納したtplot変数を引数に取り、座標変換を行う。地球物理座標間の変換であればcotransでも可能だが、衛星座標系、衛星スピン座標系などTHEMIS衛星固有の座標から/への変換をする場合はこのthm_cotransを使う必要がある。

Ex) thm_cotrans, ‘tha_fgs_dsl’, in_coord=’dsl’, out_coord=’gsm’ 
      ;;Th-A衛星の磁場データをDSL座標系からGSM座標系に変換

Altitude-Adjusted GeoMagnetic Coordinates (AACGM) 座標系への変換を行うDLM。TDASにはIDL-nativeなAACGMライブラリもあるが、DLMゆえにこちらの方が数十倍高速。座標変換する前に、どの年のAACGM係数を使うかをaacgm_load_coef, yyyy (yyyyは西暦、1975,1980,1985,1990,1995,2000,2005 のどれか)で指定する。1番近い年の係数を使う。

IDL> aacgm_load_coef, 2000            ;;2000年版の係数をロード
IDL> aacgm_conv_coord, lat0, lon0, height, lat1, lon1, err, /TO_AACGM
 	 ;;高度 height [km]の地理座標の緯度・経度 (lat0,lon0)を
	 ;;AACGM座標の緯度・経度 (lat1,lon1)に変換
 ;; lat,lonの単位は[degree]なのに注意

IDL> aacgm_conv_coord, lat0, lon0, height, lat1, lon1, err, /TO_GEO
 ;;AACGM座標の緯度・経度 (lat0,lon0)を
 ;;高度height [km]での地理緯度・経度 (lat1,lon1)に変換

IDL> mlt = aacgm_mlt( 2000, yrsec, mlon )     ;;AACGM MLTを計算
 ;; yrsec: その年の1月1日00:00UTからの通し秒 (long型整数)
 ;; 	別途time_struct等を駆使して計算する必要あり
 ;; mlon: AACGM経度 [degree]
 ;; mltはdecimal hour (e.g., 23:30  23.5)

IDL> mlon= aacgm_mlong( 2000, yrsec, mlt)  
;;AACGM MLT  AACGM経度への変換
;;経度として負の値が返ることがあるので、[0,360]の範囲に直したい場合は
;; mlon = (mlon+360) MOD 360. などとする