たぶん動く...

多分GIS系の人。あくまで個人的見解であり、所属団体を代表するものではありません。

GCOM-C L1BのHDF5 データを地図投影までしてみる。

JAXA GCOM-C SGLIのLevel 1のHDF5ファイルをgeotiffに変換と地図投影するpythonプログラムを作ってみました。

開発環境

gdal-pythonって結構依存性とかがあったりするからどんな環境で動いているのか載せておきます。
個人的にはdockerで環境構築するのが楽だと思います。
hub.docker.com
* python 3.7.4
* h5py 2.9.0
* hdf5 1.10.4
* numpy 1.16.5
* gdal 1.11.4

GCOM-C (しきさい)の概要

概要はこちら
suzaku.eorc.jaxa.jp
SGLIという光学観測センサを搭載しています。観測データから地球環境に関する物理量をたくさん得られるみたいです。

データの入手

JAXA G-Portalから入手できます。

gportal.jaxa.jp

google検索すると一番最初に出てこないのはなんでだろ。

L1Bデータの概要

GCOM-Cのデータ利用ガイドとか読むと、観測データままのL1Aから放射輝度を算出し、幾何補正をしたものになっているようです。 L1Bのプロダクトフォーマットガイド(ファイル構造の説明書的なもの)は以下です。 https://gportal.jaxa.jp/gpr/assets/mng_upload/GCOM-C/SGLI_Level1_Product_Format_Description_jp.pdf

L1Bフォーマットの説明

G-PortalからダウンロードできるHDF5形式のファイルについて画像を作るための情報をかんたんに解説します。
ものすごく簡単に言うとHDF5とは階層構造をもったバイナリファイルです。なんかしらのプログラムを書く知識がある人には、リストの中にラベルと配列が入っているようなものと言えば分かるかもしれません。
GCOM-C SGLI L1Bデータの場合以下の大きな枠組みにデータが別れます。

  1. Geometry_Data
    GCOM-Cの観測している地球表面での太陽の方向、観測した地球表面の緯度経度とかの情報が入っています。
  2. Image_Data
    L1Bの場合、どこが陸地か海か、品質の情報、そして、画像データがこの下に入っています。

データについてはこれくらい説明しとけば大丈夫かな。

プログラムのアプローチ

だいたいこれくらいでしょうか。

  1. HDF5ファイルから画像データを取り出す。
  2. 画像データをgeotiffにする。
  3. 画像データにGCPを付与する。
  4. GCPを元にgdal_translateでEPSG:4326に投影変換

プログラムを少し解説

少し流れを見ていきましょう。

lat_arr=np.array(f['Geometry_data']['Latitude'],dtype='float64')
lon_arr=np.array(f['Geometry_data']['Longitude'],dtype='float64')
#GCPのリストを作る
gcp_list=[]
for column in range(0,lat_arr.shape[0],20):
for row in range(0,lat_arr.shape[1],20): 
gcp=gdal.GCP(lon_arr[column][row],lat_arr[column][row],0,row*10,column*10)
gcp_list.append(gcp) 

HDF5内のGeometry_Dataをとりだし、画像上の位置row,columnとlat,lonの対応リスト化して保持します。
HDF5ファイル内では10ピクセル間隔で緯度経度の情報が入っています。全部抜き出してもgdaWarpの投影変換でエラーになるので、さらに飛び飛びに200ピクセル毎にGCPを設定します。

#RGBの3バンドだけ抽出する
    sgli_bands=['Lt_VN03','Lt_VN05','Lt_VN08']
    rgb_list=[]

    for sgli in (sgli_bands):
        print(sgli)
        band_image_arr=f['Image_data'][sgli]
        slope=f['Image_data'][sgli].attrs['Slope']
        offset=f['Image_data'][sgli].attrs['Offset']

        #numpyの行列にする
        band_image_arr=np.array(band_image_arr,dtype='uint16')
        #MSBを0にする
        band_image_arr=np.where(band_image_arr>=32768,band_image_arr-32768,band_image_arr)
        #15ビット目を0にする
        band_image_arr=np.where(band_image_arr>=16384,band_image_arr-16384,band_image_arr)
        #欠損値をnanで埋める
        band_image_arr=np.where(band_image_arr==16383,np.nan,band_image_arr)
        #大気上端放射輝度を求める
        Lt=slope*band_image_arr+offset
        Lt=np.array(Lt,dtype='float64')
        rgb_list.append(Lt)

