Zenn
🛰️

GCOM衛星画像の取得と変換

2023/11/16に公開
Python
衛星データ
GDAL
人工衛星
satellite
tech

JAXAが開発、打ち上げ、運用している人工衛星にGCOMシリーズという一連のプロジェクトがある。

GCOMはGlobal Change Observation Missionの略で、(環境という単語が抜けている気がするが、)地球環境の変動を観測するというミッションのもと2機の人工衛星を打ち上げている。

初号機は2012年に打ち上げられた「GCOM-W1」通称「しずく」と呼ばれる衛星で、水循環をターゲットに降水量や土壌に含まれる水分、海面水温などの観測を行っている。しずくは設計寿命の5年の2倍を超えて今も運用中だ。

二号機は2017年に打ち上げられた「GCOM-C」通称「しきさい」と呼ばれる衛星で、こちらは気候変動に焦点をあて、エアロゾル(大気中の微粒子)だったり、植生、地表面温度などを観測している。こちらも設計寿命の5年を超えたところだが、現役である。

この二つの衛星が観測したデータはJAXAが公開するG-Portalというプラットフォームから提供されており、ブラウザからのアクセスのほか、FTP経由でのダウンロードも可能となっている。

さてこのGCOM衛星のデータであるが、生データからのデータ処理の具合に応じてレベル1からレベル3までのレベル分けがされている。

プロダクトレベル
(画像はG-Portalより)

レベル1はほぼ生のデータなのであまり使う機会はないのだが、レベル2以上では海水面温度など、各地点での物理量が計算されたデータ、プロダクトが用意されており、用途に合わせてプロダクトを選択してダウンロードすることが可能になっている。

プロダクト定義
(画像はG-Portalより)

ダウンロードした画像は、各ピクセルの値として色ではなく、温度などの物理量を持っているため、地図上に重ね合わせて表示することで各地点での物理量が見てわかる。

地表面温度画像
(地球全体の地表面温度画像を表示した様子)

FTP経由で画像ダウンロード

ところで、大量の画像を使って解析を行いたい時など、いちいちブラウザ上からダウンロードするのは手間なので、プログラムでダウンロードしたいところ。

今回、地表面温度(LST: Land Surface Temperature)のレベル2(L2)データを利用したい。ところがJAXA Earth APIから公開されているのはL3のデータのみなので、FTPサーバー経由で取得する。

FTPサーバーのディレクトリ階層は以下のようになっている。

/standard/<プロジェクト名>/<衛星名>.<センサー名>/<プロダクト名>/<プロダクトバージョン>/<年>/<月>
  • standard: 標準プロダクト、この他にリクエストに応じて生成されるデータが置かれる/orderといったパスがあるが、今回は/standardで固定。
  • プロジェクト名: G-Portalは複数の衛星プロジェクトをホストしているため、衛星プロジェクトの名前が入る。今回はGCOM-C。
  • 衛星名: 人工衛星の名前、今回はGCOM-C。
  • センサー名: 人工衛星は複数のセンサー(≈カメラ)を搭載することがあり、それぞれデータを取得するため、センサー名を指定する。今回はSGLIという名前のセンサー。
  • プロダクト名: プロダクトの名称。今回は「レベル2の陸域(LAND)地表面温度(LST)統計量プロダクト」を意味するL2.LAND.LST_.Statisticsを指定する。
  • プロダクトバージョン: プロダクトのバージョン番号、今回はバージョン3。
  • 年、月: データを取得した期間

具体例としては以下のようなパスになる。

/standard/GCOM/GCOM-C.SGLI/L2.LAND.LST_.Statistics/3/2023/10

より詳細はG-Portalのユーザーマニュアルに記載がある。

では早速、Pythonのftplibパッケージを使ってFTP接続してみる。

from ftplib import FTP
ftp = FTP('ftp.gportal.jaxa.jp')
ftp.login('<G-Portalのログインユーザー名>', 'anonymous');
ftp.dir()

ログインの際、G-Portal上のログインユーザー名が必要となる。また、パスワードは固定でanonymousとなっている。

