GMT 5 Default formatをQGISで表示する

今日は、GMT 5のDefault formatでQGISで読む方法について記載します。

ここでは、GMT 5 Default format=GMTv5-netCDFとして表記します。

環境の説明

利用するQGISは、CentOS 7＋epelを使用します。QGISのバージョンは2.8.2-Wienです。

以下の理由で利用します。

QGISは、netcdfファイルを読み込むのにgdalライブラリを利用してしています。
GMT 5のデフォルトnetcdfの読み込みに、gdalがサポートしていないため、読み込めません。
Windows版のQGISでは、gdalのライブラリを加工するのに手間がかかるので、Linux版で行います。
Linux版とWindows版ではLinux版の方がかなり高速に動きます。
CentOS 6＋epleでは、QGISのバージョンが1.8であり、機能的に古いです。国土地理院地図を利用する場合、TileLayerPluginを利用すると便利ですが、利用できません。

準備

Windowsを利用している人が多いと思いますので、VMware Player、VirtualBoxなどを利用し、Linux CentOS 7をインストールしてください。
CentOS 7では、Desktop環境を利用しますので、GnomeDesktopでインストールしてください。

yum install -y epel-release
yum install -y qgis

これでgdal-libがインストールされていると思いますが、gdalライブラリは1.11.2のためGMTv5-netCDFは正しく読み込めません。gdalライブラリ2.0では修正されていますので、その修正を当てます。netcdfdataset.cppの450行目くらいです。
以下の下4行のif文を追加します。

netcdfdataset.cpp

/* -------------------------------------------------------------------- */
/*      Check for variable chunking (netcdf-4 only)                     */
/*      GDAL block size should be set to hdf5 chunk size                */
/* -------------------------------------------------------------------- */
# ifdef NETCDF_HAS_NC4
    int nTmpFormat = 0;
    size_t chunksize[ MAX_NC_DIMS ];
    status = nc_inq_format( cdfid, &nTmpFormat);
    if( ( status == NC_NOERR ) && ( ( nTmpFormat == NCDF_FORMAT_NC4 ) ||
          ( nTmpFormat == NCDF_FORMAT_NC4C ) ) ) {
        /* check for chunksize and set it as the blocksize (optimizes read) */
        status = nc_inq_var_chunking( cdfid, nZId, &nTmpFormat, chunksize );
        if( ( status == NC_NOERR ) && ( nTmpFormat == NC_CHUNKED ) ) {
             CPLError( CE_Failure, CPLE_AppDefined,
                      "setting block size to chunk size : %ld x %ld\n",
                      chunksize[nZDim-1], chunksize[nZDim-2]);
            nBlockXSize = (int) chunksize[nZDim-1];
            nBlockYSize = (int) chunksize[nZDim-2];
        }
	} 		
# endif

/* -------------------------------------------------------------------- */
/*      Force block size to 1 scanline for bottom-up datasets if        */
/*      nBlockYSize != 1                                                */
/* -------------------------------------------------------------------- */
    if( poNCDFDS->bBottomUp && nBlockYSize != 1 ) {
        nBlockXSize = nRasterXSize;
        nBlockYSize = 1;
    }
}

修正方法について

epelのソースから修正します。
修正方法についてはほかのサイトを見ていただいた方がよいと思いますので、簡単に記載します。
以下のURLからソースファイルをダウンロードします。
http://dl.fedoraproject.org/pub/epel/7/SRPMS/g/gdal-1.11.2-1.el7.src.rpm
(こっちは、リンクページ http://dl.fedoraproject.org/pub/epel/7/SRPMS/repoview/gdal.html)

以下でソースを修正して、rpmを作成してください。

rpmbuild -ba rpmbuild/SPECS/gdal.spec
rpm -ihv --force rpmbuild/RPMS/x86_64/gdal-libs-1.11.2-2.el7.centos.x86_64.rpm

※筆者は、正しくバージョンの修正を行っていないため、競合があると怒られてしました。
　--forceフラグを付けてごまかしています。また、yumのバージョン管理からも怒られるため注意してください。

読み込み速度が遅いという人は

GMTでは、GMTv5-netCDFとした場合、chunk sizeが設定されています。これは、netCDFファイルを全領域ではなく、一部の領域のデータを読み込みたいときに高速にデータを読み込めるように緯度ごとの並びではなく、緯度×経度の格子データで保存されています。

しかし、詳しいことは省略しますが、QGISが利用するgdalモジュールは、GMTv5-netCDFフォーマットにおいては、chunkによる読み込みがうまく行われません。データによっては、chunkによる保存がされていることで、同じ領域を何回も読み込むという無駄な処理がされています。(GDALではnetCDFの読み込みにおいてデータ順が反対になっているため、複数行の読み込みができず1行ずつの読み込みを行います。そのため、複数行の読み込みを行うときに、1行読み込むときに格子データを読み込み、次の1行を読み込むときにまた同じ格子データを読み込むという処理を行ってしまっています)

nccopyコマンド使って、chunkサイズを調整する。

ETOPO1のデータを利用しています。
ETOPO1のデータはGMTv4のnetcdfフォーマット(cfオプションで出力)ですので、grdreformatでGMTv5-netCDFフォーマットに変換しています。
nccopyコマンドを利用して、chunkサイズを変更します。経度方向は、経度方向すべてがよいです。また、緯度方向は、1がよいです。つまり、chunkを設定はしているが、格子ではなく、一列にしてしまっています。

通常、128～256の間でchunkサイズが設定されているので、数10～100倍近く高速になります。

gmt grdreformat ETOPO1_Ice_c_gmt4.grd ETOPO1_Ice_c_gmt4.nc=nf 
nccopy -c "lon/21600,lat/1" ETOPO1_Ice_c_gmt4.nc ETOPO1_Ice_c_gmt4_a.nc

参考ページ

Chunking Data: Why it Matters