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