上記コマンドで以下のように表示されれば接続は成功。

drwxr-xr-x   5 3498     3498         4096 May 29  2018 nrt
drwxr-xr-x 331 3498     3498       131072 Nov 13 19:25 order
drwxr-xr-x 1655 3498     3498       131072 Nov 10 14:01 process
drwxr-xr-x  17 3498     3498         4096 Mar 18  2020 standard

では画像のダウンロードということで、standard/GCOM-C/GCOM-C.SGLI/L2.LAND.LST_.Statistics/3/2023/10/GC1SG1_20231008D08D_T0828_L2SG_LST_Q_3000.h5のパスにあるファイルを取得してみる。

with open('downloaded.h5', 'wb') as f:
	ftp.retrbinary("RETR /standard/GCOM-C/GCOM-C.SGLI/L2.LAND.LST_.Statistics/3/2023/10/GC1SG1_20231008D08D_T0828_L2SG_LST_Q_3000.h5",
		       f.write)

./downloaded.h5というファイル名でデータがローカルに保存されていればダウンロードは完了。

HDF5ファイルとその変換

さて、気づいたかもしれないが、ダウンロードしたファイルの拡張子がh5となっており、これはHDF5というファイル形式。

このファイル形式についてはあまり詳しくないので、Wikipediaの説明およびQiita記事を紹介する。

Hierarchical Data Format（階層的データ形式、略称：HDF）は、大量のデータを格納および構造化するために設計した一連のファイル形式（HDF4、HDF5）。米国立スーパーコンピュータ応用研究所が開発し、非営利法人である HDF グループが対応している。HDFグループは、HDF5 技術の継続的な開発と HDF で保存しデータの継続的な提供を確保することを使命としている。

(Wikipedia — Hierarchical Data Formatより)

https://qiita.com/simonritchie/items/23db8b4cb5c590924d95

このHDF5ファイルをGeoTiffに変換したいのだが、この操作にあたり、地理空間画像を取り扱うには必須のツールであるGDALを持ち出す。

GDALは Geospatial Data Abstraction Library の略で、地理空間画像を扱う際に便利な機能が詰まったコマンドラインプログラムおよび各言語のライブラリで、OSGeo (The Open Source Geospatial Foundation) という団体がOSSとして公開している。

https://sorabatake.jp/32245/

まず始めに変換対象となるデータについてのメタ情報を確認してみる。これにはgdalinfoというツールが用意されている。

gdalinfo downloaded.h5

このコマンドを叩くと出力される情報には、以下のように座標に関する情報、

Geometry_data_Grid_interval=0.0020833334 
Geometry_data_Grid_interval_unit=deg
Geometry_data_Image_projection=EQA (sinusoidal equal area) projection from 0-deg longitude
Geometry_data_Lower_left_latitude=0 
Geometry_data_Lower_left_longitude=100 
Geometry_data_Lower_right_latitude=0 
Geometry_data_Lower_right_longitude=110 
Geometry_data_Number_of_lines=4800 
Geometry_data_Number_of_pixels=4800 
Geometry_data_Upper_left_latitude=10 
Geometry_data_Upper_left_longitude=101.54266 
Geometry_data_Upper_right_latitude=10 
Geometry_data_Upper_right_longitude=111.69693

衛星やセンサー、プロダクトに関する情報、

Global_attributes_Algorithm_developer=Japan Aerospace Exploration Agency (JAXA)
Global_attributes_Algorithm_version=3.00
Global_attributes_Dataset_description=Level-2 temporal statistics file of LST
Global_attributes_Image_end_time=20231015 03:18:00.000
Global_attributes_Image_start_time=20231008 02:57:57.600
Global_attributes_Parameter_version=000.01  Global_attributes_Product_file_name=GC1SG1_20231008D08D_T0828_L2SG_LST_Q_3000.h5
Global_attributes_Product_level=Level-2
Global_attributes_Product_name=Level-2 Tile Statistics
Global_attributes_Product_version=3000
Global_attributes_Satellite=Global Change Observation Mission - Climate (GCOM-C)
Global_attributes_Sensor=Second-generation Global Imager (SGLI)
Global_attributes_Tile_number=0828

