先輩に捧げる全国の飛び地リストを作る（市区町村編）

先輩が「飛び地を訪れる」のが趣味[^1] とのことなので、全国の飛び地を調べてみました。
「飛び地のまとめサイト」なんかも既にありますが、せっかく国土地理院さんが協力されているので、ここでは国土地理院が提供する地図データを使って飛び地を抽出してみましょう。

結果だけ見たい人は [**ココ**](https://github.com/amay077/japan_tobichi/tree/master/out) の適当なファイルをクリックしてみてください。ファイル名の連番は [都道府県コード](http://nlftp.mlit.go.jp/ksj/gml/codelist/PrefCd.html) です。

[^1]: テレビで先輩が訪れていた地域は厳密?には飛び地でなく「細長く入り組んだ境界」でしたが、[FUZ](https://comic-fuz.com/series/1856) で観たらちゃんと「飛び地」でしたね。

## 使用する地図データと対象とする行政区画

国土地理院が提供する「地球地図日本」の「行政界」データを使います。

* [地球地図日本｜国土地理院](https://www.gsi.go.jp/kankyochiri/gm_jpn.html)

このデータに含まれる行政界は「市区町村」なのでそれを対象とします（町丁目レベルの細かい飛び地は対応できません）。
また、データの精度は高くないので、市区町村レベルでは微小とされる細かな飛び地は拾えないものと思います。

1. この地球地図日本のサイトから、「第2.2版ベクタ（2016年公開）／全レイヤ」のファイルをダウンロードして解凍します。
2. 解凍してできたファイル群の内、``polbnda_jpn.shp``, ``polbnda_jpn.shx``, ``polbnda_jpn.dbf``, ``polbnda_jpn.prj`` の 4つ のファイルをZIP圧縮して、ファイル名を ``polbnda_jpn.zip`` とします（polbnda_jpn.xxx が「行政界」のポリゴンデータ(Polygon Boundary Japan の略と思われる)なので、他のは不要です）。

この ``polbnda_jpn.zip`` というファイル、 https://mapshaper.org/ というサイトを開いて、ドラッグ＆ドロップすると、その内容が見られます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/c27be93d-74f7-77d9-6f1b-73cc67246373.png)

## MapShaper で飛び地を確認

このデータが本当に使用に耐えられるものか、引き続き MapShaper で確認してみましょう。

例えば「千葉県船橋市」。この市は大きな「本体ポリゴン」の他に、3つ の「サブポリゴン」を持っています。

<table>
<thead>
<tr>
<th align="center">本体ポリゴン</th>
<th align="center">サブポリゴン</th>
</tr>
</thead>
<tbody>
<tr>
<td align="center"><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/8beadd40-0219-eb7b-badc-1b9b23b0af53.png" style="max-width:100%;"></td>
<td align="center"><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/8b9dde48-5275-85b0-8db5-fb2d88730b86.png" style="max-width:100%;"></td>
</tr>
</tbody>
</table>

左上に表示される属性情報を見てみると次のような特徴がありそうです。

* nam: 都道府県名
* laa: 市区町村名
* adm_code: 市区町村コード
* pop: 人口(ただしサブポリゴンの場合は、``-89999999`` が入る

以上を踏まえると、

* 本体ポリゴン: ``pop`` の値が 0以上
* サブポリゴン: ``pop`` の値が 0未満

という雑な条件で判定できそうです。

尚、このデータの詳細な仕様は、[Global Map Specifications](https://github.com/globalmaps/specifications/blob/master/README.md) に置いてある PDF で確認できます。

``pop`` については、

> -99999999 Unknown 
> -89999999 *1 Attributed in the main polygon

と書かれているのでまあよいでしょう。

## 「飛び地」の定義

さて、本体ポリゴンとサブポリゴン群を識別することができましたが、ところで「飛び地」とは何でしょう？
幾何学形状だけを見れば、陸地続きの離れた地域も、「島」も飛び地とみなせます。
また、他の行政界領域に空いた「穴」のみを飛び地とみなしたい人も居るかもですね。

ここでは、先輩が訪れやすいように、「島」は除外して「陸地続きの離れた地域」を「飛び地」とすることにします。
~~上の船橋市の例では南にある「島」のサブポリゴンは飛び地とみなしません。~~

処理に不具合があったので、「サブポリゴンと隣接するサブポリゴン」は飛び地とする事にしました。たぶん島の多い地域ではカウントが増えます。

データの下調べは十分に行えたので、このデータを使って飛び地を抽出するプログラムを作成します。

## 必要な環境など

### 環境

* [node.js](https://nodejs.org/ja/): v12.14.1 以上
* [TypeScript](https://www.typescriptlang.org/): v3.7.3

TypeScript/JavaScript を使うので Windows でも macOS でも Linux でも OK なはずです。

### 主な使用ライブラリ

* [shpjs](https://www.npmjs.com/package/shpjs): Shapefile(準備した ``polbnda_jpn.zip``) を読み込んで GeoJson 形式に変換する JavaScript ライブラリ
* [turfjs](https://www.npmjs.com/package/@turf/turf): 地理空間データの接合判定とか距離計算とかいろいろできちゃうすごい JavaScript ライブラリ

これらがあれば8割くらいは出来たも同然です。

## 実装

そしてできたものがこちらになります。

* [amay077/japan_tobichi: 国土地理院の地球地図日本データから飛び地を抽出するツール](https://github.com/amay077/japan_tobichi)

### 使い方

1. [地球地図日本データ](https://www.gsi.go.jp/kankyochiri/gm_jpn.html) 第2.2版ベクタ（2016年公開）の「全レイヤ」をダウンロードする
2. ファイルを解凍して ``polbnda_jpn.*`` だけを圧縮して ``polbnda_jpn.zip`` とする
3. このリポジトリを Clone する
4. ``polbnda_jpn.zip`` を ``tool/assets`` ディレクトリに置く
5. ``tool`` ディレクトリで ``npm ci`` を行う
6. ``npm run exec`` を行う
7. ・・・しばらく待つ
8. ``../out`` ディレクトリに ``*.geojson`` ファイルが出力される

## 実装の簡単な解説

すべて ``src/tobichi-extractor.ts`` で行っています。

### 1. Shapefile の読み込み

```typescript
// import * as shp from 'shpjs';

// Shapefile を読み込み
const buf = fs.readFileSync(Path.resolve(__dirname, '../assets/polbnda_jpn.zip'));
const geoJson = await shp(buf);

const sourceFeatures: Feature[] = geoJson.features
```

``shpjs`` を使って ``polbnda_jpn.zip`` を読み込みます。ファイルパスを直接指定できなかったので、自力でバイナリとして読んでそのデータを渡します。

読み込まれた ``geoJson`` は GeoJson 形式の FeatureCollection です。


### 2. 本体ポリゴンとサブポリゴンの識別

```typescript
const sourceFeatures: Feature[] = geoJson.features
  .filter(f => f.properties.pop != -99999999) // UNKNOWN は除外
  .map((f, index) => {
    f.properties._no = index; 
    return f;
  });

const mainFeatures = sourceFeatures.filter(f => f.properties.pop >= 0);
const subFeatures = sourceFeatures.filter(f => f.properties.pop < 0);
```

``pop`` が ``UNKNOWN`` なデータを除外し、Feature を識別できるように ``_no`` に採番しておきます。

「2. 本体ポリゴンとサブポリゴンの識別」で定義したように、``pop`` 値で本体かサブを判定して、それぞれ別々のリストにしておきます。

### 3. サブポリゴン群から「飛び地」とみなすポリゴンのみを抽出

```typescript
// import * as turf from '@turf/turf';

const tobichiSubFeatures = subFeatures.filter(subF => {
  return sourceFeatures.find(outerF => {
    // 自分は除外
    if (outerF.properties._no === subF.properties._no) {
      return false;
    }
    return !turf.booleanDisjoint(outerF, subF);
  }) != null;
});
```

~~本体ポリゴン群と「接している」サブポリゴンは「飛び地」とみなしています。~~
他のポリゴン群と「接している」サブポリゴンは「飛び地」とみなしています。

「接している」の判定は [``turf.booleanDisjoint``](http://turfjs.org/docs/#booleanDisjoint) を使います。これは「図形と図形が離れているか？」を調べる関数ですが、これが ``false`` であれば「接している」とみなせます(他によい関数があるかも知れません)。

~~また、本体ポリゴン群のみ判定に使っているのは、冒頭で表示した千葉県船橋市の南の島のような、「他のサブポリゴンに接している島」を除外するためです。~~
↑この条件は、青森県中泊町のように本州の端っこでサブポリゴンにしか接していないサブポリゴンが意図せず除外されてしまうのでやめました。

本体ポリゴンがサブポリゴンを内包している(本体ポリゴンの一部が他市区町村の飛び地である)場合、``turf.disjoint`` が ``true`` を返すのかは未検証ですなぜなら、このような場合、本体ポリゴンには「穴」が適切に空いているためです[^2]。``turf.disjoint`` は穴も適切に判定しているようです。

[^2]: GeoJson の Polygon は「外周」「穴(群)」の座標群を持つ事ができ、外周は座標が時計回り、「穴(群)」は反時計周りとなっています。

そしてこの「総当りの」二重ループはとても遅いですが、我慢できる範囲なので今は最適化はしません。

ここで得られた ``tobichiSubFeatures`` が、全国市区町村の「飛び地」のリストです。
あとはこれを本体ポリゴンと組にして GeoJson ファイルに出力するのが、後続の処理ですが、ここでは説明を割愛します。

## 飛び地を地図で見てみよう

処理を実行した結果の ``.geojson`` ファイル群は、``tool`` と同じレベルの ``out`` というディレクトリに出力されます。

* ``tobichi_00_all.geojson`` … すべての飛び地とその本体ポリゴンが格納された GeoJson ファイル
* ``tobichi_{nn}.geojson`` … 都道府県ごとの飛び地とその本体ポリゴンが格納された GeoJson ファイル

このファイルを地図上に表示する方法を 2つ 紹介します。

### GitHub で見る

GitHub は ``.geojson`` ファイルを自動的に Mapbox の地図に描画してくれます。

私が出力してコミットしておいた https://github.com/amay077/japan_tobichi/tree/master/out の任意のファイルをクリックすると、それぞれ地図上に表示ができます。

下図は千葉県の ``tobichi_12.geojson`` の船橋市付近を表示してみた例です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/599bdfbf-4cee-08b7-1a9c-6e765e3f3996.png)

### 地理院地図で見る

最後まで国土地理院さんにお世話になりましょう。
国土地理院が提供する「地理院地図」というサービスがあります。

* [地理院地図｜国土地理院](https://maps.gsi.go.jp/)

ブラウザで上記URLを開いて、ローカルの任意の ``tobichi_xx.geojson`` ファイルを、表示されている地図へドラッグ＆ドロップします。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/8161e20e-f3b5-1814-11da-199dbb75626f.png)

すると上図のように ``geojson`` ファイルが地図に描画されます。
``tobichi_00_all.geojson`` は少し負荷が高いので、都道府県毎のファイルの方が快適です。
こちらも飛び地ポリゴンと本体ポリゴンが描画され、同じ市区町村は同じ色になります。

複数の ``.geojson`` ファイルをドロップして表示できるし、右上の「機能」→ ツール → 作図・ファイル で、表示の ON/OFF や図形やスタイルの変更もできます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/15b32384-db02-f733-8ae7-101a2cf119af.png)

背景地図がごちゃごちゃして見づらいならば、右上の「情報」→ ベースマップ で、「白地図」や「淡色地図」を選んでもよいし、ベースマップを選んでから、グレースケールON や 透過率を設定することもできます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/e1ef798a-6533-a870-0699-8eab41d32b32.png)

衛星(航空)写真や、先輩と同じように国土地理院の地図が見たい人は、地理院地図を使うのがオススメです。

## Statistics

抽出したデータをいろいろ集計してみました。

お断り：地図データの精度により、飛び地の数や面積、距離など、実際とは異なりますので予めご了承ください。

飛び地ポリゴンの判定処理を変更したのと件数に本体ポリゴンが含まれていたので再集計しました。

### 飛び地の数

まず、全国市区町村レベルでの飛び地の総数は **144** 箇所でした。

都道府県毎で、最も飛び地を多く持つのは、

1. 27: 大阪府(11箇所)
2. 46: 鹿児島県(10箇所)
3. 06: 山形県(9箇所)
4. 12: 千葉件(9箇所)

でした。

再集計したら大阪府が1位になってしまった。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/f7580cd2-49a2-ce36-2222-d4a55cc9f07d.png)

関西国際空港って、３つに割れとったんかい！

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/883c15d0-d47f-577f-3c2a-ea30d60c948b.png)


3位の山形県を地理院地図で見てみましょう。一見するとそう飛び地が多くは見えませんが。。。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/8dd99652-63e1-169a-d928-93ca211db38e.png)

拡大してみると・・・

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/b9996916-7a00-e9fb-74de-edeccd879459.png)

山形県大江町、なかなか闇の深そうな飛び地をしておられます。写真で見ると田畑っぽいですね。。。
ちなみに市区町村毎の飛び地数でも山形県大江町が1位(6つ)でした。

一方、飛び地がひとつも無い都道府県が 10 ありました。領地争いの少ない、平和な土地だったのでしょうか。。

1. 03: 岩手県
2. 22: 静岡県
3. 25: 滋賀県
4. 29: 奈良県
5. 31: 鳥取県
6. 38: 愛媛県
7. 40: 福岡県
8. 41: 佐賀県
9. 44: 大分県
10. 45: 宮崎県

東京都もお台場などが飛び地扱いになりますね。練馬区の例の場所は精度上拾えませんが。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/6cd78b1c-e87c-a66c-d418-0d49ffc171ad.png)

### 面積の大きい／小さい飛び地

次に、「面積が最も小さい」飛び地は、宮城県仙台市 の飛び地で **2,833㎡** でした。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/f7ad7fa3-597e-b7ed-b84f-f8592f84d8c3.png)

これは小っさ！ちなみに地理院地図では距離や面積の計測もできますよ。

反対に「面積が最も大きな」飛び地は、やはり北海道、 北海道日高町 で、面積は **564,806,056㎡** でした。ケタ違い。
尚、このランキングは上位3位を北海道が独占しました。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/b6e8bff3-dfc6-b69c-31ac-f1d510bb48bf.png)

地図で見ると本体よりデカいんじゃないかと思えます。お隣の伊達市の飛び地もそうですね。

という事で次は「飛び地の方が本体よりデカい市区町村」を調べてみました。

**7件** もありました。

1. 01233:北海道伊達市 - 本:175,055,053㎡ < 飛:274,726,478㎡
2. 01601:北海道日高町 - 本:429,338,399㎡ < 飛:564,806,056㎡
3. 10203:群馬県桐生市 - 本:137,068,936㎡ < 飛:138,857,m922㎡
4. 20385:長野県南箕輪村 - 本:20,433,376㎡ < 飛:20,503,493㎡
5. 21202:岐阜県大垣市 - 本:79,689,605㎡ < 飛:123,817,042㎡
6. 27362:大阪府田尻町 - 本:2,479,536㎡ < 飛:3,343,210㎡
7. 43212:熊本県上天草市 - 本:29,751,514㎡ < 飛:80,488,305㎡

北海道の２つはやはりとして、熊本県上天草市は本体小さスギィ！

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/2dedd6b0-de5f-ca20-5e03-e41ad46e2062.png)

地図で見ると、むしろ **本体が島** なのですね。

当初は、「都道府県ポリゴン群で一番面積が広いのが本体」などと思案していたので、これはよい発見でした、今後に活きます。

### 本体から飛び地までの距離

最後に、「本体と飛び地の距離」を調べてみました。

まず、距離の遠いランキング。

1. 01601:北海道日高町 - 16km [^4]
2. 01206:北海道釧路市 - 12km
3. 34211:広島県大竹市 - 11km

[^4]: [日本一離れた飛び地「北海道沙流郡日高町」 | ココロココ 地方と都市をつなぐ・つたえる](https://cocolococo.jp/27311)

北海道勢が安定のワン・ツー。
3位の 広島県大竹市 を見てみましょう。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/e16f18e3-2756-1db6-3aa2-e31fe9c6af26.png)

「甲島」という島が、大竹市と岩国市で分け合う形になってますね。どうしてこうなった。

ちなみに距離の算出は「頂点 vs 頂点」で行ってます。嘘、[``turf.nearestPointOnLine``](http://turfjs.org/docs/#nearestPointOnLine) で行っており飛び地側は「線上の点」です。

逆に、距離の近いランキング。

1. 08546:茨城県境町 - 0km
2. 11100:埼玉県さいたま市 - 0km
3. 14100:神奈川県横浜市 - 0.027km

おい 0km ってｗｗ
プログラムのバグかと思いましたが一応見てみます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/eb6cec48-e164-6832-0c43-49f9431eb050.png)

んーなるほど！これはデータの精度のせいで、実際には離れていますね。ベースマップを見ると分かります。「境町飛地」とも書いてあるし。

つかこの飛び地、すごい豪邸に見えるのだけど、「ポツンと豪邸」じゃない？

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/0c2a0f0a-8f7a-002c-a3e7-035eb91e442c.png)

