Ocean Remote Sensing: Practical Data Analysis 04

衛星SSTから海洋熱波を検出する

OISST と GHRSST/MUR の海面水温データを用いて、2018年2月のオーストラリア・ニュージーランド周辺における Marine Heatwave（MHW）を検出します。ここでは、P90しきい値、SST偏差、MHW日数、強度、カテゴリー、5日以上持続条件を、MATLABで段階的に計算します。

Marine HeatwaveOISSTGHRSST/MURP90ClimatologyMHW category5-day persistenceNew Zealand

この解析の中心的な問い

2018年2月のニュージーランド周辺では、衛星SSTからどのようなMHWが検出されるのか。さらに、OISSTとGHRSST/MURでは、その分布・強度・カテゴリーの見え方がどのように変わるのかを調べます。

このページの構成

1. MHWの要点 2. MATLAB基礎 3. Block 0–16 4. 考察課題

進め方：各Blockのコードは、穴埋め箇所以外も含めてそのまま掲載しています。MATLABに順番に貼り付け、________ の部分だけを自分で考えて埋めます。答えはページ上には表示していません。

1. Marine Heatwave の要点

MHWは、その場所・その季節として異常に高い水温が持続する現象です。この解析では、2003–2014年の基準期間から作成した日別ClimatologyとP90を使い、2018年2月のSSTを判定します。

MHW候補日SST > P90 となる日。これはまだ「正式なイベント」ではなく、P90超過日です。

MHW intensitySST - clim。MHW候補日のみを残して平均・最大値を計算します。

MHW categorydT = P90 - clim を基準幅として、Moderate, Strong, Severe, Extreme に分類します。

5日以上持続条件1日だけのP90超過は除外し、5日以上続いたものだけをMHWイベントとして残します。

2. MATLABで必要な基礎知識

3次元配列sst(lon, lat, time) の第3次元が日付です。2018年2月なので通常28日です。

logical配列sst > p90 は true/false の配列です。trueを足すと日数になります。

NaN処理MHWではない日をNaNにすると、平均計算や図示から除外できます。

forループ5日以上持続条件では、各格子点の時系列を取り出して連続日数を調べます。

配布ファイル

OISST_AusNZ_DJF_climP90P10_base2003_2014.mat
OISST_AusNZ_Feb2018_daily.mat

MUR_AusNZ_DJF_climP90P10_base2003_2014_coarse010.mat
MUR_AusNZ_Feb2018_daily_coarse010.mat

3. 完成版スクリプトに対応したBlock 0–16

以下は完成版スクリプトと同じ流れです。図示・保存を含む全体を掲載し、学生が考えるべき箇所だけを黄色で示しています。

実習用データのダウンロード
以下の4つのMATファイルを、このHTMLページと同じフォルダに保存してください。 MATLABスクリプトも同じフォルダに置いて実行します。

プロダクト	ファイル	内容
OISST	OISST_AusNZ_DJF_climP90P10_base2003_2014.mat	2003–2014年を基準期間として作成したDJFのclimatology, P90, P10。
OISST	OISST_AusNZ_Feb2018_daily.mat	2018年2月の日別OISST SST。オーストラリア・ニュージーランド周辺のみ切り出し済み。
GHRSST/MUR	MUR_AusNZ_DJF_climP90P10_base2003_2014_coarse010.mat	2003–2014年を基準期間として作成したDJFのclimatology, P90, P10。MURは0.10°格子に平均化済み。
GHRSST/MUR	MUR_AusNZ_Feb2018_daily_coarse010.mat	2018年2月の日別MUR SST。元の0.01°MURを0.10°格子に平均化済み。

フォルダ構成の例

04.MHW/
  orsprac04.html
  OISST_AusNZ_DJF_climP90P10_base2003_2014.mat
  OISST_AusNZ_Feb2018_daily.mat
  MUR_AusNZ_DJF_climP90P10_base2003_2014_coarse010.mat
  MUR_AusNZ_Feb2018_daily_coarse010.mat
  student_MHW_Feb2018_practice.m

Block 0

ファイル名と基本設定

OISST と MUR/GHRSST は同じスクリプトで扱えるようにします。冒頭の product だけを切り替えます。ここは穴埋めなしです。

clear; close all; clc;

%% ============================================================
% 0. File names and settings
% ============================================================

% Choose product
product = "OISST";
% product = "MUR";