物理量に関する情報、

Image_data_LST_AVE_Data_description=Average of LST parameter (temporal statistics)
Image_data_LST_AVE_Error_DN=65535 
Image_data_LST_AVE_Maximum_valid_DN=65534 
Image_data_LST_AVE_Minimum_valid_DN=0 
Image_data_LST_AVE_Offset=0 
Image_data_LST_AVE_Slope=0.02 
Image_data_LST_AVE_Unit=Kelvin

そしてHDFは階層構造のため、複数のサブデータセットを持てる。最後のセクションでは、どんなデータセットを含むのかに関する情報が載っている。

Subdatasets:
  SUBDATASET_1_NAME=HDF5:"downloaded.h5"://Image_data/LST_AVE
  SUBDATASET_1_DESC=[4800x4800] //Image_data/LST_AVE (16-bit unsigned integer)
  SUBDATASET_2_NAME=HDF5:"downloaded.h5"://Image_data/LST_Date
  SUBDATASET_2_DESC=[4800x4800] //Image_data/LST_Date (8-bit character)
  SUBDATASET_3_NAME=HDF5:"downloaded.h5"://Image_data/LST_MAX
  SUBDATASET_3_DESC=[4800x4800] //Image_data/LST_MAX (16-bit unsigned integer)
  ...
