GSMaPのバイナリーデータの取得から読み取りまで

仕事の関係上、GSMaP（衛星全球降水マップ）を使う必要性があったので使ってみました。GSMaPはテキストデータとバイナリデータで提供されていますが、ある程度まとまったデータを取り扱おうとするとバイナリーデータの方が便利です。ここでは、GSMaPのデータファイル取得から数値データの読み取り方法まで一通りの手順を解説します。

1. はじめに

仕事でRRIモデルを使って過去に起こった洪水・氾濫の解析をすることとなりましたが、対象地域は日本ではなくてインドネシアでした。当然、当地ではなかなか当てになる降雨データが入手できませんでした。そこで、人工衛星データをベースとして全世界的な降雨分布を提供しているGSMaPを使うことになりました。

この記事ではGSMaPのバイナリーデータを取得し、それらから数値データを読み取る方法について解説していきたいと思います。

2. GSMaPとは

GSMaP（Global Satellite Mapping of Percipitation）「衛星全球降水マップ」は複数の降水を観測する人工衛星や静止気象衛星を組み合わせて開発した世界の雨水マップです。

データには全世界版のバイナリデータと地域版（15の地域に分割された地域のデータ）のテキストデータがあります。

バイナリデータは、時間データ（hourly）、ゲージ補正済みの時間データ（hourly_G）、日データ（daily）、ゲージ補正済みの日データ（daily_G）、月データ（monthly）、ゲージ補正済みの月データ（monthly_G）の6種類で1998年ころから揃っています。

https://sharaku.eorc.jaxa.jp/GSMaP/guide_j.html#
より参照

3. GSMaPのバイナリデータファイルの取得

バイナリデータはWebのGUI（グラフィカルユーザインターフェース）からダウンロードすることもできますが、まとまった数量のデータを取得する場合はFTPでダウンロードするのが便利です。

FTPでGSMaPのデータをダウンロードするためには、下記URLでGSMaPへのユーザ登録をする必要があります。ユーザ登録の手続きが完了するとFTPサイトのURL、ユーザIDおよびパスワードを記載したメールが送られてきます。
https://sharaku.eorc.jaxa.jp/GSMaP/registration_j.html

下記のURLにFTPクライアントアプリ（ここではWinSPC）でアクセスします。
hokusai.eorc.jaxa.jp

FTPサイトのルートの階層には、standard、riken_nowcast、realtime_ver、realtime、now、climateというディレクトリがあります。この中で過去の補正されたデータが格納されているのはstandardですので、その中のデータを選択します。

standardのディレクトリの下にはv8からv5までに分かれています。特に理由がなければ最も新しいv8を用いればいいと思います。

v8ディレクトリの下にはhourly_G、hourly、dailyという名前のディレクトリがありますが、それぞれゲージ補正済み時間雨量、時間雨量、日雨量のデータが格納されています。txtというディレクトリには地域ごとのテキストデータが格納されています。

洪水・氾濫解析ではなるべく時間解像度が高いデータが望ましく、かつゲージ補正された時間データであるhourly_Gを選択します。

hourly_Gディレクトリの下には、さらに年／月／日のディレクトリの階層があり、日のディレクトリの下には時間ごとのデータが格納されています。個々のデータはGzip圧縮されています。これらのデータの中から解析期間に合致したものをダウンロードしてやります。

4. バイナリデータからの数値データの読み取り

解凍した個々のデータはダイレクトアクセスバイナリなので以前紹介したJavaのプログラムを改良したものが使えます。

GSMapのデータは、下図のように東西方向に3600点、南北方向に1200点からなるデータです。それぞれの点には単精度浮動小数点数の降雨データ（ここでは時間降雨）が格納されています。データの並び順は、東から西（左から右）、北から南（上から下）に一直線に並んでいます。データを読み込むときは、データを先頭から終末まで順番に読んでいきます。

ここでは、バイナリデータをESRI-ASCIIフォーマットに変換することにしていますので、読み取ったデータを2次元に並べ替えてやります。

