たぶん動く...

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

GOSAT TANSO-FTS L2データをGeoJSON出力する

GOSAT(いぶき)のTANSO-FTSセンサで取得された、Level 2のHDF5ファイルからCO2、CH4のデータを取得しGeoJSONに出力するpythonプログラムを作成しました。 途中でgeopandasで変数を持つように設定しています。
NIES (国立環境研究所)がデータ提供されています。

データの取得

GOSAT Data Archive Service (GDAS)より入手ができます。ユーザ登録が必要です。HPにたどり着くまで分かりにくい....
GOSAT Data Archive Service

TANSO-FTS L2データの概要

プロダクトフォーマット説明書によると、 L2データはTANSO-FTSで観測された輝度スペクトルから算出された、CO2やCH4の鉛直プロファイル(TIR)、カラム量(SWIR)があります。 TANSO-FTSセンサが観測した、ポインティング範囲でのプロファイル、カラム量のようです。 ポインティング範囲は地上換算で10.5 kmになります。

開発環境

gdalのdockerで開発しています。以下の環境で開発しています。
Ubuntu 22.04.1 LTS
GDAL 3.4.1, released 2021/12/27
Python 3.10.6

ひとまず中身を見てみる

pythonのh5pyモジュールでSWIRL2CO2データのkey一覧と次元(shape)を取得してみます。

# coding:utf-8
import h5py

if __name__ == '__main__':

# Open File
    hdf_file = h5py.File('GOSATTFTS20220705_02C01SV0291R220705GU000.h5','r')

    for key_1st in list(hdf_file['Data'].keys()):
        for key_2st in list(hdf_file['Data'][key_1st].keys()):
            print('Data/'+key_1st+'/'+key_2st,hdf_file['Data'][key_1st][key_2st].shape)
##CLOSE HDF FILE
    hdf_file=None

出力結果

Data/auxiliaryParameter/aerosolOpticalThickness (376, 13, 3)
Data/auxiliaryParameter/dryAirTotalColumn (376,)
Data/auxiliaryParameter/surfaceAlbedo (376, 24)
Data/auxiliaryParameter/surfacePressure (376,)
Data/geolocation/airMass (376,)
Data/geolocation/footPrintLatitude (376, 36)
Data/geolocation/footPrintLongitude (376, 36)
Data/geolocation/height (376,)
Data/geolocation/landFraction (376,)
Data/geolocation/latitude (376,)
Data/geolocation/latitudeOriginal (376,)
Data/geolocation/longitude (376,)
Data/geolocation/longitudeOriginal (376,)
Data/geolocation/satelliteAttitude (376, 4)
Data/geolocation/satelliteAzimuth (376,)
Data/geolocation/satellitePosition (376, 3)
Data/geolocation/satelliteZenith (376,)
Data/geolocation/solarAzimuth (376,)
Data/geolocation/solarZenith (376,)
Data/geolocation/sunglintFlag (376,)
Data/mixingRatio/XCO2 (376,)
Data/mixingRatio/XCO2ExternalError (376,)
Data/mixingRatio/XCO2InterferenceError (376,)
Data/mixingRatio/XCO2RetrievalNoise (376,)
Data/mixingRatio/XCO2SmoothingError (376,)
Data/retrievalQuality/CO2DFS (376,)
Data/retrievalQuality/chi2 (376,)
Data/retrievalQuality/columnAveragingKernel (376, 15)
Data/retrievalQuality/iterations (376,)
Data/retrievalQuality/residualMeanSquare (376, 4)
Data/retrievalQuality/totalPostScreeningResult (376,)
Data/totalColumn/CO2TotalColumn (376,)
Data/totalColumn/CO2TotalColumnExternalError (376,)
Data/totalColumn/CO2TotalColumnInterferenceError (376,)
Data/totalColumn/CO2TotalColumnRetrievalNoise (376,)
Data/totalColumn/CO2TotalColumnSmoothingError (376,)

カラム量データなので、1次元配列のデータになっています。
位置情報が三種類用意されているようです。それぞれの内容は以下になります。