コマンドのアウトプット全文
Driver: HDF5/Hierarchical Data Format Release 5
Files: downloaded.h5
Size is 512, 512
Metadata:
  Geometry_data_Grid_interval=0.0020833334 
  Geometry_data_Grid_interval_unit=deg
  Geometry_data_Image_projection=EQA (sinusoidal equal area) projection from 0-deg longitude
  Geometry_data_Lower_left_latitude=0 
  Geometry_data_Lower_left_longitude=100 
  Geometry_data_Lower_right_latitude=0 
  Geometry_data_Lower_right_longitude=110 
  Geometry_data_Number_of_lines=4800 
  Geometry_data_Number_of_pixels=4800 
  Geometry_data_Upper_left_latitude=10 
  Geometry_data_Upper_left_longitude=101.54266 
  Geometry_data_Upper_right_latitude=10 
  Geometry_data_Upper_right_longitude=111.69693 
  Global_attributes_Algorithm_developer=Japan Aerospace Exploration Agency (JAXA)
  Global_attributes_Algorithm_version=3.00
  Global_attributes_Dataset_description=Level-2 temporal statistics file of LST
  Global_attributes_Image_end_time=20231015 03:18:00.000
  Global_attributes_Image_start_time=20231008 02:57:57.600
  Global_attributes_Parameter_version=000.01
  Global_attributes_Product_file_name=GC1SG1_20231008D08D_T0828_L2SG_LST_Q_3000.h5
  Global_attributes_Product_level=Level-2
  Global_attributes_Product_name=Level-2 Tile Statistics
  Global_attributes_Product_version=3000
  Global_attributes_Satellite=Global Change Observation Mission - Climate (GCOM-C)
  Global_attributes_Sensor=Second-generation Global Imager (SGLI)
  Global_attributes_Tile_number=0828
  Image_data_Grid_interval=0.0020833334 
  Image_data_Grid_interval_unit=deg
  Image_data_Image_projection=EQA (sinusoidal equal area) projection from 0-deg longitude
  Image_data_Lower_left_latitude=0 
  Image_data_Lower_left_longitude=100 
  Image_data_Lower_right_latitude=0 
  Image_data_Lower_right_longitude=110 
  Image_data_LST_AVE_Data_description=Average of LST parameter (temporal statistics)
  Image_data_LST_AVE_Error_DN=65535 
  Image_data_LST_AVE_Maximum_valid_DN=65534 
  Image_data_LST_AVE_Minimum_valid_DN=0 
  Image_data_LST_AVE_Offset=0 
  Image_data_LST_AVE_Slope=0.02 
  Image_data_LST_AVE_Unit=Kelvin
  Image_data_LST_Date_Data_description=Dates of SGLI observation during 8-day (or 1-month) averaging period. e.g., if Bit3 & Bit5=1 (101000) then data of 4th (or 17-20th) and 6th (or 25-28th) days are used in the statistics.
  Image_data_LST_Date_Error_DN=255 
  Image_data_LST_Date_Maximum_valid_DN=254 
  Image_data_LST_Date_Minimum_valid_DN=0 
  Image_data_LST_Date_Offset=0 
  Image_data_LST_Date_Slope=1 
  Image_data_LST_Date_Unit=day in a 8-day or period in a month, see Data_description
  Image_data_LST_MAX_Data_description=Maximum of LST parameter (temporal statistics)
  Image_data_LST_MAX_Error_DN=65535 
  Image_data_LST_MAX_Maximum_valid_DN=65534 
  Image_data_LST_MAX_Minimum_valid_DN=0 
  Image_data_LST_MAX_Offset=0 
  Image_data_LST_MAX_Slope=0.02 
  Image_data_LST_MAX_Unit=Kelvin
  Image_data_LST_MIN_Data_description=Minimum of LST parameter (temporal statistics)
  Image_data_LST_MIN_Error_DN=65535 
  Image_data_LST_MIN_Maximum_valid_DN=65534 
  Image_data_LST_MIN_Minimum_valid_DN=0 
  Image_data_LST_MIN_Offset=0 
  Image_data_LST_MIN_Slope=0.02 
  Image_data_LST_MIN_Unit=Kelvin
  Image_data_LST_Ninput_Data_description=The number of input data
  Image_data_LST_Ninput_Error_DN=65535 
  Image_data_LST_Ninput_Maximum_valid_DN=65534 
  Image_data_LST_Ninput_Minimum_valid_DN=0 
  Image_data_LST_Ninput_Offset=0 
  Image_data_LST_Ninput_Slope=1 
  Image_data_LST_Ninput_Unit=NA
  Image_data_LST_Nused_Data_description=The number of data used in the statistics
  Image_data_LST_Nused_Error_DN=65535 
  Image_data_LST_Nused_Maximum_valid_DN=65534 
  Image_data_LST_Nused_Minimum_valid_DN=0 
  Image_data_LST_Nused_Offset=0 
  Image_data_LST_Nused_Slope=1 
  Image_data_LST_Nused_Unit=NA
  Image_data_LST_QA_flag_Data_description=QA flag (Bit0-1: 00:water, 11:land)
  Image_data_LST_QA_flag_Error_DN=255 
  Image_data_LST_QA_flag_Maximum_valid_DN=253 
  Image_data_LST_QA_flag_Minimum_valid_DN=0 
  Image_data_LST_QA_flag_Offset=0 
  Image_data_LST_QA_flag_Slope=1 
  Image_data_LST_QA_flag_Unit=NA
  Image_data_LST_RMS_Data_description=Root Mean Square of LST parameter (temporal statistics)
  Image_data_LST_RMS_Error_DN=65535 
  Image_data_LST_RMS_Maximum_valid_DN=65534 
  Image_data_LST_RMS_Minimum_valid_DN=0 
  Image_data_LST_RMS_Offset=0 
  Image_data_LST_RMS_Slope=0.02 
  Image_data_LST_RMS_Unit=Kelvin
  Image_data_Number_of_lines=4800 
  Image_data_Number_of_pixels=4800 
  Image_data_Upper_left_latitude=10 
  Image_data_Upper_left_longitude=101.54266 
  Image_data_Upper_right_latitude=10 
  Image_data_Upper_right_longitude=111.69693 
  Processing_attributes_Contact_point=JAXA/Earth Observation Research Center (EORC)
  Processing_attributes_Input_files=GC1SG1_20231008D01D_T0828_L2SG_LST_Q_3000.h5 GC1SG1_20231010D01D_T0828_L2SG_LST_Q_3000.h5 GC1SG1_20231011D01D_T0828_L2SG_LST_Q_3000.h5 GC1SG1_20231012D01D_T0828_L2SG_LST_Q_3000.h5 GC1SG1_20231013D01D_T0828_L2SG_LST_Q_3000.h5 GC1SG1_20231014D01D_T0828_L2SG_LST_Q_3000.h5 GC1SG1_20231015D01D_T0828_L2SG_LST_Q_3000.h5
  Processing_attributes_Input_File_Name_001=GC1SG1_20231008D01D_T0828_L2SG_LST_Q_3000.h5
  Processing_attributes_Input_File_Name_002=GC1SG1_20231010D01D_T0828_L2SG_LST_Q_3000.h5
  Processing_attributes_Input_File_Name_003=GC1SG1_20231011D01D_T0828_L2SG_LST_Q_3000.h5
  Processing_attributes_Input_File_Name_004=GC1SG1_20231012D01D_T0828_L2SG_LST_Q_3000.h5
  Processing_attributes_Input_File_Name_005=GC1SG1_20231013D01D_T0828_L2SG_LST_Q_3000.h5
  Processing_attributes_Input_File_Name_006=GC1SG1_20231014D01D_T0828_L2SG_LST_Q_3000.h5
  Processing_attributes_Input_File_Name_007=GC1SG1_20231015D01D_T0828_L2SG_LST_Q_3000.h5
  Processing_attributes_Processing_organization=JAXA/GCOM-C science project
  Processing_attributes_Processing_result=Good
  Processing_attributes_Processing_UT=20231015 17:03:38
