LoginSignup
18
14

More than 1 year has passed since last update.

@RyoWakabayashi(若林 諒)in株式会社オーイーシー

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

Last updated at Posted at 2023-01-04

はじめに

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

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

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

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

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

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

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

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

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

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

出典

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

参考サイト

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

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

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

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

実行環境

  • macOS Ventura 13.1
  • Rancher Desktop 1.7.0

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

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

セットアップ

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

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 による行列演算を EXLA を使って高速化する設定をしています

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

行政区域データのダウンロードページから取得したい都道府県、対象年を選択します

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

スクリーンショット 2023-01-02 10.31.17.png

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

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

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

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

gml_dir = "/tmp/GML"

File.ls!(gml_dir)

スクリーンショット 2023-01-02 10.38.31.png

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

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

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

GeoJSON の読込

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

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

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

スクリーンショット 2023-01-02 21.51.28.png

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

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

スクリーンショット 2023-01-02 21.52.47.png

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

%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

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

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

スクリーンショット 2023-01-02 22.10.36.png

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

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

oita.gif

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

Smart Cell を使わない可視化

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

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

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

中心座標を取得します

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 を渡してプロパティを指定し、可視化します

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

市区町村の可視化

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

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

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

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

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

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

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

中心座標を取得します

center = get_center.(city_geojson.geometries)

地図に表示します

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

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

Evision による可視化

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

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

# 緯度経度の最大最小を求める
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)

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

{height, width} = {1280, 1280}

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

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

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

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

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 を指定します

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

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

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

スクリーンショット 2023-01-02 23.09.13.png

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

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

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

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

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

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

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

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

まとめ

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

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

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

18
14
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
18
14