【MATLAB】複数の地理空間ポリゴンの外郭を取得する方法

on 4月 20, 2023

続編で地理空間ポリゴンについて

以前の記事「【MATLAB】isinterior関数を使って地理空間ポリゴンの内外を判定する方法」で地理空間ポリゴンの内か外か判定する処理を紹介しました。今回はその続編ということで、複数のポリゴンの外郭を求める方法を書きたいと思います。

というのもネットで調べてもよく分からなかったので、時間を掛けて色々と組み合わせていったらできるようになったので備忘録を兼ねて書きます。

国勢調査の東京都のポリゴンデータを使用

まず、以前の記事で国勢調査の東京都の丁目レベルのポリゴンになっているシェープファイル (.shp) をreadgeotable で読み込んで、港区だけ抽出したということをやりました。今回もそこからスタートします。

%令和2年の国勢調査のシェープファイルを読み込み
tokyoTbl = readgeotable('r2ka13.shp');

% KEY_CODEカラムから「13103」で始まる港区だけを抽出
minatoTbl = tokyoTbl(startsWith(tokyoTbl.KEY_CODE, "13103"), :);Code language: Matlab (matlab)

これで港区だけのテーブルができました。

丁目 (ちょうもく)レベルのポリゴンなので、「芝公園1丁目」のような単位の複数のポリゴンがあります。このポリゴンの外郭の座標を得るのが今回のゴールです。

MATLAB 本体のpolyjoin とunion を活用

R2023a時点のMapping Toolbox ではストレートにやる方法がないので、MATLAB 本体の2つの集合の和集合を求めるunion という関数を使います。これがポリゴンの外郭を求めるのに使えます。

そしてunion の入力できるようにするため、R2021b で導入されたgeotable2table 関数を使って、ジオテーブルから普通のテーブルに変換していきます。

ジオテーブルだとShape という1つのカラムに緯度経度の座標がシェープの形式で含まれていましたが、テーブルに変換する際に、緯度と経度を別カラムにするようにします。

% tableに変換 (Latitude、Longitudeに緯度経度が格納される)
minatoLatLonTbl = geotable2table(minatoTbl, ["Latitude" "Longitude"]);Code language: Matlab (matlab)

続いてMapping Toolbox のpolyjoin 関数を使って緯度と経度を1列ずつに結合します。さらにMATLAB 本体のpolyshape 関数で2 次元の多角形を作ります。最後にunion 関数で和集合を求めます。

% 緯度経度を1列ずつに結合
[latJoin, lonJoin] = polyjoin(minatoLatLonTbl.Latitude', minatoLatLonTbl.Longitude');

% 2 次元の多角形を生成
pgon = polyshape(lonJoin, latJoin); % x/y と lat/lonで逆向き

% 和集合の計算
[polyout,shapeID,vertexID] = union(pgon);Code language: JavaScript (javascript)

地図にプロットして確認

これでpolyout のVertices カラムの中にポリゴンの外郭の緯度経度が計算されました。

実際に地図にプロットして確認しましょう。分かりやすいように元の港区の丁目のポリゴンを薄い水色で描画し、その上に外郭を金色で描画します。geoplot の色のオプションがGeopolygon の場合とLine の場合でEdgeColor だったりColor だったり違うのでご注意ください。

% 元の複数ポリゴンの描画
geoplot(minatoTbl, EdgeColor=[0.3 0.3 0.3], FaceColor="#66ccff"); % Geopolygonプロパティ
hold on

% 計算されたポリゴン外郭を描画
geoplot(polyout.Vertices(:,2), polyout.Vertices(:,1), LineWidth=2, Color="#cc9900") % Lineプロパティ
hold offCode language: Matlab (matlab)