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

GCOM-Cのデータ利用ハンドブックのA版が2020/12/22に公開されていました。 A版では、L2タイルの緯度経度の計算方法が新しくなりました。その方法を実装してみたので紹介したいと思います。また、しばらく使ってみて修正すべき点があったので改修しましたので紹介したいと思います。

GCOM-Cのプロダクトについて
新しい緯度経度の計算方法について
- 実装
- 精度の確認
その他
- Nodataの追加
- リサンプリング手法の指定
githubはこちら
- 参考文献

GCOM-Cのプロダクトについて

前回の記事を参照してください。

tty6335.hatenablog.com

新しい緯度経度の計算方法について

データ利用ハンドブックの初版で紹介していた計算方法と一番違うところは経度(lon)の計算方法です。

実装

とりあえずpythonで実装してみます。

        #SGLI/L2であれば固定
        #縦方向の総タイル数
        vtilenum=18
        #横方向の総タイル数
        htilenum=36
    
        #緯度方向(dlin)、経度方向(dcol)
        #それぞれの1画素の大きさ
        #d=dlin=dcol
        d=180.0/lintile/vtilenum
    
        #求めたりタイル番号の左上画素の中心の緯度[deg]は、
        #1タイルあたりの角度が10[deg]であることから、
        lat0=90.0-vtile*10-d/2
    
        #求めたいタイル番号の左上画素の中心の経度[deg]は、
        #1タイルあたりの角度が10[deg]であることから、
        lon0=-180.0+htile*10+d/2

        #gdal_translateに与えるGCPのリスト
        gcp_list=[]
        
        for lin in range(0,lintile+1,100):
                lat=lat0-lin*d
                r=np.cos(np.radians(lat))
                for col in range(0,coltile+1,100):
                        if(lin==lintile):
                                lin=lin-1
                        if(col==coltile):
                                col=col-1
                        lon=(lon0+col*d)/r
                        gcp=gdal.GCP(round(lon,6),round(lat,6),0,col+0.5,lin+0.5)

実装してみた印象としては、初版よりA版のほうが実装が簡単（計算が簡単）でした。
forの二重ループは100画素ステップで計算しています。
50画素ステップでもgdalwarpの投影変換ができますが、CPUの性能によると思いますが、3分以上時間がかかります。(GPUを有効化すると速いのかな？)

精度の確認

A版の手法は変に四捨五入したりしているところがないので、手計算で試した時点でかなりよさそうなきがしました。四隅の緯度経度を計算してみます。対象のタイルはT0529です。

左上であれば一番左上のピクセルの一番左上、右上であれば一番右上のピクセルの一番右上を求めています。

計算結果

	左上	右上	左下	右下
緯度	40.0000000	40.0000000	30.0000000	30.0000000
経度	143.5948018	156.6488747	127.0170592	138.5640646

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

HDF5ファイル内に入っている四隅の緯度経度座標

	左上	右上	左下	右下
緯度	40.0000000	40.0000000	30.0000000	30.0000000
経度	143.5948	156.64888	127.01706	138.56407

かなり初版の計算方法より合っているようですね。
計算結果と、HDF5ファイル内で記載されている緯度経度を比較すると、小数点5桁目から差がみられていますが、無視できる誤差だと考えています。というのも、小数点以下4桁目まで一致していれば、小数点5桁目以下の違いは角度にして0.0001度以下、距離にして約10m程度*1となり、GCOM-Cのプロダクトの分解能250mよりも十分に小さいと考えられます。

その他

使っているうちに修正すべき点がいくつかあったので修正しました。

Nodataの追加

gdalwarpで投影変換の結果余白域が0としてはいっていました。NDVIなど数値が0をまたぐプロダクトでは、解析に支障がでるので、NodataになるようgdalWarpのオプションを追加しました。

リサンプリング手法の指定

日本の夏場のLSTを解析していた時のこと、地表面温度が60以上のピクセルが現れてこれはおかしいだろ。とサンプリング手法をnearest neighbor法を指定しました。

オプションを追加したgdalWarpは以下です。

output=gdal.Warp(output_file, output, dstSRS='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs',srcNodata=np.nan,dstNodata=np.nan,tps = True, outputType=dtype,multithread=True,resampleAlg=gdalconst.GRIORA_NearestNeighbour)

githubはこちら

github.com

参考文献

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

*1:緯度1度を111kmとして単純計算

たぶん動く...

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