陸域をmaskする方法

地球の格子データでu（2250×2001 double）, v（2250×2001 double）の流速のgrid data（uv.mat）に、陸域をmaskしたいです。

緯度（lat; 1×2001 double）、経度（lon; 1×2250 double）は下記です。

lon=0:0.08:179.92;lat=-80:0.08:80;

陸域にも微妙に数値データが入っていることから、陸域のu,vは全てNaNにしたいです。

しかし下記だと、u,vにうまくNaNが反映されていないようです。

どこが変なのかわからず困っています。ご検討、よろしくお願いいたします。

load(['uv.mat']);
% 陸地のみ0にする
landmask = shaperead('landareas.shp'); % landareas.shpを読み込む（MATLABのマッピングツールボックスが必要）
landlat = [];
landlon = [];
for i = 1:length(landmask)
    % landmask構造体から緯度経度を抽出
    if isfield(landmask,'LatLon')
        latlon = landmask(i).LatLon;
    else
        latlon = [landmask(i).X, landmask(i).Y];
    end
    % 緯度と経度を分ける
    lat = latlon(:, 2);
    lon = latlon(:, 1);
    % 経度が-180度から180度の範囲になるように調整
    lon(lon > 180) = lon(lon > 180) - 360;
    % 緯度と経度を陸地のデータに追加
    landlat = [landlat; lat];
    landlon = [landlon; lon];
end
[in, on] = inpolygon(lon, lat, landlon, landlat); % 陸地にあるグリッドの場所を検索
u(in) = NaN; % 陸地にあるグリッドの値を0にする
v(in) = NaN; % 陸地にあるグリッドの値を0にする

[Kojiro Saito]

上記コードでinpolygonの実行結果のinが1x1のサイズの数値(1)になってしまうのが原因だと思われます。

またバージョンはR2021a固定でしょうか？というのもR2021bでreadgeotableが導入され、シェープファイルをダイレクトにテーブルとして読み取れるようになり、R2022aで出たisinterior関数によってポリゴンの中に入るかの計算ができるようになっているため、新しいバージョンだと簡便にできるためです。

[Hiroki Takeda]

申し訳ありません。R2022aでした。

そちらの方法を教えていただくことは可能でしょうか。

よろしくお願いいたします。

Answer 1

[Kojiro Saito]

isinterior関数 (Mapping Toolbox)を使用して、陸域のポリゴンの内か外かを判定し、ポリゴン内の流速データをNaNに変換すれば実現できます。

以下がコード例です。流速データのuv.matが手元に無いので、ドキュメントに含まれている風速データのwind.matを使用しています。

R2022aでも動作確認済みです。

% 風速データの読み込み
% https://jp.mathworks.com/help/matlab/import_export/matlab-example-data-sets.html
% wind.mat 北アメリカの空気の流れに関する 3 次元データ。
% データは (x,y,z) の位置成分と (u,v,w) の速度成分で構成されています。
load("wind","x","y","u","v")
lat = y(:,:,1);
lon = -x(:,:,1); % ドキュメントの説明にはないが、北米のデータで西経表示のようなのでマイナスを付けて東経表示にする
dlat = v(:,:,1);
dlon = u(:,:,1);
clear x y u v

% マスクする陸域のポリゴンの読み込み
tbl = readgeotable('landareas.shp');

% 風速の緯度・経度からgeopointshapeの作成
shape = geopointshape(lat, lon);

% 陸域の内外判定をして陸域内のデータをNaNにする
for n=1:length(tbl.Shape)
    inpoly = isinterior(tbl.Shape(n), shape);
    dlat(inpoly) = NaN;
    dlon(inpoly) = NaN;
end

% 陸域のポリゴンをプロット
axesm('mercator', 'MapLatLimit', [min(lat, [], 'all') max(lat, [],'all')],'MapLonLimit', [min(lon,[],'all') max(lon,[],'all')],'MeridianLabel', 'on','ParallelLabel','on','MLineLocation',10,'PLineLocation',10,'MLabelLocation',10,'PLabelLocation',10)
geoshow(tbl)
hold on
% 風速をプロット
quiverm(lat,lon,dlat,dlon)
hold off

