トップ 一覧 検索 ヘルプ RSS ログイン

JAXA-MODISの変更点

  • 追加された行はこのように表示されます。
  • 削除された行はこのように表示されます。
!!! JAXA recieved MODIS (level1b, cha, SST products)
CEReS では 宇宙航空研究開発機構 (JAXA) および東海大学(TSIC/TRIC) で受信,処理された Terra/Aqua 搭載の MODIS (MODerate resolution Imaging Spectroratiometer; 中分解能撮像分光放射計) の取得および内部公開(現状では千葉大学西千葉キャンパスのみ)を行っています.

!! 公開 ftp site 
*ftp://modis.chiba-u.ac.jp/pub/jaxa-modis/ (anonymous ftp; 千葉大西千葉キャンパス内であれば取得可能)
*ftp://modis.cr.chiba-u.ac.jp/pub/jaxa-modis/ (anonymous ftp; 千葉大西千葉キャンパス内であれば取得可能)
** sub-directory で,disk1/year/year(4文字)month(2文字)day(2文字)となっています.このルールから該当する(欲しい)データの日付が分かると思います.
** ftp server の中身を見れば分かりますが,データは2004年8月01日〜現在まであります.
** global データセットではありません.日本付近です.CEReS受信の AVHRRよりもやや狭い程度の受信状況となっています.
!! ファイル名ルールおよびフォーマット仕様
!JAXA/TSIC 提供プロダクトファイル名(クロロフィルa濃度,海面温度 [SST])およびファイル仕様
* A2GL***********_****_pixelサイズ_lineサイズ_chla or sst (クロロフィルかSSTかを示す).gz 
 ***で書かれているところは気にしないで,ファイル名の末尾(*.gz は gzip で圧縮されていることを示します)に着目してデータ選択をしてください.ファイル名中,pixelサイズとlineサイズは5桁の整数で示されています.

* ファイル仕様
 バイナリデータ仕様は以下のような1line 分のヘッダと2byte 整数(符号なし整数,bigendian)で構成されています.

,ヘッダ(pixel数 分のバイト,テキストファイル), データ本体(line数 x pixel 数 x 2byte)

*ヘッダ部分(テキストファイル)
 pixel数 ライン数 左上緯度 左上経度 緯度経度刻み(空間解像度) 
 物理量換算傾き(slope) 物理量換算切片(offset) パラメータ名 入力ファイル名
 残りは空白あるいは使用する必要はありません.各ヘッダ要素はカンマ(,)で区切られています.UNIXコマンドの''head'' と ''awk'' 等を組み合わせることで必要な情報を得ることが可能だと思います.

*物理量換算
以下の式で換算することができます.
**Chl-a [mg/m^3] = (digital number x slope) + offset
**SST [K] = (digital number x slope) + offset