Dataset名 内容
footPrintLatitude, footPrintLongitude 2次元 観測点(視野範囲36点)の測地緯度経度(幾何補正、オルソ補正後)
latitude, longitude 1次元 観測点(FTS視野中心)の緯度経度(幾何補正、オルソ補正後)
latitudeOriginal, longitudeOriginal 1次元 観測点(FTS視野中心)の測地緯度(幾何補正、オルソ補正前)

latitude, longitudeが補正後のデータ点として使えそうです。 各データセットの内容説明については、L2プロダクトフォーマット説明書を参照してください。

データを取り出してみる

TANSO-FTS L2データのうち、XCO2 (乾燥空気に対するCO2の体積混合比)を取得してみます。

# coding:utf-8
import h5py

if __name__ == '__main__':

# Open File
    hdf_file = h5py.File('GOSATTFTS20220705_02C01SV0291R220705GU000.h5','r')
    targetdata=hdf_file['Data']['mixingRatio']['XCO2'][()]
    print(targetdata)

#CLOSE HDF FILE
    hdf_file=None

出力結果

[410.5042  410.02228 414.63004 412.96826 413.34454 414.24316 416.7913  
 409.75916 414.24655 412.99857 414.30426 414.12793 411.68893 414.52136  
 410.4508  410.15506 414.98773 412.5016  411.52545 411.5886  411.72226  
....

単位はフォーマット説明書によると単位はppmv (parts per million by volume)は体積百万分率で、乾燥した空気分子100万個中の当該ガスの分子数のことらしいです。

ファイル出力する

これまで見てきて、Latitude, Longitudeそしてデータの取り出し方が分かりました。 GeoJSONファイルに出力してみます。 単位と観測時刻も同時に出力されるようにしました。 将来的にデータ分析に使えるようにするため、GeoDataFrameに変換してからえ出力できるようにしました。

        #緯度経度情報を取り出す
        latitude=hdf_file['Data']['geolocation']['latitude'][()]
        longitude=hdf_file['Data']['geolocation']['longitude'][()]
        lat_lon_arr=np.column_stack((longitude,latitude))
    
        mid_key=find_key(input_file,dataset_name)
        targetdata=hdf_file['Data']['mixingRatio'][dataset_name][()]

        #単位を取り出す
        unit=str(hdf_file['Data'][mid_key][dataset_name].attrs['unit'][()][0].decode("utf-8"))
        #longNameを取り出す
        longName=str(hdf_file['Data'][mid_key][dataset_name].attrs['longName'][()][0].decode("utf-8"))

        d = {dataset_name: targetdata,
                'unit':[unit]*len(targetdata),
                'longName':[longName]*len(targetdata),
                'geometry': [Point(x) for x in lat_lon_arr],
                'time':self.metadata.time()
                }
        # GeoDataFrameに変換
        gdf = geopandas.GeoDataFrame(d, crs="EPSG:4326")

        # GeoJSON出力
        gdf.to_file(driver = 'GeoJSON', filename= '20220711.geojson')

今回作成したたプログラムは、以下になります。

GitHub - TTY6335/GOSAT_h5: program which extract data from GOSAT h5 files

出力結果を見てみる

QGISで出力してみました。
SWIRは観測点が昼間且つ、雲のかかってない領域のデータになるので、一日分でも観測点は多くないですね。

所感

まず、データの入手サイトがなかなか見つからず、そこが難儀しました。 (プロジェクトホームページからすぐに飛べないってどういうこと?) 最初はGDALですべて加工しようとしたのですが、h5ファイルの中に次元が異なるデータが存在するとうまく扱えないようで断念しました。 特に投影変換とかはしてないのでまだ分かりやすいほうだと思いました。

参考文献

NIES GOSAT TANSO-FTS SWIR レベル 2 プロダクトフォーマット説明書
NIES GOSAT TANSO-FTS TIR レベル 2 プロダクトフォーマット説明書 GitHub - TTY6335/GOSAT_h5: program which extract data from GOSAT h5 files