Mapping Toolbox のisinterior 関数は便利
【MATLAB】政府統計の国勢調査のシェープファイルを読み込んで可視化の記事でMATLAB での地図データのプロットについて書きました。今回はその続編ということで、ある位置座標が地理空間ポリゴンの中にあるか外にあるかを判定する方法について紹介します。
Mapping Toolbox のisinterior関数を使うとかなり楽に実現できます。
MATLAB 本体の標準関数にもisinterior がありますが、それとは違いR2022aのMapping Toolboxで導入されたisinterior 関数なのでご留意ください。
港区のポリゴンデータを準備
例として、ある緯度経度の場所が東京都の港区の地域のポリゴンの内外にあるか判定する方法をご紹介したいと思います。
東京・港区を含むポリゴンをe-Stat からダウンロード
まずは、港区のポリゴンを入手するところからですね。前回の記事と同様に、日本政府の国勢調査のデータを使います。e-Stat のサイトから右下のランキングにある「国勢調査」をクリックして色々と選択すると出てきますが、直リンクはこちらにある、2020年の国勢調査のデータ(2023/01/04公開)の小地域 町丁・字等の「13000 東京都全域」のデータを使います。zip ファイルを解凍するとr2ka13.dbf、r2ka13.prj、r2ka13.shp、r2ka13.shxのファイルがあります。
シェープファイルの読み取りも、前回記事と同じMapping Toolbox のreadgeotable 関数を使います。
tokyoTbl = readgeotable('r2ka13.shp');Code language: Matlab (matlab)港区だけのテーブルを抽出
r2ka13.shpには東京都全域のポリゴンが入っているので、港区だけを抽出します。読み込んだテーブルにはKEY_CODEという列があり、ここが「13103xx」となっているのが港区です。13は東京都の都道府県番号、103が港区のコード、です。
ここではMATLAB 本体の文字列を操作するstartsWith 関数を使います。
minatoTbl = tokyoTbl(startsWith(tokyoTbl.KEY_CODE, "13103"), :);Code language: Matlab (matlab)これでtokyoTblの中からKEY_CODEが「13103…」のものを抽出できました。
東京都内の緯度経度をランダムに生成
続いて乱数を使ってランダムに東京都内の緯度経度のポイントを100個作ります。
latMin = 35.65870207047825;
latMax = 35.72760939676349;
lonMin = 139.67907277760844;
lonMax = 139.77870831937292;
latR = (latMax-latMin).*rand(100,1) + latMin;
lonR = (lonMax-lonMin).*rand(100,1) + lonMin;
outT = table;
outT.Shape = geopointshape(latR, lonR);Code language: Matlab (matlab)最後にoutT というテーブルを作ったのは、港区の中かどうかの判定を入れるための変数です。こちらにShape という列を作り、乱数で生成した緯度経度をgeopointshape に変換して格納しています。
港区内にあるかどうかをisinterior で判定
それではいよいよisinterior 関数を使ってポリゴンの内にあるかどうかを判定していきます。isinterior は1つ目の引数にポリゴン、2つ目の引数に内外の判定を調べたい点、を入れます。minatoTbl の各行には港区の丁目のデータが入っていますので1行ごとにポリゴンを処理していきます。
点の位置がポリゴンの中に含まれればisinterior のリターンはtrue (1) になりますので、そこを Region 列に「港区」というラベルを書いていきます。逆にポリゴンの外と判定されるとRegion の値が「missing (欠損値)」になるので、「それ以外」というラベルに変えます。
for n=1:height(minatoTbl)
inpoly = isinterior(minatoTbl.Shape(n), outT.Shape);
outT.Region(inpoly) = "港区";
end
outT.Region(ismissing(outT.Region)) = "それ以外";
outT.Region = categorical(outT.Region);Code language: Matlab (matlab)地図にプロットして確認
これで内外の判定ができましたが、実際に地図にプロットして確認してみましょう。港区の位置が分かるようにまず港区のポリゴンをgeoplot で描き、その上に港区内の点と、区外の点をoutT.Region の値で色を分けてプロットします。
geoplot(minatoTbl);
hold on
geoscatter(outT.Shape.Latitude, outT.Shape.Longitude, [], outT.Region);
hold offCode language: Matlab (matlab)実行結果はこのようになります。港区の中の点が黄色、区外の点が紫色でプロットされていて、ポリゴンの内外を判定できていることが分かりました。
Mapping Toolbox のisinterior 関数を使ったポリゴンの領域内外の判定方法をご紹介しました。
No responses yet