! JAXA プロダクト(chla, SST)に関しての情報(リンク)
* [JAXA MODIS product に関するFAQ|http://kuroshio.eorc.jaxa.jp/ADEOS/mod_nrt_new/html/05_faq.html]
* [JAXA MODIS product バイナリフォーマット説明|http://kuroshio.eorc.jaxa.jp/ADEOS/mod_nrt_new/PDF/MODISform.pdf]
!! JAXA 受信 MODIS level 1b データに関する情報
! JAXA 受信 MODIS ファイル名ルール
JAXAで受信してるMODIS hdf file (これはレシーバの仕様なんだと思います)はNASAが公開しているタイル管理のファイル名ルールとは異なりますので注意が必要です.多分こういうルールなんだと勝手に解釈しました.
# 違っていたら誰か教えてください.

 MOD ... Terra 衛星(日中午前overpass), MYD ... Aqua 衛星 (午後に overpass)
 02 or 03 ... 02:画像データ,03:軌道関連情報(幾何補正処理で必要になります).
 1KM, HKM, QKM:
 1KM ...1km 分解能のデータが入る
 (注:hdfの中を見ると,500m, 250m のchannelデータも入っていますが,全て1kmに平均化)
 HKM: 500m 分解能,QKM: 250m 分解能.
 JYYYYMMDDHHMMSS.YYYYMMDDHHMMSS.hdf.gz <- gzip で圧縮されている.
 最初は受信開始時刻情報,次が受信終了時刻情報.
 YYYY... 年(4桁.西暦),MM ... 月(2桁.1-9月でも1桁に詰めません),
 DD ... 日(2桁.月と同じく詰めません),HH...時,MM ...分,SS...秒 です.
 # J は Japan のJ なんでしょう :-)
 時刻は Z (UTC) です.JST ではありません.

例: MYD021KM.J20050406033924.20050406035008.hdf.gz
↑のルールに従えば,Aqua衛星, 画像データ,1km分解能,2005年04月06日03:39:24-2005年04月06日03:50:08 UTCの受信 hdf ファイルということになります.

! 必要とする tool のインストール
Linux での使用を前提に記載しています.それぞれのサイトに Windows や Mac,UNIX の version のものもあるかもしれません.適宜読み替えて下さい.
*[Java|http://www.java.com/ja/] MODISSwathTool の GUI が java で書かれているため必要です.
Linux の場合,rpm を入れることになるので,Linux distribution に依存しますが,install後,path を通して下さい.
  例)
  export JRE_HOME=/usr/java/default
  export PATH=/usr/java/default/bin:$PATH
  export LD_LIBRARY_PATH=/usr/java/default/lib:$LD_LIBRARY_PATH

*[MODIS Reprojection Tool - Swath (MODISSwathTool)|https://lpdaac.usgs.gov/lpdaac/tools] ユーザ登録が必要ですが,Free で利用可能です.
対話式のインストーラであるため,install に実行を与え(chmod +x install),管理者で実行.ここではインストール先を /usr/local/MRTSwath に設定し,path を通してください.
 例)
 export PATH=/usr/local/MRTSwath/bin:$PATH
  export MRTDATADIR=/usr/local/MRTSwath/data:$MRTDATADIR
*[HDF4|http://www.hdfgroup.org/] HDF4 (HDF5 と HDF4 は互換性がない?.EOS-HDFはHDF4をベースとしているため必要)
Linux distribution によっては足りないパッケージがありますので(RH系だと flex, bison を yum を使って事前に入れておく必要がある.gcc 等の開発環境は無論必要 ;-),さらに libjpeg や libjpeg-devel も必要),入れておいてください.
  pc{usr} %./configure --prefix=/usr/local/hdf4.2r3 <- HDF4.2r3 の場合.適宜読み替える
  $ make 
  # mkdir /usr/local/hdf4.2r3
  # make install

その後同じく path を通しておいてください.
  例
  export PATH=/usr/local/hdf4.2r3/bin:$PATH
  export LD_LIBRARY_PATH=/usr/local/hdf4.2r3/lib:$LD_LIBRARY_PATH
  export HDFINC=/usr/local/hdf4.2r3/include
! ModisSwathTool の使い方
インストールで示したソフトウエアが無事インストールされており,path が通っていれば,↓で示すようなGUI画面が,

pc{usr}$ ModisSwathTool 

で表れます.直感的に分かり易いGUIなので分かるかと思います.

{{ref_image ModisSwathTool.png}}

*'''Input File [Open Input File]:''' 対象となる画像データを開きます.開く前には *.gzファイルを解凍する必要があります(以下同じです).
*'''Input File Info:''' hdf ファイルとして開ければ,ファイルの内部情報がここに出ます.channel情報,等々です.
*'''Available Bands, Selected Bands:''' 前者はhdf file内に含まれる channel情報,後者は変換をかけるtarget の channelを示します.>> or << ボタンで,どのチャンネルを変換するかを選択・削除ができます.
*'''UL Corner, LR Corner:''' 緯度経度直交座標系に変換する際に左上の緯度・経度,右下の緯度・経度を入力します.分,秒表記ではありませんので注意.
*'''Spatial Subset:''' 緯度経度直交座標(LAT_LONG), 選択図法の形式(PROJ_COORDS),単純なラインーピクセル(LINE_SAMPLE)の3種のどれかを選択します.ここでは他のデータとの整合性を取るために緯度経度直交座標系を選びます.
*'''Geolocation  File:''' ファイル名ルールで記載した軌道情報 hdf ファイル(MOD03, MYD03) を選択します.これがないと変換できないようです.Input file と同じ受信 hdf を選ぶ必要があります. 
*'''Output File:''' 出力先.ファイルが無いと選択できないようなので,適宜,$ touch temp とかしてファイルを用意しておいた方が良いです.
*'''Output File Type:''' hdf 形式,geotiff 形式,raw bainary 形式の3つのうちどれかを選択.raw binaryを選択しました.hdf をサポートしている画像処理ソフトをよく使う場合,hdfの方が良いかもしれません. 
*'''Resampling Type:''' 元データをどの方法でサンプリングし直すかの設定.NN (Nearest Neighbor)とBL(Bi-Linear)とCC(Cubic Convolution)の3種類あります.元データと同程度の空間分解能を出す場合はNN(元データは消えない),より荒いデータに変換する場合には BLかCCという選択が良いと思います.NN, BL, CCの違いはリモセンの教科書に出ていますが,簡単に書くと;BB:一番近い場所のpixelを選ぶ,BL:周囲4点のpixelsを線形内挿.CC:周囲16点のpixels を距離の重み付けをして内挿.BB以外は元データが消えます.
*'''Output Projection Type: '''いっぱい選択肢が出てきますが,上のSpatial Subset で"PROJ_COORDS"を選んでいない場合は,GEOが無難でしょう.
*'''Output Data Type:''' raw binary と hdf だけ必要なoption? 出力 data 関数系を選びます.1byte 整数(符号付き[CHAR8 or INT8], 符号無し[UINT8]),2byte 整数(同じく [INT16, UINT16])を選びます.選ばなかった場合,元データと同じ(INT16)になります.
*'''Output Pixel Size:''' 緯度経度直交系を選ぶと degree (°)になります.非常に大雑把ですが,0.01 ...1km, 0.005 ... 500m, 0.0025 ... 250m と考えると楽です(厳密には緯度によって° <-> 距離の換算が変わりますので正しくないです.念のため)
'''Commands'''
*'''Load Parameter:''' 以前設定したパラメータファイルを呼び出します.
*'''Save Parameter: ここで設定した上記パラメータを保存します.後述するコマンド,一括処理時に重要となりますので,saveした方が無難です.

ここまで一通り設定をした後でRunボタンを押せば,設定項目が間違っていなければ変換が開始されます.output file で "hdf" を選択した場合は一つのhdfファイルが,"raw binary" を選んだ場合は,

[Output file]+[_選択されたchannel名].hdr (ヘッダファイル,テキスト),.dat (選択した data type でのバイナリファイル)

が選択したチャンネル数だけ出力されます.ここまでいけば成功で,成功したパラメータファイルを保存しておきます(後で使います).これが GUI な ModisSwathToolの使い方です.また,NASA で公開している各種MODIS product も"ModisTool"でやり方は本質的には同じです.

! CUI での処理(一括処理)
大量に同じ領域の切り出しをするのにいちいちModisSwathToolを立ち上げてマウスをクリックしてファイルを選んで…とやっていると日が暮れてしまいますので,一括で処理する方法(このやり方が正しいかどうかは不明です)を書いておきます.

ModisSwathTool で実際に処理に回すときには,swath2grid を呼び出して,パラメータを流し込んでいるだけです.言い換えると,ModisSwathTool は swath2grid のフロントエンドということになります.つまり,swath2grid のやり方を覚えてしまえば,いちいちModisSwathToolを使う必要はないということです.やりかたは,

$swath2grid -help

と打てば方法が出てきます.

この optionを全部いちいち打って手作業で処理をするのもばかげていますので,ModisSwathToolで作成したパラメータファイルを使います.パラメータファイル自身は単なるテキストファイルです.

一括処理の前手順としては
*パラメータファイル中の "INPUT_FILENAME"(画像 hdfファイル) , "GEOLOCATION FILE"(軌道情報hdfファイル),"OUTPUT FILENAME" (出力先ファイル) の行をコメントアウト(頭に # を付ける)する(理由: swath2grid の引数として入れ込むため,ここでファイル名が指定されていると読めなくなる).
*↑それだけです.後は他の項目がそのままで良いかを確認します.

逆に言えば,パラメータファイルのひな形があり,内部構成さえ分かれば ModisSwathToolは不要です.これでようやく手を抜けるわけです.original ではいっぱい引数を入れる必要がありますが,パラメータファイルを使えば,

 $swath2grid -if=${input_image_hdf_file} -gf=${orbit_hdf_file} 
    -of=${output_file} -pf=${parameter_file} (実際は 1行)

で一括変換ができます.$以下は適宜shell で宣言するようにすれば良いでしょう.

 (例: bashでの記載)
 #! /bin/sh
 temppath=/work/hogehoge/temp/
 outpath=/work/hogehoge/output-bins/
 prmpath=/work/hogehoge/convert-prms/
 prmconf=typical.prm
 # 元々のgzip をzcat で解凍させながら,temppath へ
 # gzを取ったファイルとして吐き出します.
 # gzip でその場出だしても良いですが,複数人で使うことを
 # 常に想定して(+nfs マウントで見ることも)zcat を使います.
 # extract orbit file hdf (MOD03.*)
        for obtf in MOD03.* ;
        do zcat ${obtf} > ${temppath}${obtf%.gz};
        done
 # extreact 500m res. Terra hdf (MOD02HKM.*)
        for hdf in MOD02HKM.*;
        do zcat ${hdf} > ${temppath}${hdf%.gz};
        done
        cd ${temppath}
        for fname in MOD02HKM.*.hdf;
 # ↓の部分が上で書いた部分を実際にさせたもの.軌道hdf と画像hdf は
 # ファイル名の頭の部分のみが違うので,整合性を取るようにしている.
        do swath2grid.comp -if=${fname} -gf='MOD03.'${fname#MOD02HKM.} \
                -of=${outpath}${fname%.hdf} -pf=${prmpath}${prmconf}
 
このシェルでは,同じdirectoryにgz なhdfが皆あれば,temppathに typical.prmで指定したoptionに従い,outpathで指定したdirectoryにMOD2HKM…_prmファイルで指定したチャンネル名.hdf 及び .datのヘッダ,バイナリを吐き出すようになります(prmファイルでoutputがraw binaryになっていた場合).

この方法が使えるようになれば,後は得られたバイナリファイルをチャンネル毎の単純なバイナリデータとして取り扱うことができます.

! 物理量換算はどうするのか?(Radiometric 補正)
hdf ファイルを生成する際に CCT値 (DN)から SI (Scaled Integer)に変更された2バイト整数です.この値をradiance やreflectance に変換するためには,hdf ファイルからヘッダ情報を読み取り,必要な係数を取り出す必要があります(そのために HDF4 が必要です).

$(HDF4のバイナリのpath)ncdump -h hogehoge.hdf

この中の変換係数(上記例の一部を抜粋)
 EV_500_RefSB:radiance_scales = 0.021505227f, 0.018739114f, 0.00529294f, 0.0025430776f, 0.00078053412f ;
 EV_500_RefSB:radiance_offsets = 316.9722f, 316.9722f, 316.9722f, 316.9722f, 316.9722f ;
 EV_500_RefSB:radiance_units = "Watts/m^2/micrometer/steradian" ;
 EV_500_RefSB:reflectance_scales = 3.2494634e-05f, 3.1683794e-05f, 3.5203753e-05f, 3.3397813e-05f, 2.7262724e-05f ;
 EV_500_RefSB:reflectance_offsets = 316.9722f, 316.9722f, 316.9722f, 316.9722f, 316.9722f ;
 EV_500_RefSB:reflectance_units = "none" ;
 EV_500_RefSB:corrected_counts_scales = 0.12619403f, 0.12619403f, 0.12619403f, 0.12619403f, 0.12619403f ;
 EV_500_RefSB:corrected_counts_offsets = 316.9722f, 316.9722f, 316.9722f, 316.9722f, 316.9722f ;
 EV_500_RefSB:corrected_counts_units = "counts" ;

SIから分光反射率に換算するためには;
*Reflectance (無次元) =  reflectance_scales * (SI - reflectance_offsets)
とし,分光放射は;
*radiance (W m^-2 microm^-1 sr^-1) = radiance_scales * (SI - radiance_offsets)
補正済みカウント値への変換は;
*DN** (counts) = corrected_counts_scales * (SI - corrected_counts_offsets)

となります.DN** は第2補正済カウント値ということらしい,です.式自身は variable = scale * (SI - offset) で統一されています.係数の取り出しはいろいろな方法がありますが,例として grep とawkの組み合わせを下に書いておきます.

 例(CCT 〜 放射量への変換係数をheaderから取り出す[500mの例])
 $ grep EV_500_RefSB:radiance_scales < MOD02HKM.J20050411004553.20050411005445.hdr | \
 awk  '{print $3, $4, $5, $6, $7}' | \
 awk -F, '{printf"%1.10f %1.10f %1.10f %1.10f %1.10f\n", $1, $2, $3, $4, $5}' > coeff.txt

とすると,

0.0215052270 0.0187391140 0.0052929400 0.0025430776 0.0007805341

と数字だけ取り出せます.後はこれらの係数をバイナリに掛け合わせるプログラムを作成すれば良いだけです.
CEReS Database Wiki