JASMES Tutorial

JASMES GCOM-C/SGLI、VIIRSデータ(netCDF) をGeoTiFF変換したい (Python)

JASMESサービスで提供しているnetCDF形式のGCOM-C/SGLIデータ(準リアル・標準)、VIIRSデータをGeoTiFFに変換するチュートリアルです。
また、作成したGeoTiFFをQGISで表示する方法も紹介します。

本チュートリアルでは、Linux環境でGCOM-C/SGLI標準データ(netCDF)を物理量に変換したGeoTiFFを作成します。
事前にJASMESデータのダウンロードが必要となり、JASMESデータ利用にはユーザ登録が必要です。詳細はこちらをご参照ください。

GeoTiFF変換方法として以下の2通りの説明ページを用意しております。当ページでは(1)のPythonスクリプトでの方法を示します。
(1) Linux環境でPythonのサンプルスクリプトを使用したGeoTiFF変換方法 (当チュートリアル)
(2) Windows環境でGDALコマンド及びパラメータを用いたGeoTiFF変換方法 (しきさいポータルFAQ)

注意

事前に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. 使用データの確認
1-5. サンプルスクリプトの使用方法

2. スクリプト内容説明

2-1. PythonでnetCDFの読み込み
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. スクリプトダウンロード

当チュートリアルで使用するサンプルスクリプトです。以下のボタンよりダウンロードしてください。

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

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

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

サンプルスクリプト (GCOM-C/SGLIデータ・VIIRSデータ) ダウンロード

1-2. 実行環境

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

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

実行環境

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

1-3. 対象データ

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

サンプルスクリプトを使用できるデータは以下をご確認ください。対象データ一覧

1-4. 使用データの確認

(a) SeaDASを利用しての確認
netCDFを確認するツールとして、SeaDAS(Windows、Mac、Linuxに対応)を使用する方法については、こちらをご確認ください。 JASMES SGLI準リアルデータ SeaDAS読み込み編
(b) ncdumpを利用しての確認
本チュートリアルでは、netCDFのコマンド ncdumpでの確認方法を以下に示します。「ncdump -h」を実行し、データのヘッダ情報を確認します。
「long_name」に記載の通り「netCDFに格納されたDN値 * scale_factor + add_offset」で物理量に変換できます。
チュートリアル内の処理で使用する情報はハイライトで示した項目です。

$ ncdump -h [格納ディレクトリ]/GC1SG1_20230101D01D_J0000_L2SG_CHLAQ_3000.nc
 netcdf GC1SG1_20230101D01D_J0000_L2SG_CHLAQ_3000 {
 dimensions:
	longitude = 10800 ;
	latitude = 10400 ;
	start_time = 1 ;
	end_time = 1 ;

 variables:
	float longitude(longitude) ;
		longitude:long_name = "longitude" ;
		longitude:units = "degree_east" ;
	float latitude(latitude) ;
		latitude:long_name = "latitude" ;
		latitude:units = "degree_south" ;
	double start_time(start_time) ;
		start_time:long_name = "observation start time" ;
		start_time:units = "UTC" ;
		start_time:standard_name = "yyyymmdd.hhmmss" ;
	double end_time(end_time) ;
		end_time:long_name = "observation end time" ;
		end_time:units = "UTC" ;
		end_time:standard_name = "yyyymmdd.hhmmss" ;
 	ushort CHLA(latitude, longitude) ;
		CHLA:long_name = "Chllorophyll-a concentration (CHLA) = DN * scale_factor + add_offset [mg m^-3]" ;
		CHLA:units = "mg m^-3" ;
		CHLA:scale_factor = 0.0016f ;
		CHLA:add_offset = 0.f ;
		CHLA:valid_min = 0 ;
		CHLA:valid_max = 65533 ;
		CHLA:missing_value = 65535 ;
		CHLA:no_observation_DN = 65534 ;
		CHLA:mask_for_statistics = 351 ;
 ...
// global attributes:
		:title = "SGLI equal latitude-longitude Japan data" ;
		:id = "GC1SG1_20230101D01D_J0000_L2SG_CHLAQ_3000.nc" ;
		:date_created = "2023-01-02T05:29:04Z" ;
		:product_version = "3000" ;
		:algorithm_developer = "Japan Aerospace Exploration Agency (JAXA)" ;
		:product_level = "Level-2" ;
		:product_name = "Level-2 Japan product (JASMES)" ;
		:satellite = "Global Change Observation Mission - Climate (GCOM-C)" ;
		:sensor = "Second-generation Global Imager (SGLI)" ;
		:pixel_number = 10800 ;
		:line_number = 10400 ;
		:upper_left_latitude = 49.99875 ;
		:upper_left_longitude = 123.00125 ;
		:upper_right_latitude = 49.99875 ;
		:upper_right_longitude = 149.99875 ;
		:lower_left_longitude = 123.00125 ;
		:lower_left_latitude = 24.00125 ;
		:lower_right_latitude = 24.00125 ;
		:lower_right_longitude = 149.99875 ;
		:grid_interval = 0.0025 ;
		:grid_interval_unit = "deg" ;
		:image_projection = "EQA (sinusoidal equal area) projection from 0-deg longitude" ;

 }

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

