JASMES Tutorial

JASMES MODISデータをGeoTiFF変換したい (Pythonを使ったGeoTiFF変換)

JASMESサービスで提供しているバイナリ形式のMODISデータをGeoTiFFに変換するチュートリアルです。
バイナリ形式からGeoTiFFに変換することで、QGIS等を利用してデータの閲覧・編集・分析等をすることが出来ます。

注意

事前にJASMESデータのダウンロードが必要となり、JASMESデータ利用にはユーザ登録が必要です。

サンプルスクリプトはPython3を使用します。
本チュートリアルでは、Linux環境での実行例を示しますが、サンプルスクリプトについてはWindows環境でも使用できます。
詳細は本チュートリアルの「1.サンプルスクリプトの使用方法」をご確認ください。

また利用者が本ホームページの内容(スクリプト含む)を用いて行う一切の行為についてJAXAは責任を負うものではありません。

JASMESユーザ登録の詳細は以下のページをご参照ください。JASMESユーザはJASMES関連サービスで公開している全てのデータを利用できます。
本チュートリアルでは、MODISデータ(*_le)を物理量に変換したGeoTiFFを作成します。
JASMESユーザ登録 JASMESデータ取得方法 JASMESデータ一覧

0. はじめに
1. サンプルスクリプトの使用方法

1-1. スクリプトダウンロード
1-2. 実行環境
1-3. 対象データ
1-4. スクリプト実行方法

2. スクリプト内容説明

2-1. Pythonでバイナリファイルの読み込み
2-2. PythonでSlope値、Offset値を使用して物理量に変換
2-3. PythonでGeoTiFFの出力 (GDALライブラリ使用)

3. GeoTiFF利用例 QGISでの開き方

0. はじめに

当チュートリアルは目的に応じて以下のようにご利用ください。

・GeoTiFF作成をしたい場合
「1. サンプルスクリプトの使用方法」に沿って、スクリプト実行をしてください。

・サンプルスクリプトの処理の内容を知りたい場合
「2. スクリプト内容説明」をご参照ください。

・GeoTiFFの利用方法としてQGISでの開き方を知りたい場合
「3. GeoTiFF利用例 QGISでの開き方」をご参照ください。

1. サンプルスクリプトの使用方法

1-1. スクリプトダウンロード

当チュートリアルで使用するサンプルスクリプトです。以下のボタンよりダウンロードしてください。
読み込みファイルに応じてスクリプトを複数用意しております。どのスクリプトを使用するかは、「1-3.対象データ」の一覧をご確認ください。
当チュートリアルは、(1)のスクリプトを使用した例となりますが、基本的には読み込み部分以外のGeoTiFF変換方法は共通で、スクリプト使用方法も同様となります。
スクリプト毎の特記事項は、スクリプトダウンロード時に付属したREADME_*.txtをご確認ください。

サンプルスクリプトについて

サンプルスクリプトはPython3を使用し、Pythonライブラリはctypes、NumPy、osgeo(gdal、osr)を使用します。
本チュートリアルでは、Linux環境での実行例を示しますが、サンプルスクリプトについてはWindows環境でも使用できます。

サンプルスクリプトの改変は自由で再配布も可能です。

サンプルスクリプト(1) (ファイル名「*_le」) ダウンロードサンプルスクリプト(2) (ファイル名「*_8b」) ダウンロードサンプルスクリプト(3) (複数バンド用) ダウンロードサンプルスクリプト(4) (lst用) ダウンロードサンプルスクリプト(5) (wf、wst用) ダウンロードサンプルスクリプト(6) (rgb用) ダウンロード

1-2. 実行環境

当チュートリアルではLinux環境での実行方法を記載しています。
またサンプルスクリプト内の処理では、Python3を使用しています。

サンプルスクリプト動作確認時の実行環境、使用するPythonライブラリは以下です。

実行環境

Linux：Red Hat Enterprise Linux(8.4)
Python：Python 3.6.8
Python標準ライブラリ：ctypes
Python外部ライブラリ：NumPy, osgeo(gdal, osr) ※要インストール

1-3. 対象データ

JASMESデータを取得する場合は、こちらからユーザ登録が必要となります。ユーザ登録後のデータの取得方法はこちらをご参照ください。

サンプルスクリプトを使用できるプロダクトと対応するスクリプト番号は以下の通りです。
JASMES提供データの一覧に、それぞれのプロダクトの領域、解像度、統計期間(Daily、Monthly等)の詳細がありますので、必要に応じてご参照ください。

No. Webページ名衛星/センサ圏プロダクトスクリプト備考

(1) (2) (3) (4) (5) (6)

