JASMES Tutorial

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

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

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

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

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

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

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

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

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

1-2. 実行環境

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

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

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等のパラメータが複数格納	-	-	○	-	-	-	複数バンド

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-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)

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

・QGISを使ったGCOM-C(しきさい)プロダクトの画像化手順 (G-Portal)
 当チュートリアルで既に物理量換算とGeoTiFF変換は実施済みのため、STEP5以降と「カラーランプの変換方法」をご参照ください。
QGISを使ったGCOM-C(しきさい)プロダクトの画像化手順 (G-Portal)
・クロロフィルa濃度の表示(QGIS) (しきさいポータル)
 対数目盛（ログスケール）での表示方法の記載があります。CHLA等対数目盛の方が見やすい場合もありますので、必要に応じてご覧ください。
 当チュートリアルでは既に物理量換算を実施済みのため、「DN * Slope + Offset」の計算は不要です。
クロロフィルa濃度の表示(QGIS) (しきさいポータル)