GDALでGCOM-Cの可視領域のデータからTrue Colorイメージを再構成してみた

この記事は、FOSS4G Advent Calendar 2019 の12月7日の記事です。 昨日の記事は、 northprintさんの「mapboxのMaps SDK for Androidでお好みのラスタータイルを表示する」 でした。

TL;DR

今回は2017年の12月に打ち上げられた、気候変動観測衛星「しきさい」(GCOM-C)の可視領域の数値データを、 GDAL2.4で再構成することで、True Color Imageに再構成してみました。

結果は、こんな感じになりました。欠損値は透明にしてあります。

(Original data provided by JAXA)

あと再構成プログラムをC++で書きましたが、Pythonのほうがよかったかも…

「しきさい」の概要

気候変動観測衛星「しきさい」は、地球環境変動観測ミッション(GCOM)の陸域版の衛星です。 GCOMは地球から放射される様々な光を長期間観測することで、地球の気候変動の監視と メカニズムの解明を目的としたプロジェクトです。水循環を観測するGCOM-Wと、気候変動を観測する GCOM-Cという二つの衛星で観測を行っています。 GCOM-Cに搭載されている、SGLIという多波長光学放射計は、可視、近赤外、近紫外領域の 19の観測波長帯を持っています。

今回は可視領域の観測データを使って目視に近い画像である、True Colorイメージを作成してみます。

データの入手方法

