Elixir Livebook で農林水産省の筆ポリゴンデータを可視化する（田畑の形を地図上にプロットする）

## はじめに

皆さん、お米食べてますか？

美味しいお米を生産してくれる田んぼは日本の宝です

というわけで、今回は農林水産省の筆（ふで）ポリゴンデータを使って、田んぼや畑を地図上にプロットしていきます

https://www.maff.go.jp/j/tokei/porigon/

:::note info
ちなみに、「筆」は登記するときの土地の単位です
:::

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

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

https://livebook.dev/

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

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

## 出典

「筆ポリゴンデータ（2022年度公開）」（農林水産省）（<https://www.maff.go.jp/j/tokei/porigon/>）を加工して作成

## 参考サイト

宙畑の 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://open.fude.maff.go.jp/)から画面下部の「筆ポリゴンデータダウンロードページ」のボタンをクリックします

![スクリーンショット 2023-01-03 0.43.52.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/36910e6d-5f02-6c22-2b93-59af0591b4f1.png)

表示されるアンケートに回答します

遷移先のダウンロードページで公開年度、都道府県、市区町村を選択し、「ダウンロード」をクリックします

![スクリーンショット 2023-01-03 0.47.02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/3a023513-3993-7fea-dc12-137ff39636aa.png)

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

今回は大分県大分市 2022 年度のデータを `/tmp/2022_442011.json` に展開したものとして進めます

拡張子は単に `.json` ですが、中身は　GeoJSON になっています

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

https://geojson.org/

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

## GeoJSON の読込

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

```elixir
geojson_data =
  json_file
  |> File.read!()
  |> Jason.decode!()
  |> Geo.JSON.decode!()
```

![スクリーンショット 2023-01-03 0.54.15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/8407dd65-62c3-2218-000a-d239d6f22c91.png)

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

:::note info
データの仕様は[こちら](https://opendata.fude.maff.go.jp/%E7%AD%86%E3%83%9D%E3%83%AA%E3%82%B4%E3%83%B3%E3%83%87%E3%83%BC%E3%82%BF%E3%81%AE%E4%BB%95%E6%A7%98%EF%BC%88OD%EF%BC%89.pdf)
:::

```elixir
%Geo.GeometryCollection{
  geometries: [
    # 田畑の地理情報
    %Geo.Polygon{
      coordinates: [
        # 一つのポリゴン
        [
          [
            {<経度>, <緯度>},
            {<経度>, <緯度>},
            {<経度>, <緯度>},
            ...
          ]
        ]
      ],
      srid: nil,
      properties: %{
        "edit_year" => <調製年度>,
        "history" => <履歴>,
        "issue_year" => <公開年度>,
        "land_type" => <耕地の種類:地目コード（100：田、200：畑）>,
        "last_polygon_uuid" => <前年筆ポリゴンID>,
        "local_government_cd" => <地方公共団体コード>,
        "point_lat" => <重心点座標緯度>,
        "point_lng" => <重心点座標経度>,
        "polygon_uuid" => <筆ポリゴンID>,
        "prev_last_polygon_uuid" => <前前年筆ポリゴンID>
      }
    },
    %Geo.Polygon{
      ...
```

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

```elixir
Enum.count(geojson_data.geometries)
```

件数をカウントしてみます

64,400 件ありました

田んぼだけだと 37,817 件

```elixir
geojson_data.geometries
|> Enum.filter(& &1.properties["land_type"] == 100)
|> Enum.count()
```

畑だけだと 26,583 件

```elixir
# 畑の件数
geojson_data.geometries
|> Enum.filter(& &1.properties["land_type"] == 200)
|> Enum.count()
```

件数が大きすぎるため、緯度経度で絞り込みます

```elixir
target_fields =
  geojson_data.geometries
  |> Enum.filter(fn field ->
    field.properties["point_lng"] >= 131.42 &&
    field.properties["point_lng"] <= 131.46 &&
    field.properties["point_lat"] >= 33.13 &&
    field.properties["point_lat"] <= 33.15 &&
    field.properties["land_type"] == 100
  end)
  |> then(fn geometries ->
    %Geo.GeometryCollection{geometries: geometries}
  end)

Enum.count(target_fields.geometries)
```

## 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: terrain(non-commercial)
- CENTER: 131.443, 33.131
- ZOOM: 16
- Source: target_fields
- Type: line

セルを実行すると、航空写真の地図上に田んぼの形が表示されます

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

![farm.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/65bff3be-9431-73eb-a0da-6555a132a0f7.gif)

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

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

Smart Cell を使わず、最初から当該区域を中心にズームして、田んぼを塗りつぶしてみましょう

中心座標を取得します

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

latitudes =
  target_fields.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: 14.5, style: :terrain)
|> MapLibre.add_geo_source("geojson", target_fields)
|> MapLibre.add_layer(
  id: "fill",
  source: "geojson",
  type: :fill,
  # 半透明の緑にする
  paint: [fill_color: "#00ff00", fill_opacity: 0.5]
)
```

![スクリーンショット 2023-01-03 1.42.01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/bb6c6b68-ebcd-3589-90cd-a26aec80ed2b.png)

## Evision による可視化

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

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

```elixir
# 緯度経度の最大最小を求める
coordinates =
  target_fields.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)

{min_longitude, max_longitude, min_latitude, max_latitude}
```

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

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

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

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

- 横方向:
  - 正規化前: 131.435353795 から 131.460065396
  - 正規化後: 0 から 1280
- 縦方向:
  - 正規化前: 33.135901549 から 33.129839302
  - 正規化後: 0 から 2560

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

```elixir
normalized_points =
  target_fields.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-03 1.46.21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/965e0f24-48dc-4fc8-9b80-18637eaaac99.png)

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

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

![スクリーンショット 2023-01-03 1.47.17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/ce4d27f5-42f2-cf68-73b8-54ccc939fa89.png)

縦 1280 ピクセル、横 2560 ピクセルに田んぼの形を多角形として表現できました

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

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

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

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

```elixir
MapLibre.new(center: center, zoom: 14.5, style: :terrain)
# 画像をレイヤーとして地図に重ねる
|> MapLibre.add_source(
  "field_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: "field_mask",
  type: :raster,
  layout: %{
    "visibility" => "visible"
  }
)
```

![スクリーンショット 2023-01-03 1.47.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/c8430437-b763-fd8b-4e5e-0876325116d2.png)

## まとめ

画像データ化できた筆ポリゴンデータに衛星データから計算した植生（どれくらい植物が生えているか）を重ねることで、田畑の栽培状況を見ることができます

別のオープンデータを使えば都道府県や市区町村も地図上にプロットできます

https://qiita.com/RyoWakabayashi/items/1a12439d2971444ba33b