Elixir Livebook で国土交通省の行政区域データを可視化する（都道府県・市区町村の形を地図上にプロットする）

## はじめに

みなさん、都道府県・市区町村の位置と形、分かりますか？

我が家では日本地図パズルで娘とよく勉強したものです

https://www.kumonshuppan.com/kumontoy/kumontoy-syousai/?code=54496

都道府県・市区町村の形が分かると何が良いのか、、、

都道府県・市区町村の地理情報=座標=緯度経度があれば、ある地点が県や市の中なのか外なのかが分かるということです

例えば降水量、土壌分析、植生、夜間の光量など、位置と紐付いている情報を都道府県や市区町村の形に切り取ることで、各県内や各市内の最小・最大・平均など基礎統計量を取得することができます

それを使えば県や市を数値的に比較したり、時間的な推移を可視化できるわけです

というわけで、今回は国土交通省のGIS（地理情報システム）ホームページから行政区域データを取得し、都道府県・市区町村の形を地図上にプロットしていきます

https://nlftp.mlit.go.jp/index.html

例によって可視化には Elixir Livebook を使います

手軽に綺麗な可視化ができるからです

https://livebook.dev/

実装したノートブックはこちら

https://github.com/RyoWakabayashi/elixir-learning/blob/main/livebooks/gis/gis_polygon.livemd

## 出典

「国土数値情報（行政区域データ）」（国土交通省）（<https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_1.html>）を加工して作成

## 参考サイト

宙畑の Python 実装を参考にしています

https://sorabatake.jp/28126/

## 「衛星データ利用環境整備・ソリューション開発支援事業」

この記事は経済産業省の「衛星データ利用環境整備・ソリューション開発支援事業」に関連して、衛星データを Elixir で処理する方法を研究するために書いています

https://sdu.go.jp/

近日、本記事で扱う地理情報と、 Tellus から取得する衛星データの紐付けについて記事を投稿します

https://www.tellusxdp.com/ja/

## 実行環境

- macOS Ventura 13.1
- Rancher Desktop 1.7.0

Livebook 0.8.0 の Docker　イメージを元にしたコンテナで動かしました

https://hub.docker.com/layers/livebook/livebook/0.8.0/images/sha256-6ef3fdd439fd72af9ae01d58e947664a46a8493f031e179cc0e91c08ddac02e5?context=explore

コンテナ定義はこちらを参照

https://github.com/RyoWakabayashi/elixir-learning

## セットアップ

Livebook のセットアップセル（先頭のセル）に以下を入力して実行します

```elixir
Mix.install(
  [
    {:nx, "~> 0.4"},
    {:evision, "~> 0.1"},
    {:exla, "~> 0.4"},
    {:geo, "~> 3.4"},
    {:jason, "~> 1.4"},
    {:kino, "~> 0.8"},
    {:kino_maplibre, "~> 0.1.3"}
  ],
  config: [
    nx: [
      default_backend: EXLA.Backend
    ]
  ]
)
```

以下の依存モジュールをインストールしています

