Powered by AppSignal & Oban Pro
Would you like to see your link here? Contact us

国土交通省 GIS 行政区域データ

livebooks/gis/gis_polygon.livemd

国土交通省 GIS 行政区域データ

Mix.install(
  [
    {:nx, "~> 0.9"},
    {:evision, "~> 0.2"},
    {:exla, "~> 0.9"},
    {:geo, "~> 3.5"},
    {:jason, "~> 1.4"},
    {:kino, "~> 0.14"},
    {:kino_maplibre, "~> 0.1"}
  ],
  config: [nx: [default_backend: EXLA.Backend]]
)

行政区域データの取得

国土交通省の行政区域データをダウンロードする

出典:「国土数値情報(行政区域データ)」(国土交通省)()を加工して作成

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

GeoJSON の読込

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

Smart Cell による可視化

MapLibre.new(center: {137.5, 36.0}, zoom: 4)
|> MapLibre.add_geo_source("geojson_data", geojson_data)
|> MapLibre.add_layer(
  id: "geojson_data_line_1",
  source: "geojson_data",
  type: :line,
  paint: [line_color: "#000000", line_opacity: 1]
)

Smart Cell を使わない可視化

prefecture_geojson = %Geo.GeometryCollection{
  geometries: [Enum.at(geojson_data.geometries, 0)]
}
longitudes =
  prefecture_geojson.geometries
  |> Enum.at(0)
  |> Map.get(:coordinates)
  |> List.flatten()
  |> Enum.map(&elem(&1, 0))

latitudes =
  prefecture_geojson.geometries
  |> Enum.at(0)
  |> Map.get(:coordinates)
  |> List.flatten()
  |> Enum.map(&elem(&1, 1))

center = {
  (Enum.min(longitudes) + Enum.max(longitudes)) / 2,
  (Enum.min(latitudes) + Enum.max(latitudes)) / 2
}
MapLibre.new(center: center, zoom: 7.5)
|> MapLibre.add_geo_source("prefecture_geojson", prefecture_geojson)
|> MapLibre.add_layer(
  id: "prefecture_line",
  source: "prefecture_geojson",
  type: :fill,
  paint: [fill_color: "#000000"]
)

市区町村の可視化

行政区域コード一覧

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"]
)

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)

{min_longitude, max_longitude, min_latitude, max_latitude}
{height, width} = {1280, 1280}
# 緯度経度の最小から最大をピクセル数の0から幅・高さに正規化する
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)
# 空画像を用意する
# 全て黒の不透明
empty_mat =
  [0, 0, 0, 255]
  |> Nx.tensor(type: :u8)
  |> Nx.tile([height, width, 1])
  |> Evision.Mat.from_nx_2d()
# ポリゴンを透明色で塗りつぶす
img = Evision.fillPoly(empty_mat, normalized_points, {0, 0, 0, 0})
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"
  }
)