Image_dataの配列をuint16で読み出すのがミソです。おそらく指定しないで読み込むとint16で読み込まれ、215以上の数値はマイナスになります。
そして最後にslopeを掛けてoffsetを足すことで大気上端放射輝度が求まります。

   output = gdal.GetDriverByName('GTiff').Create(output_file,col_size,row_size,band,dtype)
    output.GetRasterBand(1).WriteArray(rgb_list[2])
    output.GetRasterBand(2).WriteArray(rgb_list[1])
    output.GetRasterBand(3).WriteArray(rgb_list[0])
    wkt = output.GetProjection()
    output.SetGCPs(gcp_list,wkt)
    #GCPを使ってEPSG4326に投影変換
    output = gdal.Warp(output_file, output, dstSRS='EPSG:4326',tps = True,outputType=dtype)
    output = None    

gdalWarpを使って投影変換をします。 GCPのリストを与えるだけでなく、tps=Trueを入れておかないときちんと計算してくれません。
1分くらいかかるかな。

2020/09/21追記

gdalのバージョンによってはプログラムを実行するといつまでも終わらない現象が発生するようです。
jupyter notebookで実行するとエラーメッセージを取得できました。

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-54-c63ed80f64c5> in <module>
      7 output.SetGCPs(gcp_list,wkt)
      8 #与えたGCPを使ってEPSG4326に投影変換
----> 9 output = gdal.Warp(output_filename, output, dstSRS='EPSG:4326',tps = True,outputType=dtype)
     10 output.FlushCache()

~\.conda\envs\gdal\lib\site-packages\osgeo\gdal.py in Warp(destNameOrDestDS, srcDSOrSrcDSTab, **kwargs)
    623 
    624     if _is_str_or_unicode(destNameOrDestDS):
--> 625         return wrapper_GDALWarpDestName(destNameOrDestDS, srcDSTab, opts, callback, callback_data)
    626     else:
    627         return wrapper_GDALWarpDestDS(destNameOrDestDS, srcDSTab, opts, callback, callback_data)

~\.conda\envs\gdal\lib\site-packages\osgeo\gdal.py in wrapper_GDALWarpDestName(*args)
   3408 def wrapper_GDALWarpDestName(*args):
   3409     """wrapper_GDALWarpDestName(char const * dest, int object_list_count, GDALWarpAppOptions warpAppOptions, GDALProgressFunc callback=0, void * callback_data=None) -> Dataset"""
-> 3410     return _gdal.wrapper_GDALWarpDestName(*args)
   3411 class GDALVectorTranslateOptions(_object):
   3412     """Proxy of C++ GDALVectorTranslateOptions class."""

TypeError: in method 'wrapper_GDALWarpDestName', argument 4 of type 'GDALWarpAppOptions *'

エラーをググってみると同じような質問をしている人がいました。

gis.stackexchange.com
wrapper_GDALWarpDestName とあるのでdstSRSで指定している'EPSG:4326'がないという表示のようです。
また、別のページに詳しくエラーの原因がありました。

stackoverflow.com

どうもproj.dbが見つからず、EPSG:4326の定義が見つからないようです。
上のリンクのようにproj.dbのパスを通すのもアリだと思いますが、windows環境ではめんどくさいので、EPSG:4326の定義をそのまま貼り付けたほうが早いです。

output = gdal.Warp(output_filename, output, dstSRS='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs',tps = True,outputType=dtype)

githubも更新しました。

2020/10/10追記

出力されるtiffファイルのnodataの値を-999に変更しました。multithreadを有効化して高速化を図りました。

githubはこちらです。

GitHub - TTY6335/SGLI_L1: JAXA GCOM-C SGLI Level1 hdf5 to EPSG:4326 geotiff converter
githubwikiのページってグーグル検索で引っかからないのかな。