- [Nx](https://github.com/elixir-nx/nx/tree/main/nx): 行列演算
- [Evision](https://github.com/cocoa-xu/evision): 画像処理
- [EXLA](https://github.com/elixir-nx/nx/tree/main/exla): 行列演算の高速化
- [Geo](https://github.com/bryanjos/geo): 地理情報システム
- [Jason](https://github.com/michalmuskala/jason): JSONパーサー
- [Kino](https://github.com/livebook-dev/kino): Livebook での入出力
- [Kino MapLibre](https://github.com/livebook-dev/kino_maplibre): 地図出力

また、 Nx による行列演算を EXLA を使って高速化する設定をしています

## 行政区域データのダウンロード

[行政区域データのダウンロードページ](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_1.html)から取得したい都道府県、対象年を選択します

大正から令和まで、かなり広い期間のデータがダウンロードできるようになっています

![スクリーンショット 2023-01-02 10.31.17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/aaf41d88-e99d-7af5-590c-c0bd07034883.png)

初回はアンケートへの回答などを求められるので、それに回答してからダウンロードします

ダウンロードした ZIP ファイルを展開します

今回は大分県令和4年のデータを `/tmp/GML` に展開したものとして進めます

展開したファイルの中身を確認して見ましょう

```elixir
gml_dir = "/tmp/GML"
```

```elixir
File.ls!(gml_dir)
```

![スクリーンショット 2023-01-02 10.38.31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/19963899-d7b6-c4cf-5f6a-f28bd2662843.png)

`.geojson` というファイルがあります

GeoJSON は地理情報を JSON 形式で表現するための形式です

https://geojson.org/

Elixir では Geo モジュールを使うことで扱いやすくすることができます

## GeoJSON の読込

geojson ファイルを開いて見ましょう

まず `.geojson` ファイルのファイル名を取得します

```elixir
geojson_file =
  gml_dir
  # ファイル一覧取得
  |> File.ls!()
  # `.geojson` で終わるもののうち先頭を取得
  |> Enum.find(& String.ends_with?(&1, ".geojson"))
```

![スクリーンショット 2023-01-02 21.51.28.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/b0fd3f56-2b3a-b893-0383-2d740c0fea0f.png)

ファイルを Jason でパースした後、更に Geo.JSON でパースします

```elixir
geojson_data =
  [gml_dir, geojson_file]
  |> Path.join()
  |> File.read!()
  |> Jason.decode!()
  |> Geo.JSON.decode!()
```

![スクリーンショット 2023-01-02 21.52.47.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/2759998f-2515-a9e3-0b9d-6f4e4b499d19.png)

パースした結果は以下のような形式になっています

```elixir
%Geo.GeometryCollection{
  geometries: [
    # 県全体の地理情報
    %Geo.MultiPolygon{
      coordinates: [
        # 一つのポリゴン=多角形
        [
          [
            {<経度>, <緯度>},
            {<経度>, <緯度>},
            {<経度>, <緯度>},
            ...
          ]
        ],
        # 一つのポリゴン
        [
          [
            {<経度>, <緯度>},
            {<経度>, <緯度>},
            {<経度>, <緯度>},
            ...
          ]
        ],
        ...
      ]
      srid: nil,
      properties: %{
        "N03_001" => "<県名>",
        "N03_002" => nil,
        "N03_003" => nil,
        "N03_004" => nil,
        "N03_007" => nil
      }
    },
    # 市の地理情報
    %Geo.Polygon{
      coordinates: [
        # 一つのポリゴン
        [
          [
            {<経度>, <緯度>},
            {<経度>, <緯度>},
            {<経度>, <緯度>},
            ...
          ]
        ]
      ],
      srid: nil,
      properties: %{
        "N03_001" => "<県名>",
        "N03_002" => nil,
        "N03_003" => nil,
        "N03_004" => "<市名>",
        "N03_007" => "<行政コード>"
      }
    },
    %Geo.Polygon{
      ...
```

地理情報は、このように緯度と経度で表される地球上の点を複数列挙し、その点を結ぶ多角形（ポリゴン）で表現されます

## Smart Cell による可視化

Livebook ではこの GeoJSON データを Smart Cell に渡すだけで簡単に可視化できます

セルを追加するときに `+Smart` を開き、 `Map` を選択します

![スクリーンショット 2023-01-02 22.04.20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/c45b4f1e-9b01-5313-7dd9-f432e7a15e75.png)

追加された　Smart Cell　に値を指定します

- MAP STYLE: default
- CENTER: 137.5, 36 (日本のおおよその中心座標)
- ZOOM: 4
- Source: geojson_data
- Type: line

![スクリーンショット 2023-01-02 22.10.36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/19a6c46d-edc0-096d-7756-b92ee2a7dca8.png)

セルを実行すると、地図上に都道府県・市区町村の形が表示されます

この地図はズームしたり表示範囲を動かしたり、対話的（インタラクティブ）に動作します

![oita.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/b43ab903-4a70-66f0-4bc5-eb019680f7a6.gif)

この黒い線の折れている点全ての座標が GeoJSON の中に書かれているわけです

## Smart Cell を使わない可視化

Smart Cell を使わず、最初から当該都道府県を中心にズームして、都道府県内を塗りつぶしてみましょう

都道府県だけを使うので、先頭のデータだけを取り出します

```elixir
prefecture_geojson =
  %Geo.GeometryCollection{
    geometries: [Enum.at(geojson_data.geometries, 0)]
  }
```

中心座標を取得します

```elixir
longitudes =
  prefecture_geojson.geometries
  |> Enum.map(& &1.coordinates)
  |> List.flatten()
  |> Enum.map(& elem(&1, 0))

latitudes =
  prefecture_geojson.geometries
  |> Enum.map(& &1.coordinates)
  |> List.flatten()
  |> Enum.map(& elem(&1, 1))

center =
  {
    (Enum.min(longitudes) + Enum.max(longitudes)) / 2,
    (Enum.min(latitudes) + Enum.max(latitudes)) / 2
  }
```

MapLibre に GeoJSON を渡してプロパティを指定し、可視化します

```elixir
MapLibre.new(center: center, zoom: 7.5)
|> MapLibre.add_geo_source("prefecture_geojson", prefecture_geojson)
|> MapLibre.add_layer(
  id: "prefecture",
  source: "prefecture_geojson",
  type: :fill,
  paint: [fill_color: "#000000"]
)
```

![スクリーンショット 2023-01-02 22.25.02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/97b6a9ef-2da3-b7c3-2a92-0bc24f5bb3d8.png)

## 市区町村の可視化

特定の市区町村だけを取り出す場合、行政コードを指定します

行政コードは以下のサイトを参照してください

https://nlftp.mlit.go.jp/ksj/gml/codelist/AdminAreaCd.html

例えば大分県大分市の行政コードは 44201 です

GeoJSON から大分市のデータだけを抽出します

```elixir
city_geojson =
  %Geo.GeometryCollection{
    geometries: Enum.filter(geojson_data.geometries, & &1.properties["N03_007"] == "44201")
  }
```

中心座標を求める関数を定義しておきます

```elixir
get_center = fn geometries ->
  coordinates =
    geometries
    |> Enum.map(& &1.coordinates)
    |> List.flatten()

  longitudes = Enum.map(coordinates, & elem(&1, 0))
  latitudes = Enum.map(coordinates, & elem(&1, 1))

  {
    (Enum.min(longitudes) + Enum.max(longitudes)) / 2,
    (Enum.min(latitudes) + Enum.max(latitudes)) / 2
  }
end
```

中心座標を取得します

```elixir
center = get_center.(city_geojson.geometries)
```

地図に表示します

```elixir
MapLibre.new(center: center, zoom: 9)
|> MapLibre.add_geo_source("city_geojson", city_geojson)
|> MapLibre.add_layer(
  id: "city",
  source: "city_geojson",
  type: :fill,
  paint: [fill_color: "#000000"]
)
```

![スクリーンショット 2023-01-02 23.08.28.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/7313c5b2-0c73-430d-2c48-3cb692aae451.png)

離れ小島までかなり細かく描画されています

## Evision による可視化

では Evision を使って、この多角形を画像データにしてみましょう

まず、対象行政区域の端の緯度経度を取得します

```elixir
# 緯度経度の最大最小を求める
coordinates =
  city_geojson.geometries
  |> Enum.map(& &1.coordinates)
  |> List.flatten()

longitudes = Enum.map(coordinates, & elem(&1, 0))
latitudes = Enum.map(coordinates, & elem(&1, 1))

min_longitude = Enum.min(longitudes)
max_longitude = Enum.max(longitudes)
min_latitude = Enum.min(latitudes)
max_latitude = Enum.max(latitudes)
```

画像にするときのサイズを指定します

```elixir
{height, width} = {1280, 1280}
```

緯度経度の最小から最大をピクセル数の0から幅・高さに正規化します

大分県の場合は以下のようになります

- 横方向:
  - 正規化前: 131.41873417604154 から 131.96288810590158
  - 正規化後: 0 から 1280
- 縦方向:
  - 正規化前: 33.29024244133575 から 33.069713315217484
  - 正規化後: 0 から 1280

:::note warn
縦方向は北緯なので大小が逆転します
:::

```elixir
normalized_points =
  city_geojson.geometries
  |> Enum.map(& &1.coordinates)
  |> Enum.map(fn coordinate ->
    coordinate
    |> Enum.at(0)
    |> Enum.map(fn {x, y} ->
      [
        trunc((x - min_longitude) * width / (max_longitude - min_longitude)),
        # 縦方向は緯度が北緯の場合逆転するため、高さから引く
        trunc(height - (y - min_latitude) * height / (max_latitude - min_latitude))
      ]
    end)
    |> Nx.tensor(type: :s32)
  end)
```

多角形を描画するための空画像を用意します

透明度を使うため、 BGRA の4つに 0 を指定します

```elixir
empty_mat =
  # 全て黒の不透明
  [0, 0, 0, 255]
  |> Nx.tensor(type: :u8)
  |> Nx.tile([height, width, 1])
  |> Evision.Mat.from_nx_2d()
```

![スクリーンショット 2023-01-02 23.24.01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/29879461-93c9-bc93-82da-9da4c40b454a.png)

`Evision.fillPoly` に元画像、多角形の配列、色を渡すことで画像に多角形を描画できます

```elixir
img = Evision.fillPoly(empty_mat, normalized_points, {0, 0, 0, 0})
```

![スクリーンショット 2023-01-02 23.09.13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/a218a322-d81f-c0d3-51dd-c6075d149c12.png)

縦横 1280 ピクセル内に大分市の形を多角形として表現できました

白いところは透明になっています

これを地図上にプロットしてみます

画像を地図にプロットするため、 BASE64 にします

```elixir
img_base64 =
  Evision.imencode(".png", img)
  |> Base.encode64()
  |> then(&"data:image/png;base64,#{&1}")
```

今度は航空写真上にプロットしてみましょう

```elixir
MapLibre.new(center: center, zoom: 10, style: :terrain)
# 画像をレイヤーとして地図に重ねる
|> MapLibre.add_source(
  "city_mask",
  type: :image,
  url: img_base64,
  coordinates: [
    [min_longitude, max_latitude],
    [max_longitude, max_latitude],
    [max_longitude, min_latitude],
    [min_longitude, min_latitude]
  ]
)
|> MapLibre.add_layer(
  id: "overlay",
  source: "city_mask",
  type: :raster,
  layout: %{
    "visibility" => "visible"
  }
)
```

![スクリーンショット 2023-01-02 23.30.21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/8d8dc5d9-7129-8785-ef1d-4dbdcfd431ee.png)

道路や海岸線を見ると、かなり正確なことが分かります

## まとめ

これで都道府県や市区町村の形を画像データ化することができました

例えば衛星データから得られる植生（どれくらい植物が生えているか）をこの画像でマスクすれば、大分市内の植生の推移、などが数値化できることになります

別のオープンデータを使えば田んぼや畑も地図上にプロットできます

https://qiita.com/RyoWakabayashi/items/4d5ff1a2ae8644d3a7b2