if product == "OISST"
    climFile  = "OISST_AusNZ_DJF_climP90P10_base2003_2014.mat";
    dailyFile = "OISST_AusNZ_Feb2018_daily.mat";
    outDir    = "out_student_OISST_MHW_Feb2018";
    productLabel = "OISST";
    productTag   = "OISST";
elseif product == "MUR"
    climFile  = "MUR_AusNZ_DJF_climP90P10_base2003_2014_coarse010.mat";
    dailyFile = "MUR_AusNZ_Feb2018_daily_coarse010.mat";
    outDir    = "out_student_MUR_MHW_Feb2018";
    productLabel = "GHRSST/MUR";
    productTag   = "MUR";
else
    error("Unknown product. Use OISST or MUR.");
end

if ~isfolder(outDir)
    mkdir(outDir);
end

% Minimum duration for MHW event
% Hobday-style MHW event requires at least 5 consecutive days.
minDuration = 5;

% Minimum value for dT = P90 - climatology.
% This avoids unrealistically large categories where P90-clim is very small.
dThrMin = 0.3;  % degree C

% Plot range
lonlim = [100 180];
latlim = [-60 -10];

% Zoom region around New Zealand / Tasman Sea
zoomLon = [140 180];
zoomLat = [-50 -25];

fprintf("=== MHW practice: February 2018 ===\n");
fprintf("Product          : %s\n", productLabel);
fprintf("Climatology file : %s\n", climFile);
fprintf("Daily SST file   : %s\n", dailyFile);

Block 1

Climatology / P90 / P10 を読み込む

基準期間2003–2014から作成済みのClim, P90, P10を読み込みます。ここではまだMHW計算はしません。

%% ============================================================
% 1. Load climatology / P90 / P10 data
% ============================================================

if ~isfile(climFile)
    error("Cannot find climFile: %s", climFile);
end

C = load(climFile);

% In the climatology file:
%   lonR, latR : longitude and latitude
%   clim       : daily climatology
%   p90        : daily 90th percentile threshold
%   p10        : daily 10th percentile threshold
%   sd_dates   : dates for seasonal-day axis, usually Dec-Feb