1 JASMES Terra・Aqua/MODIS 大気 aot エアロゾル光学的厚さ ○ - - - - -

2 JASMES Terra・Aqua/MODIS 大気 aerosol エアロゾルの光学的厚さ - - ○ - - - 複数バンド

3 JASMES Terra・Aqua/MODIS 大気 alph オングストローム指数 ○ - - - - -

4 JASMES Terra・Aqua/MODIS 大気 cfr 雲フラグ - - - - - - HDFファイルのため、対象外

5 JASMES Terra・Aqua/MODIS 大気 dpar 直達PAR比率 ○ - - - - -

6 JASMES Terra・Aqua/MODIS 大気 par 光合成有効放射(PAR) ○ - - - - -

7 JASMES Terra・Aqua/MODIS 大気 ptw 水蒸気量 ○ - - - - -

8 JASMES Terra・Aqua/MODIS 大気 tip 直達透過光量 ○ ○ - - - -

9 JASMES Terra・Aqua/MODIS 大気 rpar PARおよび太陽照度に>よる表面反射率 ○ - - - - -

10 JASMES Terra・Aqua/MODIS 大気 swr 短波放射量 ○ ○ - - - -

11 JASMES Terra・Aqua/MODIS 大気 ta1 エアロゾルの光学的厚さ[WK= 412.46 nm] ○ - - - - -

12 JASMES Terra・Aqua/MODIS 大気 taua エアロゾルの光学的厚さ ○ - - - - -

13 JASMES Terra・Aqua/MODIS 大気 uva UV-A ○ ○ - - - -

14 JASMES Terra・Aqua/MODIS 大気 uvb UV-B ○ ○ - - - -

15 JASMES Terra・Aqua/MODIS 陸 lst 地表面温度 - - - ○ - -

16 JASMES Terra・Aqua/MODIS 陸 ndvi 植生分布 ○ - - - - -

17 JASMES Terra・Aqua/MODIS 陸 wf 森林火災 - - - - ○ - データにヘッダ無し

18 JASMES Terra・Aqua/MODIS 陸 wst 植生乾燥指数 - - - - ○ - データにヘッダ無し

19 JASMES Terra・Aqua/MODIS 海洋 chla クロロフィルa濃度 ○ - - - - -

20 JASMES Terra・Aqua/MODIS 海洋 olst 海面水温(+地表面温度) ○ - - - - -

21 JASMES Terra・Aqua/MODIS 海洋 sst 海面水温 ○ - - - - -

22 JASMES Terra・Aqua/MODIS 海洋 apg 水分子以外の粒子と溶存物質の吸収係数 ○ - - - - -

23 JASMES Terra・Aqua/MODIS 海洋 bbp 水分子以外の粒子の後方散乱係数 ○ - - - - -

24 JASMES Terra・Aqua/MODIS 雪氷 snwcfr 雪氷分布 ○ - - - - -

25 JASMES Terra・Aqua/MODIS - rgb RGB - - - - - ○ RGB合成スクリプト

26 JASMES Terra・Aqua/MODIS - multi PAR等のパラメータが複数格納 - - ○ - - - 複数バンド