以前の記事（https://hydrocirculate.net/read_binary/）で解説したダイレクトバイナリファイルのJavaの読み取りプログラムをGSMaP用に改良したものを載せておきます。

package net.gofort;
import java.io.*;
import java.nio.ByteBuffer;
import java.nio.ByteOrder;
public class Read_GSMaP {
    public static void main(String[] args) {

        int NCOLS = 3600;
        int NROWS = 1200;
        String inFile = "gsmap_gauge.20141217.0000.v7.0001.0.dat";
        String outFile = "test_gsmap_out.asc";
        double DX = 0.1;
        double DY = 0.1;
        double XLLCENTER = 0.05;
        double YLLCENTER = -59.95;
        try {
            DataInputStream dis = new DataInputStream(new BufferedInputStream(new FileInputStream(inFile)));
            ByteBuffer buf = ByteBuffer.allocate(4);
            float[][] meshVal = new float[NCOLS][NROWS];
            for (int j = 0; j < NROWS; j++) {
                for (int i = 0; i < NCOLS; i++) {
                    buf.order(ByteOrder.LITTLE_ENDIAN);
                    buf.putFloat(dis.readFloat());
                    buf.order(ByteOrder.BIG_ENDIAN);
                    meshVal[i][j] = buf.getFloat(0);
                    buf.clear();
                }
            }
            PrintWriter pw = new PrintWriter(new BufferedWriter(new FileWriter(outFile)));
            pw.print("NCOLS   ");
            pw.println(Integer.toString(NCOLS));
            pw.print("NROWS   ");
            pw.println(Integer.toString(NROWS));
            pw.print("XLLCENTER   ");
            pw.println(Double.toString(XLLCENTER));
            pw.print("YLLCENTER   ");
            pw.println(Double.toString(YLLCENTER));
            pw.print("DX   ");
            pw.println(Double.toString(DX));
            pw.print("DY   ");
            pw.println(Double.toString(DY));
            pw.println("NODATA_VALUE  9.999E+20");
            for (int j = 0; j < NROWS; j++) {
                for (int i = 0; i < NCOLS; i++) {
                    pw.print(String.format("%6.2f  ", meshVal[i][j]));
                }
                pw.println();
            }
            pw.flush();
            pw.close();
            dis.close();
        } catch (Exception ex) {
            ex.printStackTrace();
        }
    }
}

改変したところは、

• 8行目、9行目のそれぞれ東西方向と南北方向のグリッドの数をGSMaPのデータに合わせて変更
• つづく、10行目と11行目の入力ファイルと出力ファイル名をそれぞれ変更
• 12行目から15行目のグリッドサイドおよび南西端の座標を変更
• 43行目の南北方向の数値の読み取り方向を北から南に変更

です。

5. サンプルプログラムの使い方

このプログラムはInteliJという統合開発環境を用いて開発しているので、IntelliJから実行することができます。IntelliJから実行するためには、read_GSMaP.imlというファイルをintelliJから開いてIntelliJのプロジェクトとしてこのプログラムを読み込んでやります。IntelliJのプロジェクトとして読み込むことができれば、IntelliJ上で実行することができます。

IntelliJで実行するとトップディレクトリにtest_gsmap_out.ascというテキストファイルが出力されます。

IntelliJを使わない場合は、Javaで直接実行することができます。.Gzipファイルを解凍したトップディレクトリで、以下のコマンドでコンパイルします。
javac ./src/net/gofort/Read_GSMaP.java

つづいて、以下のコマンドで実行します。
java -cp ./src net/gofort/Read_GSMaP

プログラムを実行するとカレントディレクトリに下記のテキストファイルができます。
test_gsmap_out.asc

このファイルをGISアプリで表示すると下図のようになります。GSMaPのデータの範囲は北緯60度から南緯60度までとなっているので、地図の上下の部分にデータが表示されていないところが確認できます。

6. まとめ

この記事では、GSMaPのバイナリデータの取得から数値を読み取る方法までを解説しました。ポイントをまとめると下記のようになります。

以上、後まで読んでいただきありがとうございました。