Subdatasets:
  SUBDATASET_1_NAME=HDF5:"downloaded.h5"://Image_data/LST_AVE
  SUBDATASET_1_DESC=[4800x4800] //Image_data/LST_AVE (16-bit unsigned integer)
  SUBDATASET_2_NAME=HDF5:"downloaded.h5"://Image_data/LST_Date
  SUBDATASET_2_DESC=[4800x4800] //Image_data/LST_Date (8-bit character)
  SUBDATASET_3_NAME=HDF5:"downloaded.h5"://Image_data/LST_MAX
  SUBDATASET_3_DESC=[4800x4800] //Image_data/LST_MAX (16-bit unsigned integer)
  SUBDATASET_4_NAME=HDF5:"downloaded.h5"://Image_data/LST_MIN
  SUBDATASET_4_DESC=[4800x4800] //Image_data/LST_MIN (16-bit unsigned integer)
  SUBDATASET_5_NAME=HDF5:"downloaded.h5"://Image_data/LST_Ninput
  SUBDATASET_5_DESC=[4800x4800] //Image_data/LST_Ninput (16-bit unsigned integer)
  SUBDATASET_6_NAME=HDF5:"downloaded.h5"://Image_data/LST_Nused
  SUBDATASET_6_DESC=[4800x4800] //Image_data/LST_Nused (16-bit unsigned integer)
  SUBDATASET_7_NAME=HDF5:"downloaded.h5"://Image_data/LST_QA_flag
  SUBDATASET_7_DESC=[4800x4800] //Image_data/LST_QA_flag (8-bit character)
  SUBDATASET_8_NAME=HDF5:"downloaded.h5"://Image_data/LST_RMS
  SUBDATASET_8_DESC=[4800x4800] //Image_data/LST_RMS (16-bit unsigned integer)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

なお、Pythonから利用する場合は以下のようになる。

from osgeo import gdal
dataset = gdal.Open('downloaed.h5')
dataset.GetMetadata()

このサブデータセットの中で、今回は8日間の平均値を計算したデータである、HDF5:"downloaded.h5"://Image_data/LST_AVEをGeoTiffに変換していく。

変換にあたっては、GDALからTranslateというAPIが用意されており、この機能を使う。先にそのコマンドを示す。

gdal_translate \
    -of Gtiff \
    -a_srs ESRI:53008 \
    -a_ullr 11119492.53430028 1111949.2664455874 12231441.930901462 0.0 \
    HDF5:"downloaded.h5"://Image_data/LST_AVE \
    translated.tif

Pythonだと