No.	Webページ名	衛星/センサ	圏	プロダクト	スクリプト	備考
(1)	(2)	(3)	(4)	(5)	(6)
1	JASMES	Terra・Aqua/MODIS	大気	aot	エアロゾル光学的厚さ	○	-	-	-	-	-
2	JASMES	Terra・Aqua/MODIS	大気	aerosol	エアロゾルの光学的厚さ	-	-	○	-	-	-	複数バンド
3	JASMES	Terra・Aqua/MODIS	大気	alph	オングストローム指数	○	-	-	-	-	-
4	JASMES	Terra・Aqua/MODIS	大気	cfr	雲フラグ	-	-	-	-	-	-	HDFファイルのため、対象外
5	JASMES	Terra・Aqua/MODIS	大気	dpar	直達PAR比率	○	-	-	-	-	-
6	JASMES	Terra・Aqua/MODIS	大気	par	光合成有効放射(PAR)	○	-	-	-	-	-
7	JASMES	Terra・Aqua/MODIS	大気	ptw	水蒸気量	○	-	-	-	-	-
8	JASMES	Terra・Aqua/MODIS	大気	tip	直達透過光量	○	○	-	-	-	-
9	JASMES	Terra・Aqua/MODIS	大気	rpar	PARおよび太陽照度に>よる表面反射率	○	-	-	-	-	-
10	JASMES	Terra・Aqua/MODIS	大気	swr	短波放射量	○	○	-	-	-	-
11	JASMES	Terra・Aqua/MODIS	大気	ta1	エアロゾルの光学的厚さ[WK= 412.46 nm]	○	-	-	-	-	-
12	JASMES	Terra・Aqua/MODIS	大気	taua	エアロゾルの光学的厚さ	○	-	-	-	-	-
13	JASMES	Terra・Aqua/MODIS	大気	uva	UV-A	○	○	-	-	-	-
14	JASMES	Terra・Aqua/MODIS	大気	uvb	UV-B	○	○	-	-	-	-
15	JASMES	Terra・Aqua/MODIS	陸	lst	地表面温度	-	-	-	○	-	-
16	JASMES	Terra・Aqua/MODIS	陸	ndvi	植生分布	○	-	-	-	-	-
17	JASMES	Terra・Aqua/MODIS	陸	wf	森林火災	-	-	-	-	○	-	データにヘッダ無し
18	JASMES	Terra・Aqua/MODIS	陸	wst	植生乾燥指数	-	-	-	-	○	-	データにヘッダ無し
19	JASMES	Terra・Aqua/MODIS	海洋	chla	クロロフィルa濃度	○	-	-	-	-	-
20	JASMES	Terra・Aqua/MODIS	海洋	olst	海面水温(+地表面温度)	○	-	-	-	-	-
21	JASMES	Terra・Aqua/MODIS	海洋	sst	海面水温	○	-	-	-	-	-
22	JASMES	Terra・Aqua/MODIS	海洋	apg	水分子以外の粒子と溶存物質の吸収係数	○	-	-	-	-	-
23	JASMES	Terra・Aqua/MODIS	海洋	bbp	水分子以外の粒子の後方散乱係数	○	-	-	-	-	-
24	JASMES	Terra・Aqua/MODIS	雪氷	snwcfr	雪氷分布	○	-	-	-	-	-
25	JASMES	Terra・Aqua/MODIS	-	rgb	RGB	-	-	-	-	-	○	RGB合成スクリプト
26	JASMES	Terra・Aqua/MODIS	-	multi	PAR等のパラメータが複数格納	-	-	○	-	-	-	複数バンド

1-4. サンプルスクリプトの使用方法

入力ファイルパスと出力ファイルパスを引数として、サンプルスクリプトを実行します。
例ではファイル名「*_le」用のスクリプトを使用していますが、「*_8b」用でも使用方法は同様です。
その他それぞれのサンプルスクリプトの冒頭に使用方法の記載がありますので合わせてご確認ください。

$ ./MODIS_le_bin2geotiff.py <Input> <Output> <Area>

引数1 Input :
入力となるバイナリファイルをフルパスで指定します。

引数2 Output :
出力するGeoTiFFファイル名を指定します。ファイル名のみ指定するとカレントディレクトリに出力されます。

引数3 Area :
入力するファイルと対応した領域を指定します。指定可能な値はGlobal、Japan、Thaiです。
「*_le」用、「*_8b」用以外についてはサンプルスクリプトの冒頭の使用方法の記載をご確認ください。

コマンド実行例は以下の通りです。出力ファイル名は入力ファイルに「.tif」をつけたものとしていますが、任意のファイル名で問題ありません。
[DIR] : データの入出力ディレクトリ

$ ./MODIS_le_bin2geotiff.py [DIR]/MDS02SSH_A20230101Av1_v811_7200_3601_CHLA_le [DIR]/MDS02SSH_A20230101Av1_v811_7200_3601_CHLA_le.tif

2. スクリプト内容説明

ここから、サンプルスクリプトの各種処理について説明します。

2-1. Pythonでバイナリファイルの読み込み

「1-3. 対象データ」で示したJASMES MODISデータのフォーマットに合わせ、ヘッダの指定をしています。
JASMES MODISデータは、ヘッダ部分にピクセルサイズ等の必要な情報を掲載するフォーマットとなっているため、ヘッダ読み込み処理を行う必要があります。
こちらはJASMES MODISデータ特有の内容となりますのでご注意ください。

# Header format
class headerFormat(c.Structure):
  _fields_ = [
    ("npixel" , c.c_char*6 ),
    ("nline"  , c.c_char*6 ),
    ("lon_min", c.c_char*8 ),
    ("lat_max", c.c_char*8 ),
    ("reso"   , c.c_char*8 ),
    ("slope"  , c.c_char*12),
    ("offset ", c.c_char*12)
  ]

引数(入力バイナリファイルパス、出力GeoTiFFファイルパス、対象領域)の取得をしています。

argvs   = sys.argv

# argument ERROR
if len(argvs) != 4:
  print("argument error: 1:input file, 2:output file, 3:area (Global / Japan / Thai)")
  sys.exit(9)