まずは、GCOM-Cの観測データを入手しましょう。 JaxaのG-Portal(https://gportal.jaxa.jp/gpr/) というページでユーザー登録をすると、FTPで データをダウンロードすることができます。また、G-PortalのWebサイトにログインすると データの検索やダウンロード、フォーマット変換ができるようです。 私はデータ変換を自動化する必要がある別の目的があるので、あえて自分でやることにします。 詳しくは 「G-Portal 地球観測衛星データ提供システムユーザ向け取扱説明書」と 「気候変動観測衛星「しきさい」(GCOM-C) データ利用ハンドブック」を参照してください。

それでは、FTPサーバにログインして、/standard/GCOM-C/GCOM-C.SGLI 以下にある、L3.LAND.RV03、L3.LAND.RV05、L3.LAND.RV07 から同じ日付のデータを ダウンロードします。

今回は2019年11月の1ヶ月間に観測された平均値を使用しました。ダウンロードしたファイルは以下のとおりです。

  • GC1SG1_20191101D01M_D0000_3MSG_RV03F_1001.h5
  • GC1SG1_20191101D01M_D0000_3MSG_RV05F_1001.h5
  • GC1SG1_20191101D01M_D0000_3MSG_RV07F_1001.h5

ファイル名の後ろのほうにある、RV03、RV05、RV07はプロダクトIDとよばれる観測から得られた物理量の種類を 示すIDです。RVxxのIDを持つプロダクトは陸域の大気補正済みの反射率で、RV03が青、RV06が緑、RV07が赤の光に 相当します。

各データに関しては 「SGLI高次プロダクト フォーマット説明書」 を参照してください。

データの再構成

では、ダウンロードしたデータを再構成していきます。ダウンロードしたファイルはHDF5という形式で保存されています。 HDF5はデータやパラメータを階層的に保存できる形式で、観測値だけでなく観測した時間や、処理したアルゴリズム、パラメータ等も いっしょに一つのファイルに格納できるすぐれものです。このHDF5形式のファイルからデータを取り出して、GDALを使って GTiff形式に変換します。

可視画像は、R、G、Bの3つピクセルから構成されています。GTiffでは、それぞれ別のバンドに格納すれば良いので、 取り出したデータを各バンドにコピーするだけでOKです。簡単ですね。

以下のプログラムは、HDF5とGDALのライブラリを使用しています。DebianやUbuntuを使用している方は、

$ sudo apt install libgdal-dev libhdf5-dev

で必要なライブラリをインストールすることができます。

C++なので、あまりちょっととっつきにくいかもしれませんが、流れはコメントに書いたとおりです。 同じ流れでPythonで書くことができるはずです。

int main(void)
{
	// GDALの初期化
	GDALAllRegister();
	
	// Dataを格納する配列
	uint16_t data[WIDTH][HEIGHT];

	// File名が長いので、定数化
	static const string RED   = "GC1SG1_20191101D01M_D0000_3MSG_RV03F_1001.h5";
	static const string GREEN = "GC1SG1_20191101D01M_D0000_3MSG_RV05F_1001.h5";
	static const string BLUE  = "GC1SG1_20191101D01M_D0000_3MSG_RV07F_1001.h5";


	// GDALのGTiffドライバを呼び出す
	GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	// 呼び出したドライバで、BANDを4つもったデータセットを作成
	GDALDataset *poDstDS = poDriver->Create( "out.tiff", WIDTH, HEIGHT, 4, GDT_Byte, nullptr);

	// 1番目のバンドを取り出して
	band = poDstDS->GetRasterBand(1);
	// このバンドを赤に指定する
	band->SetColorInterpretation(GCI_RedBand);
	// HDF5からデータを読んでから
	readData(RED, data, "/Image_data/Rs_VN07_AVE");
	// 取り出したバンドにコピーする。
	copy2band(data, band);

	// 2番目は緑。赤と同様
	GDALRasterBand *band = poDstDS->GetRasterBand(2);
	band->SetColorInterpretation(GCI_GreenBand);
	readData(GREEN, data, "/Image_data/Rs_VN05_AVE");
	copy2band(data, band);

	// 3番目は青。赤と同様
	band = poDstDS->GetRasterBand(3);
	band->SetColorInterpretation(GCI_BlueBand);
	readData(BLUE, data, "/Image_data/Rs_VN03_AVE");
	copy2band(data, band);

	// 4番目はアルファチャンネル。欠損値を透明にする
	band = poDstDS->GetRasterBand(4);
	band->SetColorInterpretation(GCI_AlphaBand);
	genAlphaBand(data, band, WIDTH, HEIGHT);

	// 最後にデータセットをクローズしてできあがり
	GDALClose(poDstDS);
    return 0;

HDFからデータを取り出すreadDataは、 ただデータを取り出して、引数に指定されたアドレスにコピーしているだけです。が、 使用しているHDFライブラリでは、DataSpaceにデータを格納するのがお作法のようです。 コピーするためには、コピー元とコピー先の範囲を指定必要があります。 変数spとmemspaceのselectAll関数を呼ぶことで、DataSpace全体を指定しています。

void readData(std::string filename, void *data, std::string dataset)
{
    H5File file(filename.c_str(), H5F_ACC_RDONLY);
    DataSet ds = grp.openDataSet(dataset.c_str());

    DataSpace sp = ds.getSpace();
    const int rank = sp.getSimpleExtentNdims();
    hsize_t count[] = {WIDTH, HEIGHT, 1};
    hsize_t offset[] = {0, 0, 0};
    sp.selectAll();

    DataSpace memspace(rank, count);
    memspace.selectAll();
    ds.read(data, PredType::NATIVE_UINT16, memspace, sp);
}

取り出したデータをbandにコピーするcopy2bandは、以下のとおりです。 データをスキャンして、欠損値を意味する65535を0にして、 14ビットのデータを6bit右シフトして8bitに丸めてたあと、bandにコピーしています。

void copy2band(void *data, GDALRasterBand *band)
{
	// dataを2次元配列としてアクセスしたいので、memcpyで配列にコピー
	uint16_t val[WIDTH][HEIGHT];
	memcpy(val, data, sizeof(uint16_t) * width * height);
    
	uint8_t pixel[WIDTH][HEIGHT];

	// すべてのデータをスキャンしてデータサイズを8ビットに変換
	for(int y = 0; y < HEIGHT; y++)
		for (int x = 0; x < width; x++)
			pixel[x][y] = (val[x][y] == 65535) ? 0 : (val[x][y] >> 6);
	// バンドにコピー
	band->RasterIO(GF_Write, 0, 0, WIDTH, HEIGHT, pixel, WIDTH, HEIGHT, GDT_Byte, 0, 0);
}

最後にアルファチャンネルを生成しているgenAlphaBand関数はこんな感じにしました。

void genAlphaBand(void *data, GDALRasterBand *band)
{
        uint16_t val[WIDTH][HEIGHT];
        memcpy(val, data, sizeof(uint16_t) * WIDTH * HEIGHT);
        uint8_t tmp[WIDTH][HEIGHT];
        for(int y = 0; y < HEIGHT; y++)
                for (int x = 0; x < WIDTH; x++){
                        tmp[x][y] = (val[x][y] == 65535) ? 0 : 255;
                }
        band->RasterIO(GF_Write, 0, 0, WIDTH, HEIGHT, tmp, WIDTH, HEIGHT, GDT_Byte, 0, 0);
}

ほぼcopy2bandと同じです。違うのは欠損値だった場合は値を0、そうでないときは255を使用しているだけです。 ここもそうですが、全体的にリファクタリングの余地がありますね。

結果はTL;DRのところに載せたとおりです。 実行する時は、確実にスタック領域が不足するので、ulimitコマンドでスタック領域を無制限にすることを忘れずに。 さもないとSegmentation Faultで落ちます。

$ ulimit -s unlimited

今回は、GCOM-Cの可視光のデータを使ってTrue Color Imageを再構成してみました。 GCOM-Cは可視光以外にも、海水面温度や地表面温度、海の色、エアロゾル、雲頂高度といった興味深い値を観測し、 その数値データを数日遅れで入手することができます。

興味がある方は、ぜひ他のデータも可視化してみてください。

明日は、Kanahiroさんです。