2016年2月3日水曜日

SRTM GL1 (1秒メッシュ) の欠損値処理

他大の学生さんから解析手法について問い合わせがあったのでついでにブログ投稿。


地形データを利用する様々な分野の多くの研究者によって使われているSRTM DEMの前処理についてです。

2000年2月のスペースシャトルエンデバーにより計測されたSRTM DEM は3秒メッシュ (約90 m) 解像度 (SRTM3) については当初から全世界のDEMが公開されていましたが、1秒メッシュ(約30 m) 解像度は米国内の領域しか公開されていませんでした。

2014年には米国以外の領域の1秒メッシュ解像度DEMも公開され (SRTM GL1) ダウンロードができるようになっていますが、利用者の利便性を考えてなのか、本来計測できていない部分の標高値も他のDEM (ASTER GDEM、GMTED、NED) を使って補間してしまっています。

手っ取り早く高解像度のDEMを使いたい人にとっては問題ないのですが、地形や氷河体積の時間変化などを明らかにしたい人にとっては、2000年時点 (スナップショット) での標高値ではない標高 (例えばASTER GDEM はグリッドによって2000年~2011年と異なる) が混ざるため困ります。

そこで標高の時間変化計算に用いるためにはSRTM GL1から外来DEMのグリッドを除去する処理が不可欠となります。SRTM GL1には対応するタイルごとにSRTMGL1Nというフラグデータがあり、SRTM GL1のグリッドごとに、外来DEMを使ったかどうか使った場合はどの種類かがわかる情報があります。この情報を元にオリジナル のSRTMを使っているグリッドのみを抽出すればOKです。

QGISでの読み込み

以上の処理はRやMATLABなどのプログラミングでももちろんできますが、ここではQGISを使ってやり方を紹介します。

SRTM GL1Nの仕様のページを見ると、値は8 bit (0-255) で格納されており、201以上の値を使えば大丈夫そうです。SRTM GL1のフォーマット (*.hgt) は単純な raw binary形式ですが、QGISはSRTMのファイル名(ex. NaaEbbb.hgtの場合は北緯aa度東経bbb度を北西端)から自動で認識して正しい位置に表示してくれます。

SRTM GL1Nも単純な raw binary 形式となっていますがこちらはQGISが自動認識はしてくれないためヘッダーファイル (*.hdr) を用意して位置情報などのデータに関する上をQGISに教えてやる必要があります。


hdrファイルとは

テキスト形式でデータについての情報を記述したものです。 raw binaryファイルと拡張子より前の部分の名前を揃える必要があります (ex. データがN27E086.numの場合はN27E086.hdr)。

詳しくはこちらのページのEHdr -- ESRI .hdr labelled を参照
http://www.gdal.org/frmt_various.html


N27E086.num の場合は、N27E086.hdrという名前で以下の内容が記述してあるテキストファイルを用意すればOKです。


BYTEORDER      I
LAYOUT         BIL
NROWS          3601
NCOLS          3601
NBANDS         1
NBITS          8
BANDROWBYTES   7202
TOTALROWBYTES  7202
PIXELTYPE      UNSIGNEDINT
ULXMAP         86
ULYMAP         28
XDIM           0.000277777777777778
YDIM           0.000277777777777778
NODATA         0


一応こちらからダウンロードできます。適当なエディタで開いてみてください。
https://dl.dropboxusercontent.com/u/870568/blogger/N27E086.hdr


準備ができたら N27E086.num をQGISにドラッグ&ドロップすれば表示されるはずです。あとはQGISのラスターメニューの変換コマンド (または右クリックで別名保存) でGeoTIFF形式に変換すると他のソフトウェアでも簡単に読み込むことができます。
SRTM GL1のファイル名と同じだとややこしいので N27E086_num.tifとして保存します。


SRTM GL1からオリジナルのSRTMのみの抽出

QGISにSRTM GL1 (N27E086.hgt) とSRTM GL1N (N27E086_num.tif) の2つを読み込んで置いて下さい。

まずラスタ計算機 (Raster calculater) で以下の内容を図のようにウィンドウ下部の演算欄に入力します。

"N27E086@1" * ("N27E086_num@1" >= 201)



フラグデータで201以上の値が入っている場所のSRTM GL1のみを抽出するという意味です。ウィンドウ右上の出力ファイル名も指定しておいて下さい。
出力ファイル名を N27E086_temp.tif とします。

変換して出力されたファイルは以下の図のようにSRTMじゃない部分 (つまり欠損値) に数値の0が入っています。



標高変化の計算のためには0じゃなくてNo Dataだと認識されて欲しいのでもうひと手間かけます。

N27E086_temp.tif を右クリックで別名保存で、ウィンドウ下部のNo data valueで0を指定します。出力ファイル名を指定してOKをクリックします。
ここでは出力ファイル名をN27E086_void.tif とします。



以下のような感じでできあがりです。

0 件のコメント:

コメントを投稿