Python3を使用してnetCDFをGeoTiFFに変換します。
入力するnetCDFによって物理量名が異なるため、「ncdump -h」の結果を参照して引数を指定してください。

$ ./Sample_nc2geotiff.py <Input> <Output> <Product>

引数1 Input :
入力となるnetCDFをフルパスで指定します。

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

引数3 Product :
入力するファイルと対応した物理量を指定します。「ncdump -h」の結果を参照してください。

2. スクリプト内容説明

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

2-1. PythonでnetCDFの読み込み

netCDFをPythonで読み込みます。「netCDF4.Dataset」を使用します。

# read netCDF4
nc = netCDF4.Dataset(INPUT, 'r')
nc.set_auto_maskandscale(False)

読み込んだncから必要な値を取得しています。
slope / offset / valid_max は物理量変換時、upper_left_longitude / upper_left_longitude / grid_interval はGeoTiFF出力時に使用します。

netCDFによって、netCDFに格納されているデータの項目名が異なることがあるので必要に応じて修正してください。
(例えばSGLI STD、VIIRSデータは「latitude」、SGLI NRTデータは「Latitude」で大文字小文字が異なる等です。)

「ncdump -h」の結果を見て以下のように取得しています。
・latitude(行)とlongitude(列)のsizeは「dimensions」に記載があるので「nc.dimensions['latitude'].size」で取得。
・読み込むdataは、variablesのCHLAとなります。2で記載したPRDを使用して「nc.variables[PRD]」で取得。
・slope / offset / valid_maxは、nc.variables[PRD]内の項目のため、「nc.variables[PRD].scale_factor」のように取得。
・upper_left_longitude / upper_left_latitude / grid_intervalは、global attributesのため、「nc.upper_left_longitude」のように取得。

### latitude, longitude
#   SGLI STD, VIIRS
dlat = nc.dimensions['latitude'].size
dlon = nc.dimensions['longitude'].size

#   SGLI NRT
#dlat = nc.dimensions['Latitude'].size
#dlon = nc.dimensions['Longitude'].size

data = nc.variables[PRD]

### max, slope, offset
#   SGLI STD, VIIRS
valid_max = nc.variables[PRD].valid_max
slope     = nc.variables[PRD].scale_factor
offset    = nc.variables[PRD].add_offset

#   SGLI NRT
#valid_max = nc.variables[PRD].Maximum_valid_DN
#slope     = nc.variables[PRD].scale_factor
#offset    = nc.variables[PRD].add_offset

### grid_interval
#   SGLI STD, VIIRS
grid_interval = nc.grid_interval

#   SGLI NRT
#grid_interval = 0.0025

### upper_left
#   SGLI STD, VIIRS (pixel center -> top left corner of the pixel)
upper_left_longitude = nc.upper_left_longitude - (grid_interval * 0.5)
upper_left_latitude  = nc.upper_left_latitude + (grid_interval * 0.5)

#   SGLI NRT
#upper_left_longitude = nc.Upper_left_longitude
#upper_left_latitude  = nc.Upper_left_latitude


# pixel size
dx = grid_interval
dy = grid_interval*-1

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

読み込んだnetCDFの値はDN値となるため、「DN * slope + offset」で物理量に変換して、配列に格納します。
また、valid_maxより大きい値は、「no_observation_DN」や「missing_value」等となるため、nanに置き換えしています。
netCDFのDN値はushortですが、物理量変換でデータ型がfloatとなりますのでデータ型も変換しています。

# value = DN * slope + offset
val = data*slope+offset
otemp      = np.zeros((dlat,dlon),dtype=np.float32)
otemp[:,:] = val[:,:]
otemp[data>valid_max] = np.nan
otemp = otemp.astype(np.float32)

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

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

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

# 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)
outRaster.SetGeoTransform((upper_left_longitude, dx, 0, upper_left_latitude, 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) (しきさいポータル)