JAXA GCOM-C SGLIのLevel 1のHDF5ファイルをgeotiffに変換と地図投影するpythonプログラムを作ってみました。
- 開発環境
- GCOM-C (しきさい)の概要
- データの入手
- L1Bデータの概要
- L1Bフォーマットの説明
- プログラムのアプローチ
- プログラムを少し解説
- 2020/09/21追記
- 2020/10/10追記
- githubはこちらです。
開発環境
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という光学観測センサを搭載しています。観測データから地球環境に関する物理量をたくさん得られるみたいです。
データの入手
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データの場合以下の大きな枠組みにデータが別れます。
- Geometry_Data
GCOM-Cの観測している地球表面での太陽の方向、観測した地球表面の緯度経度とかの情報が入っています。 - Image_Data
L1Bの場合、どこが陸地か海か、品質の情報、そして、画像データがこの下に入っています。
データについてはこれくらい説明しとけば大丈夫かな。
プログラムのアプローチ
だいたいこれくらいでしょうか。
プログラムを少し解説
少し流れを見ていきましょう。
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'がないという表示のようです。
また、別のページに詳しくエラーの原因がありました。
どうも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
githubのwikiのページってグーグル検索で引っかからないのかな。