たぶん動く...

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

トップ > gdal > GCOM-C L2のHDF5 データを地図投影までしてみる。

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

gdal python

JAXA GCOM-C SGLIのLevel 1のHDF5ファイルをgeotiffに変換と地図投影するpythonプログラムを作ってみました。
位置精度が高いものができたと思うので紹介します。 この情報はGCOM-C 候変動観測衛星「しきさい」(GCOM-C) データ利用ハンドブック 初版に基づいて作成しています!!!古いです。
A版に基づいた内容はこちら

データの概要について

データの形式とか、入っている順番とかはだいたいL1Bデータと同じな気がします。詳しくは、 tty6335.hatenablog.com を参照してください。

EQA図法

これが大変。考えた人凄いと思うよ。

f:id:TTY6335:20200809155042p:plain
図はSGLIユーザーガイドより引用

地球の表面を正弦波で囲った領域に投影して、その中を南北方向18個、東西方向に36個に分割している。で、GCOM-CのL2のデータはこの一つ一つの区切り単位で頒布されているんだ。

各ピクセルの緯度経度の決定

これが難しい。日本付近のデータを使ってちょっと詳しく見ていくよ。

https://gportal.jaxa.jp/gpr/img/br/GC1SG1_20200727A08D_T0529_L2SG_LST_Q_2000/GC1SG1_20200727A08D_T0529_L2SG_LST_Q_2000.jpg

タイル番号T0529ちょうど日本列島の南側が入っている区画になるよ。世界地図で見ている日本とはだいぶ形が違うことがわかるかと思います。
きちんと緯度経度を当てはめるのが大変。

ユーザーガイドの方法で緯度経度を当てはめる。

まずはユーザーガイドで紹介されている方法を実装してみます。

 #南極から北極までの総画素数
 NL_0=int(round(180.0/d))  
 #赤道における東西方向の総画素数
 NP_0=2*NL_0 
 #gdal_translateに与えるGCPのリスト
 gcp_list=[]

 for lin in range(0,lin_tile+1,50):
     for col in range(0,col_tile+1,50): 
         if(lin==lin_tile): 
             lin=lin-1
         if(col==col_tile): 
             col=col-1
         lin_total=lin+v_tile*lin_tile  
         col_total=col+h_tile*col_tile  
         lat=90.0-(lin_total+0.5)*d  
         NP_i=int(round(NP_0*math.cos(math.radians(lat))))  
         lon=360.0*(col_total+0.5-NP_0/2)/NP_i  
         gcp=gdal.GCP(round(lon,6),round(lat,6),,0,col+0.5,lin+0.5)

こんな感じでしょうか。gdalWarpで地図投影をする都合上、すべての画素でGCPをとっても地図投影できないので、間引きしてあります。画像が出力されるまで3分くらいかかりました。

また、以下の計算は隅のピクセルの角の緯度経度を求めています。例えば左上であれば一番左上のピクセルの一番左上、右上であれば一番右上のピクセルの一番右上を求めています。

計算結果
左上 右上 左下 右下
緯度 40.0000000 40.0000000 30.0000000 30.0000000
経度 143.5953220 156.6494420 127.0172200 138.5642400

合ってるのかな? HDF5ファイルの中に四隅の緯度経度が入っているからそれを見てみる。

HDF5ファイル内に入っている四隅の緯度経度座標
左上 右上 左下 右下
緯度 40.0000000 40.0000000 30.0000000 30.0000000
経度 143.5948 156.64888 127.01706 138.56407

緯度はあってそうですが、経度が合わなさそうですね。精度的には小数点第三位までしか合わないような感じがします。 右下での位置のずれが顕著ですね。

経度方向のずれを修正する。

たぶんこのユーザーガイドに載っている方法では、NP_iを四捨五入しているせいでずれが出ていると考察します。

        #南極から北極までの総画素数
        NL_0=int(round(180.0/d))
        #赤道における東西方向の総画素数
        NP_0=2*NL_0
        #gdal_translateに与えるGCPのリスト
        gcp_list=[]

        for lin in range(0,lin_tile+1,50):
            for col in range(0,col_tile+1,50):
                if(lin==lin_tile):
                    lin=lin-1
                if(col==col_tile):
                    col=col-1
                lin_total=lin+v_tile*lin_tile
                col_total=col+h_tile*col_tile
                lat=90.0-(lin_total+0.5)*d
                NP_i=NP_0*math.cos(math.radians(lat))
                lon=360.0*(col_total+0.5-NP_0/2)/NP_i
                gcp=gdal.GCP(round(lon,6),round(lat,6),0,col+0.5,lin+0.5)
計算結果
左上 右上 左下 右下
緯度 40.0000000 40.0000000 30.0000000 30.0000000
経度 143.5948020 156.6488750 127.0170590 138.5640650

だいぶあっているのではないでしょうか。もう一度、HDF5ファイル内に記載されている四隅の緯度経度を載せます。

HDF5ファイル内に入っている四隅の緯度経度座標
左上 右上 左下 右下
緯度 40.0000000 40.0000000 30.0000000 30.0000000
経度 143.5948 156.64888 127.01706 138.56407

小数点以下5桁までの精度で考えれば、かなりいいんじゃないんでしょうか。WGS84に合うように計算すればもっと精度が出るんじゃないかな。

とりあえず地図投影までした画像が以下になります。

f:id:TTY6335:20200809181658p:plain
EPSG4326に地図投影したLST

使い方

最後になりましたが、使い方です。必要なモジュールとかはL1Bの時と同じです。

python h5_2_tiff.py input.h5 Image_data output.tif

HDF5のファイル名は変えてあっても変換できます。地図投影には2分以上かかります。GCPを間引きすればもっと速くなるかもしれないね。
実行する環境によってはgdalWarpの実行中にメモリエラーを起こすかもしれません。

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も更新しました。

githubはこちら

github.com

参考文献

気候変動観測衛星「しきさい」(GCOM-C) データ利用ハンドブック https://gportal.jaxa.jp/gpr/assets/mng_upload/GCOM-C/GCOM-C_SHIKISAI_Data_Users_Handbook_jp.pdf
SGLI高次プロダクト フォーマット説明書
https://gportal.jaxa.jp/gpr/assets/mng_upload/GCOM-C/SGLI_Higher_Level_Product_Format_Description_jp.pdf