### おまけ１

広島県広島市に穴が空いてるのに、この穴は隣の市の飛び地じゃないの？

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/3b772c64-6562-efbb-8f37-a4299c90ac7f.png)

と思ったら、[広島県安芸郡府中町](https://www.town.fuchu.hiroshima.jp/) はこの **穴でしかない** のですね。四面楚歌じゃん。

### おまけ２

山梨県の富士山付近に、小さな飛び地があります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/499fc56e-fdd8-786f-bdb7-a18a72dfd328.png)

富士河口湖町の飛び地のようです。富士山へ向かう道の途中にも見えます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/8b3bfaad-85a0-bd65-5b7f-732edf2d1972.png)

拡大してみると、「小室浅間神社（おむろせんげんじんじゃ）」という神社があるようです。

この神社だけが富士河口湖町なのか？と思いましたが、[神社のホームページ](http://www.fgo.jp/~yabusame/)を見ると 富士吉田市 と書いてあるし、飛び地とは関係ないのかも知れません、謎です。。。

## まとめ

国土地理院の地図データを使って、プログラムで「飛び地」を抽出してみました。

飛び地の発生理由はさまざまですが、大きな要因の一つは「平成の大合併」と呼ばれる 2005年頃 の市区町村合併ラッシュです。

その前後の様子が↓の資料で確認できます。

* [全国都道府県別・市町村合併新旧一覧図(平成15年以降)｜国土地理院](https://www.gsi.go.jp/KOKUJYOHO/gappei_index.htm)

飛び地となっている地域と併せて、この資料を見てみるとなかなか面白いです。

今回は市区町村レベルで行いましたが、もっと詳細な地図データを使えば、町丁目界での抽出も可能かなと思います。ただしデータ量がバクハツするので線形に処理している箇所で空間インデックスを使うなどの工夫が必要です。

### より詳細な境界データ

#### 国土数値情報

詳細なデータとしては、国土地理院も所属する国土交通省が [国土数値情報 行政区域データ](http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_3.html) というデータを提供しています。

#### 国土地理院ベクトルタイル？

あるいは、地理院地図のベクトル地図版およびベクトルタイルデータが配信実験中です。

* [地理院地図｜ベクトルタイルとその提供実験について](https://maps.gsi.go.jp/development/vt.html)

詳細な行政界のベクトルタイルが配信されたなら、そこから町丁目界のポリゴンを取得することが可能になるかも知れません。[^3]

#### 統計GISの境界データ

[統計GISデータダウンロード | 政府統計の総合窓口](https://www.e-stat.go.jp/gis/statmap-search?page=1&type=2&aggregateUnitForBoundary=A&toukeiCode=00200521&toukeiYear=2015&serveyId=A002005212015&coordsys=1&format=shape) では、かなり精度の高いの町丁目界のポリゴンデータがダウンロードできるようです。こちらも国土地理院の提供で、細かな飛び地を拾える可能性があります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/e987c4e9-0330-f3a0-2026-2e979824df62.png)

また、同じデータが [H27国勢調査（小地域）境界データ - データセット - NIAES VIC](https://niaesvic.dc.affrc.go.jp/dataset/h27-census-polygon) で **GeoJSON 変換済** で提供されています。[情報ありがとうございます！](https://twitter.com/wata909/status/1222777564558548992)


[^3]: ベクトルタイルとは言え分析用ではなく「表示用」と思われるため、面(ポリゴン)データは提供されないかも知れません。

私は、実際に街歩きをするよりも、このようなデータをこねくり回す方が好きなので、この方面からお手伝いできたらよいなと思います。

## バグってた→直した

### 青森県中泊町

[青森県](https://github.com/amay077/japan_tobichi/blob/4bcc49778d0e9c0f2c9e815663aee40e003f74fc/out/tobichi_02.geojson) 中泊町は飛び地のはずですよね。。。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/60459310-6cca-54bf-dba9-8276f7ef7b7b.png)

飛び地判定ロジックを直したので正しく [拾えました](https://github.com/amay077/japan_tobichi/blob/e0279eab4ee787bdad6509441b29ee46be87d1bd/out/tobichi_02.geojson)。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/8227/2a3c63af-2b29-913d-418b-41ef13f9c09f.png)