lst_ave_dataset = dataset.GetSubDatasets()[0][0]
gdal.Translate(
        'translated.tif',
        lst_ave_dataset,
        format="GTiff",
        outputSRS="ESRI:53008",
        outputBounds=[11119492.53430028, 1111949.2664455874, 12231441.930901462, 0.0],
        noData=65535,
        creationOptions=['COMPRESS=LZW'])

この変換コマンドの中で、元画像の座標および投影法 (CRS: Coordinate Reference System、参照座標系という) を指定している部分があるので、それについて詳しくみていく。

サンソン図法と座標計算

先も少し触れたのだが、地理空間情報を扱う際にはその情報がどこにどういう形で投影されているのか、というメタ情報が必須になる。これを参照座標系、CRSと呼ぶのだが、地球が球体なのに対して、投影は平面に行われるため、様々な投影法がある。

今回取り扱うGCOM-Cのデータはサンソン図法、という投影法で書かれており、EPSGによるとESRI:53008というIDで認識される。

サンソン図法
(サンソン図法で書かれた地球 By Strebe - Own work, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=16115343)

先にgdalinfoで表示した画像のメタ情報の中には、対象領域の座標が含まれているので以下に抜き出す。

'Geometry_data_Lower_left_latitude': '0 ',
'Geometry_data_Lower_left_longitude': '100 ',
'Geometry_data_Lower_right_latitude': '0 ',
'Geometry_data_Lower_right_longitude': '110 ',
'Geometry_data_Upper_left_latitude': '10 ',
'Geometry_data_Upper_left_longitude': '101.54266 ',
'Geometry_data_Upper_right_latitude': '10 ',
'Geometry_data_Upper_right_longitude': '111.69693 ',

この中で今回必要なのは左上(Upper_left)と右下(Lower_right)の座標になるが、これらは緯度・経度の値となっているので、これをサンソン図法上の値に変換する必要がある。

座標変換の式はSinusoidal Projectionを参考にすると、地球半径をr=6,371,000[m]として、

x = r * 経度 * cos(緯度)
y = r * 緯度

となる。

ここで得られる値の単位は[m]になるので、ここに係数をかけて調整する。ESRI:53008では0.017453292519943295が係数として示されているのでこれらの値を掛けて座標を得る。

具体的に先のデータで計算すると、

x_upper_left = 6,371,000 * 101.54266 * cos(10 * π / 180) * 0.017453292519943295 = 11119492.53430028
y_upper_left = 6,371,000 * 10 * 0.017453292519943295 = 1111949.2664455874 
x_lower_right = 6,371,000 * 110 * cos(10 * π / 180) * 0.017453292519943295 = 12231441.930901462
y_lower_right = 6,371,000 * 0 * 0.017453292519943295 = 0.0

となり、先に示したgdal_translateコマンド内の座標の値と一致する。

再変換

ここまででも地図上で重ね合わせて確認することはできるが、最後に地図上で見る際に確認しやすいEPSG:4326という投影法に変換し直す。

この作業にはgdalwarpコマンドを利用する。

gdalwarp \
    -of Gtiff \
    -t_srs EPSG:4326 \
    translated.tif \
    warped.tif

Pythonなら

gdal.Warp(
        'warped.tif',
        'translated.tif',
        dstSRS='EPSG:4326'
)

得られた画像を地図上に重ね合わせて表示してみると以下のようになり、無事に変換できていることがわかる。


(マレーシア周辺に重なっている)


(よりクローズアップしてみた様子。島もバッチリ重なっている。)

まとめ

最近縁があって人工衛星画像を触る機会があり、いろいろ調べていたことをまとめようとは思っていたのだが、なんだかんだ後回しにしてしまっていたので、なんだか長文となってしまった。
結果として衛星データを始めて触る人には舌足らずだし、GCOM-Cデータの変換をしたいという人には冗長な内容となってしまったのは反省。

もう少し分かるようになったらもうちょっと上手くまとめられるかな、と思うのでわかったところからもう少しまめにまとめていきたい。

GCOM情報まとめ

最後に、GCOMに関する各種資料が各所に散らばっていて探しづらいのでここにまとめておく。

G Portal

GCOM-C

JAXA Earth API

サンソン図法

Discussion