[Hiroki Takeda]

ありがとうございます。

はい、イメージ通りのご回答で、できそうです。

私の持っているuv.matは下記のような変数を持っています。

lat 1x2001 double(lat=-80:0.08:80)

lon 1x2250 double(lon=0:0.08:179.92)

u 2250x2001 double

v 2250x2001 double

上記のような処理を行うと、matlabが落ちてしまいます（計算処理能力？）。

実質は日本周辺だけが欲しいので、緯度経度を限定し、下記のように操作をしています。

clear; close;
load('uv.mat');
tbl = readgeotable('landareas.shp');
[Lat, Lon] = meshgrid(lat, lon);
shape = geopointshape(Lat, Lon);
% 緯度と経度が指定された範囲内かどうかを判定
lat_mask = (Lat >= 30) & (Lat <= 45);
lon_mask = (Lon >= 125) & (Lon <= 148);
% 範囲外の点をNaNにする
dlat(~lat_mask | ~lon_mask) = NaN;
dlon(~lat_mask | ~lon_mask) = NaN;
% 陸域の内外判定をして陸域内のデータをNaNにする
for n=1:length(tbl.Shape)
    inpoly = isinterior(tbl.Shape(n), shape);
    dlat(inpoly) = NaN;
    dlon(inpoly) = NaN;
end
filename = ['uv_mask.mat'];
save(filename,'lon','lat','dlat','dlon');

この処理を軽くする方法はございますでしょうか（上記の方法はそもそも合っていますでしょうか）。

ご検討、どうぞよろしくお願いいたします。

[Kojiro Saito]

確かに全世界のデータをそのまま適用すると、isinteriorのところで「要求された 1323x4502250 (44.4GB) 配列は、最大配列サイズの基本設定 (31.8GB) を超えています。これにより、MATLAB は反応しなくなる可能性があります。」と出てしまい、メモリーエラーでこけてしまいます。

データが0.08度刻みでメッシュが細かいので、日本周辺の緯度経度だけに絞って処理しないといけないですね。

下記が変更したサンプルです。uとvはダミーデータを入れています。

lat=-80:0.08:80;
lon=0:0.08:179.92;
u = ones(2250, 2001); % 2250x2001 double
v = 2*ones(2250, 2001); % 2250x2001 double

latMin = 30;
latMax = 46;
lonMin = 125;
lonMax = 148;
tbl = readgeotable('landareas.shp');
[Lat, Lon] = meshgrid(latMin:0.08:latMax, lonMin:0.08:lonMax);
shape = geopointshape(Lat, Lon);

% 緯度と経度が指定された範囲内かどうかを判定
lat_mask = (Lat >= latMin) & (Lat <= latMax);
lon_mask = (Lon >= lonMin) & (Lon <= lonMax);

% 範囲内のデータのみ抽出
u = reshape(u(lat_mask & lon_mask), size(Lat));
v = reshape(v(lat_mask & lon_mask), size(Lat));

% 範囲内のポリゴンのみ抽出。geoclipはR2022aから使用可能
clipped = geoclip(tbl.Shape, [latMin latMax], [lonMin lonMax]); 
idxLand = clipped.NumRegions ~= 0;
tbl = tbl(idxLand,:);

% 陸域の内外判定をして陸域内のデータをNaNにする
for n=1:length(tbl.Shape)
    inpoly = isinterior(tbl.Shape(n), shape);
    u(inpoly) = NaN;
    v(inpoly) = NaN;
end

filename = ['uv_mask.mat'];
save(filename,'lon','lat','u','v');

% 地図にプロット
axesm('mercator', 'MapLatLimit', [latMin latMax],'MapLonLimit', [lonMin lonMax])
geoshow(tbl)
hold on
quiverm(Lat, Lon, u, v)
hold off

このような結果になります。