else:
  binfile = argvs[1]
  output  = argvs[2]
  area    = argvs[3]

バイナリファイル読み込みをしています。
ファイルに応じてヘッダ情報をformat_dataに読み込んだ後、ヘッダ分をスキップしてデータ部をbinarrに読み込んでいます。

# 1 line * 2byte
HEADER_SIZE = X_SIZE*2
format_data = headerFormat()

# Binary file open
with open(binfile, 'rb') as f:
  # Read header
  f.readinto(format_data)
  # Header skip
  f.seek(HEADER_SIZE)
  # Read binary data as unsiged-2byte integer
  binarr = np.fromfile(f, dtype=np.uint16)

format_dataに読み込んだヘッダの必要な値の取得をしています。
取得した値は以下の通りです。
・SLOPE / OFFSETは物理量変換時(DN * SLOPE + OFFSET)に使用。
・X_SIZE / Y_SIZE / UpperLeft_LON / UpperLeft_LAT / RESO はGeoTiFF出力時に使用。

# Header value set
X_SIZE        = np.uint16(format_data.npixel)
Y_SIZE        = np.uint16(format_data.nline)
UpperLeft_LON = np.float32(format_data.lon_min)
UpperLeft_LAT = np.float32(format_data.lat_max)
RESO          = np.float32(format_data.reso)
SLOPE         = np.float32(format_data.slope)
OFFSET        = np.float32(format_data.offset)

# pixel center -> top left corner of the pixel
UpperLeft_LON = UpperLeft_LON - (RESO * 0.5)
UpperLeft_LAT = UpperLeft_LAT + (RESO * 0.5)

2-2. PythonでSlope値、Offset値を使用して物理量に変換

読み込んだ値はDN値となるため、「DN * SLOPE + OFFSET」で物理量に変換して、配列に格納します。
また65535は「Error_DN」となるため、nanに置き換えしています。(「*_8b」の場合は255が「Error_DN」)
物理量変換でデータ型がfloatとなりますのでデータ型も変換しています。

# Reshape pixel x line data
orgdat = binarr.reshape(Y_SIZE,X_SIZE).copy()
# value = DN * SLOPE + OFFSET
otemp  = orgdat * SLOPE + OFFSET
# convert error value to nan
otemp[orgdat==65535] = np.nan

2-3. PythonでGeoTiFFの出力 (GDALライブラリ使用)

物理量変換したpython配列をPythonのGDALライブラリでGeoTiFF形式で出力します。
まずはGeoTiFFを作成するGDALのドライバを作成します。

出力GeoTiFF名、配列の列数、配列の行数、バンド数、データ型を指定しています
・driver.Create(出力ファイル, 配列の列数, 配列の行数, バンド数, データ型)

# output GeoTiFF
driver = gdal.GetDriverByName('GTiff')

# outputfile, cols, rows, band num, data type
cols      = otemp.shape[1]
rows      = otemp.shape[0]
outRaster = driver.Create(output, cols, rows, 1, gdal.GDT_Float32)

SetGeoTransformに座標値を渡して設定をします。
・outRaster.SetGeoTransform((左上の経度, ピクセルサイズ(x), 回転, 左上の緯度, 回転, ピクセルサイズ(y)))

# upper_left_longitude, pixel size(x), rotate, upper_left_latitude, rotate, pixel size(y)
dx = RESO
dy = RESO*-1
outRaster.SetGeoTransform((UpperLeft_LON, dx, 0, UpperLeft_LAT, 0, dy))

バンドを指定します。

# band
outband = outRaster.GetRasterBand(1)
outband.WriteArray(otemp)

投影法はEPSG:4326を指定しています。

# EPSG
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromEPSG(4326)
outRaster.SetProjection(outRasterSRS.ExportToWkt())

ファイルに書き出しをしてGeoTiFFを出力します。

outband.FlushCache()
outRaster = None

print("### END output : "+output+" ###")
sys.exit(0)

3. GeoTiFF利用例 QGISでの開き方

当チュートリアルで作成したGeoTiFFは、無料のGISソフトであるQGISで、そのまま開いてデータの閲覧、編集、分析等が出来ます。
以下にQGISの使用方法の記載がありますので、必要に応じてご確認ください。

・クロロフィルa濃度の表示(QGIS) (しきさいポータル)
対数目盛（ログスケール）での表示方法の記載があります。CHLA等対数目盛の方が見やすい場合もありますので、必要に応じてご覧ください。
当チュートリアルでは既に物理量換算を実施済みのため、「DN * Slope + Offset」の計算は不要です。
クロロフィルa濃度の表示(QGIS) (しきさいポータル)