大学時代は気象学会所属で大気データの解析をする研究だったので、この分野をMATLABでやったらどんな感じだろうということでやってみました。
GOSAT（愛称：いぶき）とは
温室効果ガス観測技術衛星「GOSAT (Greenhouse gases Observing SATellite)」は、JAXA・環境省・国立環境研究所の共同プロジェクトとして2009年に打ち上げられた衛星で、主に二酸化炭素（CO₂）とメタン（CH₄）などの温室効果ガスを観測しています。
今回扱う「L3 CH₄データ」は、GOSATが観測したメタン濃度の結果をもとに、2.5度メッシュで月平均のカラム濃度をグリッド化しているデータです。
この投稿では、HDF5形式で提供されるL3 CH₄データをMATLABで読み込み、可視化する手順を紹介します。衛星観測データに興味がある方や、MATLABでの衛星データ処理に挑戦したい方の参考になれば幸いです。
こういうのができます！
GOSAT L3 CH₄データの入手方法
GOSATのデータは、国立環境研究所が運営する公式データ提供サイトGOSAT Data Archive Service（GDAS） から無料で入手できます。
公開されているのはCO₂とCH₄のざっくり以下のデータです。
L1…生のスペクトル・輝度データ
L2…観測から導出されたCO₂・CH₄の濃度など
L3…L2データを空間・時間でグリッド化したもの
L4A…フラックス（排出・吸収量）の推定結果
L4B…L4Aを使った濃度分布のシミュレーション結果
L1データは私のような素人には扱えず、L2データも可視化が大変なので、今回はCH₄のL3データを使います。
日本の衛星だけあってデータ説明書の日本語版があるのがとても嬉しい。
なお、今はCOGデータなる便利なものがTellusで公開されているようですが、今回はGDASから入手したデータを使います。
（パッと見ただけで便利すぎてこの記事いらない気もしてきましたが…）
解凍
まずはとってきたtarファイルを解凍します。
MATLABにはuntarというtar解凍コマンドがあるのでそれを使います。
% 対象フォルダのパス
sourceFolder = "tarファイルがあるパスを入力";
% .tar または .tar.gz ファイルの一覧を取得
tarFiles = dir(fullfile(sourceFolder, "*.tar*"));
% 解凍したファイルを集めるフォルダ
outputFolder = fullfile(sourceFolder, "all_files");
if ~exist(outputFolder, "dir")
mkdir(outputFolder);
end
% 解凍ループ
for k = 1:length(tarFiles)
tarPath = fullfile(sourceFolder, tarFiles(k).name);
% 解凍実行
untar(tarPath, outputFolder);
end
複数のtarファイルをダウンロードした場合も、同じフォルダにtarファイルを集めて上記を実行すれば解凍できます。
all_filesというフォルダにSWIRL3CH4というフォルダができ、その下にHDF5(.h5)ファイルが並んでいるはずです。
データの確認＆前処理
次はh5ファイルのデータ構造を確認します。
h5infoという関数でHDF5ファイルの構造がわかるのですが、下記を実行したところデータ構造が入っているのが
info.Groups(2).Groups
の下だったのでちょっとわかりにくいかも。
% データファイルが入っているフォルダ
dataFolder = "h5ファイルがあるパスを入力";
% 読み込むファイル名
fileName = "なんでもいいのでh5ファイルを一つ選んでファイル名を入力";
% フルパスを作成
flfname = fullfile(dataFolder, fileName)
info = h5info(flfname);
上記の実行結果のinfoを掘っていくか、公式ドキュメントを確認すると、全球メタン濃度マップを描くのに必要なデータが以下にあることがわかりました。
メタン濃度データ(補正済み):/Data/mixingRatio/XCH4BiasCorrected
緯度：/Data/geolocation/latitude
経度：/Data/geolocation/longitude
データを変数に格納しましょう。
ch4 = double(h5read(flfname, "/Data/mixingRatio/XCH4BiasCorrected")); % カラム平均濃度
lat = double(h5read(flfname, "/Data/geolocation/latitude"));
lon = double(h5read(flfname, "/Data/geolocation/longitude"));
このままマップを描いてしまうと
無効値のせいでこんなマップになってしまうので前処理が必要です。
ドキュメントを読むと無効値=-9999と書かれていますが、コマンドで探すときはsummary関数が便利です。
summary(ch4(:))
NumMissing 0
Min -9999
Median -9999
Max 1.9657
Mean -5.3555e+03
Std 4.9879e+03
今回は-9999が無効値なので、-9999をNaNに変換して再度マップ化します。
ch4(ch4 == -9999) = NaN; % 無効値を NaN に変換
figure;
ax = axesm("MapProjection","robinson","Frame","on","Grid","on");
surfm(lat, lon, ch4);
colorbar;
title("GOSAT L3 CH_4 (ppmv)")
load coastlines % 海岸線データが読み込まれる
geoshow(coastlat, coastlon, "Color", "k", "LineWidth", 1) % 海岸線を重ねる
できた！
ここで重要なのはaxesm関数とsurfm関数。
axesmで投影方法などを指定してベースのマップを作り、その上にsurfm関数でlat(緯度), lon(経度)に応じたメッシュに、ch4（濃度データ)の値に応じた色を付けています。
この2つだけであれば、あとは細かい調整をするだけです。
細かい調整
文字を見やすく
タイトルとカラーバーの文字が小さいので大きくします。
% タイトル文字を太字＋サイズ16に指定
title("GOSAT L3 CH_4 (ppmv)", "FontWeight", "bold", "FontSize", 16)
cb.FontSize = 12; % 目盛りの文字の大きさを指定
カラーバーを変える
デフォルトのカラーバー（parula）は赤系がなくて淡泊な感じなのでカラーバーを "jet" に変えてみます。
ついでに、今後のためにclimでカラーバーの範囲も固定しておきます。
colormap("jet"); % カラーマップをjetにする
clim([1.80 1.96]); % カラーバーの範囲を固定
カラフルな図になりました。
どんなカラーバーが用意されているかは、こちらのチートシートを見ると手っ取り早いです。
最新のカラーバーはcolormaplistで確認できます。2025時点の最新ではnebula(青～ピンク)が追加されています。
マップの形を変更
今はaxesm関数で "robinson" （ロビンソン図法）が指定されています。
"eqdcylin"（等距円筒図法）に変えてみましょう。
figure;
ax = axesm("MapProjection","eqdcylin","Frame","on","Grid","on"); % 等距円筒図法を指定
surfm(lat, lon, ch4); % CH4濃度を描画
cb = colorbar;
title("GOSAT L3 CH_4 (ppmv)", "FontWeight", "bold", "FontSize", 16) % タイトルを設定
cb.FontSize = 12; % カラーバーのフォントサイズを12に
colormap("jet"); % カラーマップをjetにする
clim([1.80 1.96]); % カラーバーの範囲を固定
load coastlines % 海岸線データが読み込まれる
geoshow(coastlat, coastlon, "Color", "k", "LineWidth", 1) % 海岸線を重ねる
平面的な図になりました。グリッドデータはこちらのほうがわかりやすいですね。
どんな投影方法があるかは、maps関数で確認できます。
1データ可視化のまとめ
Step1. とってきたTarファイルを解凍
※Tarファイルは１つのフォルダにまとめておく
% 対象フォルダのパス
sourceFolder = "tarファイルがあるパスを入力";
% .tar または .tar.gz ファイルの一覧を取得
tarFiles = dir(fullfile(sourceFolder, "*.tar*"));
% 解凍したファイルを集めるフォルダ
outputFolder = fullfile(sourceFolder, "all_files");
if ~exist(outputFolder, "dir")
mkdir(outputFolder);
end
% 解凍ループ
for k = 1:length(tarFiles)
tarPath = fullfile(sourceFolder, tarFiles(k).name);
% 解凍実行
untar(tarPath, outputFolder);
end
Step2. データひとつを指定して可視化
以下のコードで、dataFolderにh5ファイルが入っているフォルダパスを、fileNameに可視化したいファイル名を指定して実行すれば可視化できます。
% データファイルが入っているフォルダ
dataFolder = "h5ファイルがあるパスを入力";
% 読み込むファイル名
fileName = "可視化したいh5ファイルを一つ選んでファイル名を入力";
% フルパスを作成
flfname = fullfile(dataFolder, fileName)
% データを変数に格納
ch4 = double(h5read(flfname, "/Data/mixingRatio/XCH4BiasCorrected")); % カラム平均濃度
lat = double(h5read(flfname, "/Data/geolocation/latitude"));
lon = double(h5read(flfname, "/Data/geolocation/longitude"));
ch4(ch4 == -9999) = NaN; % 無効値を NaN に変換
% マップを作成（可視化）
figure;
ax = axesm("MapProjection","eqdcylin","Frame","on","Grid","on"); % 等距円筒図法を指定
surfm(lat, lon, ch4); % CH4濃度を描画
cb = colorbar;
title("GOSAT L3 CH_4 (ppmv)", "FontWeight", "bold", "FontSize", 16) % タイトルを設定
cb.FontSize = 12; % カラーバーのフォントサイズを12に
colormap("jet"); % カラーマップをjetにする
clim([1.80 1.96]); % カラーバーの範囲を固定
load coastlines % 海岸線データが読み込まれる
geoshow(coastlat, coastlon, "Color", "k", "LineWidth", 1) % 海岸線を重ねる
１年分を一気にマップ化
Step1. とってきたTarファイルを解凍
→上記同様なので省略
Step2. １年分のデータを並べて可視化
以下のコードで変数yearに解析したい年、dataFolderにh5ファイルが入ったフォルダを指定すれば１年分のデータを可視化できます。
ファイル名が
GOSATTFTS[年][月]….h5
になっていたので、filenamesに１年分のデータ（12個）のファイル名を格納し、順次読んでいっているだけです。
プロットを並べるのにtiledlayout関数を使っています。
year = 2024; % ← ここを解析したい年に変えればOK
dataFolder = "h5ファイルがあるパスを入力"; % データファイルのフォルダパス
pattern = sprintf("GOSATTFTS%d*.h5", year);
% FileDatastoreでHDF5ファイルを管理
fds = fileDatastore(fullfile(dataFolder, pattern), ...
'ReadFcn', @readH5Data, ...
'FileExtensions', '.h5');
% 今回は12か月分(小さい)なので全ファイルを読み込み
dataAll = readall(fds);
% タイルレイアウトの準備
figure;
tiledlayout(3, 4, "TileSpacing", "tight", "Padding", "tight");
% カラースケールの固定範囲
cmin = 1.80;
cmax = 1.96;
% 海岸線の読み込み
load coastlines;
% 12か月分をプロット
for m = 1:12
nexttile;
axesm("MapProjection", "eqdcylin", "Frame", "on", "Grid", "on");
surfm(dataAll{m}.lat, dataAll{m}.lon, dataAll{m}.ch4);
title(sprintf("%02d月", m), "FontWeight", "bold", "FontSize", 14);
clim([cmin cmax]);
geoshow(coastlat, coastlon, "Color", "k", "LineWidth", 1);
end
% カラーマップとカラーバー
colormap("jet");
cb = colorbar;
cb.Layout.Tile = "east"; % 右端に配置
cb.FontSize = 12;
% 全体タイトル
sgtitle(sprintf("%d年の月別メタン濃度分布 (ppmv)", year), ...
"FontWeight", "bold", "FontSize", 16);
%% --- HDF5読み込み関数 ---
function data = readH5Data(filename)
ch4 = double(h5read(filename, "/Data/mixingRatio/XCH4BiasCorrected"));
lat = double(h5read(filename, "/Data/geolocation/latitude"));
lon = double(h5read(filename, "/Data/geolocation/longitude"));
ch4(ch4 == -9999) = NaN;
data = struct('ch4', ch4, 'lat', lat, 'lon', lon);
end
１年分マップ化できました！季節変化が見えていますね。
次回はもう少し深く解析していきたいと思います。