lon_clim = single(C.lonR(:).');
lat_clim = single(C.latR(:));

clim_djf = single(C.clim);
p90_djf  = single(C.p90);
p10_djf  = single(C.p10);

if isfield(C, "sd_dates")
    sd_dates = C.sd_dates;
else
    % If sd_dates does not exist, assume Dec 1 to Feb 28.
    sd_dates = (datetime(2001,12,1):datetime(2002,2,28))';
end

if isfield(C, "baseStart") && isfield(C, "baseEnd")
    baseStart = C.baseStart;
    baseEnd   = C.baseEnd;
else
    baseStart = 2003;
    baseEnd   = 2014;
end

fprintf("\nLoaded climatology data.\n");
fprintf("  lon size = %d\n", numel(lon_clim));
fprintf("  lat size = %d\n", numel(lat_clim));
fprintf("  clim size = [%d %d %d]\n", size(clim_djf,1), size(clim_djf,2), size(clim_djf,3));
fprintf("  baseline = %d-%d\n", baseStart, baseEnd);

Block 2

2018年2月のDaily SSTを読み込む

学生に渡すDaily SSTを読み込みます。OISSTは0.25度、MURは0.10度へ粗視化済みです。

%% ============================================================
% 2. Load daily SST data for February 2018
% ============================================================

if ~isfile(dailyFile)
    error("Cannot find dailyFile: %s", dailyFile);
end

D = load(dailyFile);

lon = single(D.lon(:).');
lat = single(D.lat(:));
targetDates = D.targetDates;
sst = single(D.sst);

[nlon, nlat, ntime] = size(sst);

fprintf("\nLoaded daily SST data.\n");
fprintf("  product = %s\n", D.product_name);
fprintf("  lon size = %d\n", nlon);
fprintf("  lat size = %d\n", nlat);
fprintf("  time size = %d\n", ntime);
fprintf("  first date = %s\n", datestr(targetDates(1), "yyyy-mm-dd"));
fprintf("  last  date = %s\n", datestr(targetDates(end), "yyyy-mm-dd"));

if isfield(D, "spatial_resolution")
    fprintf("  spatial resolution = %s\n", D.spatial_resolution);
end

Block 3

Daily SST と基準値のグリッド一致を確認する

MHW計算では、sst, clim, p90, p10 が完全に同じ格子でなければなりません。ここは穴埋めなしです。

%% ============================================================
% 3. Check whether grids match
% ============================================================

% MHW calculation requires daily SST, climatology, P90 and P10
% to be on exactly the same grid.

if numel(lon) ~= numel(lon_clim)
    error("Longitude size differs between daily SST and climatology.");
end

if numel(lat) ~= numel(lat_clim)
    error("Latitude size differs between daily SST and climatology.");
end

lonDiff = max(abs(double(lon(:)) - double(lon_clim(:))));
latDiff = max(abs(double(lat(:)) - double(lat_clim(:))));

fprintf("\nGrid check:\n");
fprintf("  max lon difference = %.10f\n", lonDiff);
fprintf("  max lat difference = %.10f\n", latDiff);

if lonDiff > 1e-5 || latDiff > 1e-5
    error("Daily SST grid and climatology grid do not match.");
end

Block 4

2月分のClim / P90 / P10を取り出す

Clim/P90/P10ファイルはDJF 90日分なので、2018年2月のDaily SSTに対応する2月28日分だけを取り出します。

%% ============================================================
% 4. Extract February climatology / P90 / P10
% ============================================================

% The climatology file contains DJF.
% Here we extract only February because targetDates are February 2018.

idxFeb = find(month(sd_dates) == 2);

if numel(idxFeb) ~= ntime
    error("Number of February days in climatology does not match daily SST.");
end

clim = clim_djf(:,:,idxFeb);
p90  = p90_djf(:,:,idxFeb);
p10  = p10_djf(:,:,idxFeb);

fprintf("\nExtracted February thresholds.\n");
fprintf("  clim size = [%d %d %d]\n", size(clim,1), size(clim,2), size(clim,3));
fprintf("  p90  size = [%d %d %d]\n", size(p90,1),  size(p90,2),  size(p90,3));

Block 5

入力データを図で確認する

ここから図示が始まります。SST, Clim, P90の空間分布を最初に確認します。ここは完成コードです。

%% ============================================================
% 5. Quick check: mean SST, climatology and P90
% ============================================================

meanSST  = mean(sst, 3, "omitnan");
meanClim = mean(clim, 3, "omitnan");
meanP90  = mean(p90, 3, "omitnan");

coast = load("coastlines");

figure("Color","w","Position",[100 100 1500 500]);
tiledlayout(1,3,"TileSpacing","compact","Padding","compact");

nexttile;
h = imagesc(double(lon), double(lat), double(meanSST'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(meanSST')));
set(gca, "Color", [1 1 1]);
xlim(lonlim); ylim(latlim);
grid on; box on;
colorbar;
caxis([0 32]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title(productLabel + " mean SST, Feb 2018");
xlabel("Longitude"); ylabel("Latitude");

nexttile;
h = imagesc(double(lon), double(lat), double(meanClim'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(meanClim')));
set(gca, "Color", [1 1 1]);
xlim(lonlim); ylim(latlim);
grid on; box on;
colorbar;
caxis([0 32]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Mean climatology, February");
xlabel("Longitude"); ylabel("Latitude");

nexttile;
h = imagesc(double(lon), double(lat), double(meanP90'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(meanP90')));
set(gca, "Color", [1 1 1]);
xlim(lonlim); ylim(latlim);
grid on; box on;
colorbar;
caxis([0 32]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Mean P90 threshold, February");
xlabel("Longitude"); ylabel("Latitude");

sgtitle("Step 1: Input data check, " + productLabel);

exportgraphics(gcf, fullfile(outDir, "Fig01_input_check_" + productTag + "_Feb2018.png"), "Resolution", 200);

Block 6

MHW day と intensity を計算する

ここがMHW検出の中心です。SST偏差、P90超過、日数、強度を計算します。穴埋めはMHW定義に直結する行だけです。

%% ============================================================
% 6. MHW day and intensity
% ============================================================

% SST anomaly relative to climatology
anom = ________ - ________;

% MHW day:
% A grid cell is a MHW day when daily SST exceeds daily P90.
mhw_day = ________ > ________;

% Cold-spell day:
% A grid cell is a cold-spell day when daily SST is below daily P10.
cold_day = ________ < ________;

% Remove invalid data
valid = isfinite(sst) & isfinite(clim) & isfinite(p90);
mhw_day(~valid) = false;

validCold = isfinite(sst) & isfinite(clim) & isfinite(p10);
cold_day(~validCold) = false;

% Number of MHW days in February
mhw_days = sum(________, 3);

% Number of cold-spell days in February
cold_days = sum(________, 3);

% MHW intensity:
% Usually intensity is SST - climatology during MHW days.
mhw_intensity = ________;
mhw_intensity(~mhw_day) = NaN;

% Exceedance above P90:
% This shows how much SST exceeds the threshold.
mhw_excess = ________ - ________;
mhw_excess(~mhw_day) = NaN;

% Summary maps
mean_anom = mean(anom, 3, "omitnan");
max_anom  = max(anom, [], 3, "omitnan");

mean_mhw_intensity = mean(mhw_intensity, 3, "omitnan");
max_mhw_intensity  = max(mhw_intensity, [], 3, "omitnan");

mean_mhw_excess = mean(mhw_excess, 3, "omitnan");
max_mhw_excess  = max(mhw_excess, [], 3, "omitnan");

fprintf("\nMHW day calculation finished.\n");
fprintf("  Max MHW days = %.0f days\n", max(mhw_days(:), [], "omitnan"));
fprintf("  Max mean MHW intensity = %.2f degC\n", max(mean_mhw_intensity(:), [], "omitnan"));
fprintf("  Max MHW excess above P90 = %.2f degC\n", max(max_mhw_excess(:), [], "omitnan"));

Block 7

MHW日数・強度を図示する

Block 6で計算したMHW days, intensity, excessを地図化します。図示コードも全て示します。ここでは、どの変数を地図にするのかを考えるため、プロットする配列名の一部を穴埋めにします。

%% ============================================================
% 7. Figure: MHW days and intensity
% ============================================================

figure("Color","w","Position",[100 100 1500 850]);
tiledlayout(2,3,"TileSpacing","compact","Padding","compact");

nexttile;
h = imagesc(double(lon), double(lat), double(meanSST'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(meanSST')));
set(gca, "Color", [1 1 1]);
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colorbar;
caxis([8 30]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Mean SST");
xlabel("Longitude"); ylabel("Latitude");

nexttile;
h = imagesc(double(lon), double(lat), double(meanP90'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(meanP90')));
set(gca, "Color", [1 1 1]);
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colorbar;
caxis([8 30]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Mean P90");
xlabel("Longitude"); ylabel("Latitude");

nexttile;
h = imagesc(double(lon), double(lat), double(________'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(________')));
set(gca, "Color", [1 1 1]);
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colorbar;
caxis([-3 3]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Mean SST anomaly");
xlabel("Longitude"); ylabel("Latitude");

nexttile;
h = imagesc(double(lon), double(lat), double(________'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(________')));
set(gca, "Color", [1 1 1]);
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colorbar;
caxis([0 28]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("MHW days: SST > P90");
xlabel("Longitude"); ylabel("Latitude");

nexttile;
h = imagesc(double(lon), double(lat), double(________'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(________')));
set(gca, "Color", [1 1 1]);
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colorbar;
caxis([0 5]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Mean MHW intensity");
xlabel("Longitude"); ylabel("Latitude");

nexttile;
h = imagesc(double(lon), double(lat), double(________'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(________')));
set(gca, "Color", [1 1 1]);
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colorbar;
caxis([0 3]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Mean exceedance above P90");
xlabel("Longitude"); ylabel("Latitude");

sgtitle("Step 2: " + productLabel + " MHW days and intensity, February 2018");

exportgraphics(gcf, fullfile(outDir, "Fig02_MHW_days_intensity_" + productTag + "_Feb2018.png"), "Resolution", 200);

Block 8

日別MHW intensityのスナップショットを描く

日平均ではなく、特定日のMHW強度を表示します。イベントが日ごとにどう見えるかを確認します。ここでは、日別に切り出す3次元配列名を穴埋めにします。

%% ============================================================
% 8. Daily snapshots of MHW intensity
% ============================================================

snapshotDates = [datetime(2018,2,1), datetime(2018,2,10), ...
                 datetime(2018,2,20), datetime(2018,2,28)];

figure("Color","w","Position",[100 100 1400 850]);
tiledlayout(2,2,"TileSpacing","compact","Padding","compact");

for k = 1:numel(snapshotDates)

    it = find(targetDates == snapshotDates(k), 1);

    A = ________(:,:,it);

    nexttile;
    h = imagesc(double(lon), double(lat), double(A'));
    set(gca, "YDir", "normal");
    set(h, "AlphaData", isfinite(double(A')));
    set(gca, "Color", [1 1 1]);

    xlim(zoomLon);
    ylim(zoomLat);
    grid on; box on;
    colorbar;
    caxis([0 5]);

    hold on;
    plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);

    title("MHW intensity: " + string(datestr(snapshotDates(k), "yyyy-mm-dd")));
    xlabel("Longitude");
    ylabel("Latitude");
end

sgtitle("Step 3: Daily MHW intensity snapshots, " + productLabel);

exportgraphics(gcf, fullfile(outDir, "Fig03_daily_MHW_intensity_" + productTag + "_Feb2018.png"), "Resolution", 200);

Block 9

MHWカテゴリーを計算する

P90とClimの差 dT を基準に、Moderate/Strong/Severe/Extreme を計算します。ここも重要な穴埋めです。

%% ============================================================
% 9. MHW category calculation
% ============================================================

% MHW categories are defined by the distance between climatology and P90.
%
% dT = P90 - climatology
%
% Category 1: Moderate
%   clim + 1*dT <= SST < clim + 2*dT
%
% Category 2: Strong
%   clim + 2*dT <= SST < clim + 3*dT
%
% Category 3: Severe
%   clim + 3*dT <= SST < clim + 4*dT
%
% Category 4: Extreme
%   clim + 4*dT <= SST

dT = ________ - ________;

% Avoid too small dT
dT_eff = dT;
dT_eff(dT_eff < dThrMin) = dThrMin;

% Ratio of anomaly to threshold distance
ratio = (________ - ________) ./ ________;

% Initialize category array
% 0 = no MHW
% 1 = Moderate
% 2 = Strong
% 3 = Severe
% 4 = Extreme
cat_day = zeros(nlon, nlat, ntime, "uint8");

% Temporary category
tmp = floor(ratio);

% Only MHW days are categorized
tmp(~mhw_day) = 0;

% Limit category values
tmp(mhw_day & tmp < 1) = 1;
tmp(mhw_day & tmp > 4) = 4;

cat_day = uint8(tmp);

fprintf("\nMHW category calculation finished.\n");
fprintf("  Max daily category = %d\n", max(cat_day(:)));

Block 10

5日以上持続条件を適用する

1点ごとの時系列を取り出し、5日以上連続してP90を超えたイベントだけを残します。穴埋めは2箇所だけです。

%% ============================================================
% 10. Apply 5-day persistence condition
% ============================================================

% cat_day includes all days with SST > P90.
% However, a formal MHW event requires persistence.
%
% Here, only events lasting 5 consecutive days or longer are retained.
% Short gaps are not merged in this beginner-friendly version.

cat_5day = zeros(size(cat_day), "uint8");

fprintf("\nApplying %d-day persistence condition...\n", minDuration);

for ix = 1:nlon

    if mod(ix, 20) == 0
        fprintf("  lon index %d / %d\n", ix, nlon);
    end

    for iy = 1:nlat

        % Extract category time series at one grid cell
        xcat = squeeze(cat_day(ix, iy, :));

        % Convert category into true/false MHW days
        xhot = xcat > ________;

        % If there is no MHW day, skip this grid cell
        if ~any(xhot)
            continue;
        end

        % Find continuous runs of MHW days.
        dx = diff([false; xhot(:); false]);

        runStart = find(dx == 1);
        runEnd   = find(dx == -1) - 1;

        runLength = runEnd - runStart + 1;

        for irun = 1:numel(runStart)

            if runLength(irun) >= ________

                i1 = runStart(irun);
                i2 = runEnd(irun);

                cat_5day(ix, iy, i1:i2) = xcat(i1:i2);
            end
        end
    end
end

fprintf("5-day persistence condition finished.\n");
fprintf("  Max category after 5-day filter = %d\n", max(cat_5day(:)));

Block 11

カテゴリー別日数を集計する

MHW日数、カテゴリー別日数、Strong以上の日数、最大カテゴリーを計算します。ここでは、どのカテゴリーを数えているのかを意識するために、要所を穴埋めにします。

%% ============================================================
% 11. Summary maps for category
% ============================================================

% MHW days before 5-day condition
mhw_days_raw = sum(________ >= 1, 3);

% MHW days after 5-day condition
mhw_days_5day = sum(________ >= 1, 3);

% Days by each category
cat1_days = sum(cat_5day == ________, 3);
cat2_days = sum(cat_5day == ________, 3);
cat3_days = sum(cat_5day == 3, 3);
cat4_days = sum(cat_5day == 4, 3);

% Days with category 2 or higher
cat_ge2_days = sum(cat_5day >= ________, 3);

% Maximum category in February
max_cat_5day = max(________, [], 3);

fprintf("\nCategory summary:\n");
fprintf("  Max raw MHW days   = %.0f\n", max(mhw_days_raw(:), [], "omitnan"));
fprintf("  Max 5-day MHW days = %.0f\n", max(mhw_days_5day(:), [], "omitnan"));
fprintf("  Max Cat >= 2 days  = %.0f\n", max(cat_ge2_days(:), [], "omitnan"));

Block 12

カテゴリー用カラーパレットを定義する

カテゴリー図で使う色を定義します。白、黄色、オレンジ、赤、濃赤紫の順です。

%% ============================================================
% 12. Category color palette
% ============================================================

% Same idea as previous GHRSST-style MHW category plots:
% 0 = white
% 1 = yellow
% 2 = orange
% 3 = red
% 4 = dark red / purple

catCmap = [
    1.00 1.00 1.00;   % 0 no MHW
    1.00 0.90 0.20;   % 1 Moderate
    1.00 0.55 0.00;   % 2 Strong
    0.85 0.15 0.05;   % 3 Severe
    0.45 0.00 0.25    % 4 Extreme
];

Block 13

日別MHWカテゴリーを図示する

Block 9で作った日別カテゴリーをスナップショットとして描きます。ここでも、どの配列を描くべきかを穴埋めにします。

%% ============================================================
% 13. Figure: Daily MHW category snapshots
% ============================================================

figure("Color","w","Position",[100 100 1400 850]);
tiledlayout(2,2,"TileSpacing","compact","Padding","compact");

for k = 1:numel(snapshotDates)

    it = find(targetDates == snapshotDates(k), 1);

    A = ________(:,:,it);

    nexttile;
    imagesc(double(lon), double(lat), double(A'));
    set(gca, "YDir", "normal");

    xlim(zoomLon);
    ylim(zoomLat);

    colormap(gca, catCmap);
    caxis([-0.5 4.5]);

    cb = colorbar;
    cb.Ticks = 0:4;
    cb.TickLabels = {'No MHW','Moderate','Strong','Severe','Extreme'};
    ylabel(cb, "MHW category");

    hold on;
    plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);

    grid on; box on;
    xlabel("Longitude");
    ylabel("Latitude");
    title("Category: " + string(datestr(snapshotDates(k), "yyyy-mm-dd")));
end

sgtitle("Step 4: Daily MHW category, " + productLabel + ", February 2018");

exportgraphics(gcf, fullfile(outDir, "Fig04_daily_MHW_category_" + productTag + "_Feb2018.png"), "Resolution", 200);

Block 14

MHWカテゴリー要約図を作る

5日条件後のMHW日数、最大カテゴリー、Cat≥2日数、5日条件で除外された日数を描きます。どの集計量を各パネルに入れるかを考えるため、要所を穴埋めにします。

%% ============================================================
% 14. Figure: Category summary
% ============================================================

figure("Color","w","Position",[100 100 1500 850]);
tiledlayout(2,2,"TileSpacing","compact","Padding","compact");

nexttile;
h = imagesc(double(lon), double(lat), double(________'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(________')));
set(gca, "Color", [1 1 1]);
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colorbar;
caxis([0 28]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("MHW days after 5-day condition");
xlabel("Longitude"); ylabel("Latitude");

nexttile;
imagesc(double(lon), double(lat), double(________'));
set(gca, "YDir", "normal");
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colormap(gca, catCmap);
caxis([-0.5 4.5]);
cb = colorbar;
cb.Ticks = 0:4;
cb.TickLabels = {'No MHW','Moderate','Strong','Severe','Extreme'};
ylabel(cb, "Maximum category");
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Maximum MHW category");
xlabel("Longitude"); ylabel("Latitude");

nexttile;
h = imagesc(double(lon), double(lat), double(________'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(________')));
set(gca, "Color", [1 1 1]);
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colorbar;
caxis([0 28]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Days with category >= 2");
xlabel("Longitude"); ylabel("Latitude");

nexttile;
removed_by_5day = ________ - ________;

h = imagesc(double(lon), double(lat), double(________'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(________')));
set(gca, "Color", [1 1 1]);
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colorbar;
caxis([0 10]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Removed by 5-day condition");
xlabel("Longitude"); ylabel("Latitude");

sgtitle("Step 5: " + productLabel + " MHW category summary, February 2018");

exportgraphics(gcf, fullfile(outDir, "Fig05_MHW_category_summary_" + productTag + "_Feb2018.png"), "Resolution", 200);

Block 15

カテゴリー別日数を図示する

Moderate, Strong, Severe, Extremeがそれぞれ何日あったかを描きます。カテゴリー番号と変数名の対応を考えるため、プロットする配列名を穴埋めにします。

%% ============================================================
% 15. Figure: Days by category
% ============================================================

figure("Color","w","Position",[100 100 1500 850]);
tiledlayout(2,2,"TileSpacing","compact","Padding","compact");

nexttile;
h = imagesc(double(lon), double(lat), double(________'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(________')));
set(gca, "Color", [1 1 1]);
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colorbar;
caxis([0 28]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Moderate days, Cat 1");
xlabel("Longitude"); ylabel("Latitude");

nexttile;
h = imagesc(double(lon), double(lat), double(________'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(________')));
set(gca, "Color", [1 1 1]);
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colorbar;
caxis([0 28]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Strong days, Cat 2");
xlabel("Longitude"); ylabel("Latitude");

nexttile;
h = imagesc(double(lon), double(lat), double(________'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(________')));
set(gca, "Color", [1 1 1]);
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colorbar;
caxis([0 28]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Severe days, Cat 3");
xlabel("Longitude"); ylabel("Latitude");

nexttile;
h = imagesc(double(lon), double(lat), double(________'));
set(gca, "YDir", "normal");
set(h, "AlphaData", isfinite(double(________')));
set(gca, "Color", [1 1 1]);
xlim(zoomLon); ylim(zoomLat);
grid on; box on;
colorbar;
caxis([0 28]);
hold on;
plot(coast.coastlon, coast.coastlat, "k", "LineWidth", 0.8);
title("Extreme days, Cat 4");
xlabel("Longitude"); ylabel("Latitude");

sgtitle("Step 6: Number of days by MHW category, " + productLabel + ", February 2018");

exportgraphics(gcf, fullfile(outDir, "Fig06_days_by_category_" + productTag + "_Feb2018.png"), "Resolution", 200);

Block 16

結果を保存する

最後に、計算した変数をMATファイルとして保存します。ここまで実行すれば完成です。

%% ============================================================
% 16. Save results
% ============================================================

outMat = fullfile(outDir, productTag + "_MHW_Feb2018_student_result.mat");

save(outMat, ...
    "lon", "lat", "targetDates", ...
    "sst", "clim", "p90", "p10", ...
    "anom", "mhw_day", "cold_day", ...
    "mhw_days", "cold_days", ...
    "mhw_intensity", "mhw_excess", ...
    "mean_anom", "max_anom", ...
    "mean_mhw_intensity", "max_mhw_intensity", ...
    "mean_mhw_excess", "max_mhw_excess", ...
    "dT", "dT_eff", "cat_day", "cat_5day", ...
    "mhw_days_raw", "mhw_days_5day", ...
    "cat1_days", "cat2_days", "cat3_days", "cat4_days", ...
    "cat_ge2_days", "max_cat_5day", "removed_by_5day", ...
    "minDuration", "dThrMin", ...
    "baseStart", "baseEnd", ...
    "-v7.3");

fprintf("\nSaved result file:\n  %s\n", outMat);

disp("Done.");

4. 最後の考察課題

2018年2月のMHWは、どの海域に多く分布していましたか。
MHW days と Maximum MHW category は同じ意味ですか。違う場合、何が違いますか。
Days with category >= 2 が多い場所はどこですか。
Removed by 5-day condition が大きい場所は、どのような場所だと解釈できますか。
OISSTとGHRSST/MURで、MHWの分布はどこが似ていますか。
OISSTとGHRSST/MURで、MHWの分布はどこが違いますか。
その違いは、空間解像度、補間、前線構造、沿岸域、海氷・雲の影響のどれに関係しそうですか。

戻る

海洋リモートセンシング講義ページへ戻る