下の図のような、最も道のりが近い県庁毎に日本地図を塗分けます。俗に言うボロノイ図というやつでしょうか。
なお、ここでいう道のりは、OpenStreetMapのWayの内、highwayタグが含まれているものを対象としています。海路の扱いにより離島は含まれたり含まれていなかったりしますが、政治的な意図はありませんので悪しからず。

※東京都庁と埼玉県庁、愛媛県庁と高知県庁は不幸にも同じ色で塗分けられてしまいました。

また、各県庁からの道のりの可視化も行います。

グラフ（地図）の読み込みや検索は、なんだか最近グラフの検索をすることが多いので作ってみたgithub.com/takoyaki-3/goraphを使っています。

地図のダウンロード

初めに、地図をダウンロードします。
転送速度が遅く、ダウンロードにかなり時間がかかりますが、辛抱強く待ちます。

$ curl -o japan-latest.osm.pbf https://download.geofabrik.de/asia/japan-latest.osm.pbf

県庁の緯度経度情報の作成

続いて、県庁の緯度経度をCSVにまとめます。
県庁緯度経度データはみんなの知識ちょっと便利帳に掲載されているものを使用しました。

lat,lon
43.06417,141.34694
40.82444,140.74
39.70361,141.1525
38.26889,140.87194
39.71861,140.1025
38.24056,140.36333
37.75,140.46778
36.34139,140.44667
36.56583,139.88361
36.39111,139.06083
35.85694,139.64889
35.60472,140.12333
35.68944,139.69167
35.44778,139.6425
37.90222,139.02361
36.69528,137.21139
36.59444,136.62556
36.06528,136.22194
35.66389,138.56833
36.65139,138.18111
35.39111,136.72222
34.97694,138.38306
35.18028,136.90667
34.73028,136.50861
35.00444,135.86833
35.02139,135.75556
34.68639,135.52
34.69139,135.18306
34.68528,135.83278
34.22611,135.1675
35.50361,134.23833
35.47222,133.05056
34.66167,133.935
34.39639,132.45944
34.18583,131.47139
34.06583,134.55944
34.34028,134.04333
33.84167,132.76611
33.55972,133.53111
33.60639,130.41806
33.24944,130.29889
32.74472,129.87361
32.78972,130.74167
33.23806,131.6125
31.91111,131.42389
31.56028,130.55806
26.2125,127.68111

モジュールのインポート

ここからはプログラムです。全体像はGithubにもアップしています。

使用するライブラリをインポートします。

モジュール	用途
fmt	入出力
os	ファイルの読み込み書き
log	エラー出力
strconv	文字列と浮動小数点の変換
github.com/takoyaki-3/goraph	グラフモジュール
github.com/takoyaki-3/goraph/geometry/h3	最寄りノードを見つけるためのインデックス用
github.com/takoyaki-3/goraph/loader/osm	OSMの読み込み用
github.com/takoyaki-3/goraph/search	グラフ検索用
github.com/takoyaki-3/gomc	CSVファイルの読み込み
encoding/csv	実行結果CSV出力用

package main

import (
	"fmt"
	"os"
	"log"
	"strconv"
	"github.com/takoyaki-3/goraph"
	"github.com/takoyaki-3/goraph/geometry/h3"
	"github.com/takoyaki-3/goraph/loader/osm"
	"github.com/takoyaki-3/goraph/search"
	"github.com/takoyaki-3/gomc"
	"encoding/csv"
)

グラフ（地図データ）の読み込みから出力

さらに、地図データの読み取りから結果の出力まで、一気に行います。

グラフの読み込み

今回はOpenStreetMapのデータを読み込むので、github.com/takoyaki-3/goraph/loader/osmのLoad関数を用いて、グラフとして読み込みます。

func main(){

	// Graph load
	fmt.Println("Graph load")
	g := osm.Load("japan-latest.osm.pbf")

h3インデックスの作成

県庁の最寄りのノードを探すため、uberの開発しているオープンソースであるh3を用いたインデックスを作成します。
最寄りノードを検索する際、グラフに含まれるすべてのノードを検索対象とすると非常に長い時間がかかるので、h3インデックスを用いて検索対象を絞り込んでいます。

	// Make index
	fmt.Println("Make index")
	h3indexes := h3.MakeH3Index(g, 9)

ボロノイ図のbaseの読み込み

さらに、先ほど作成した県庁の緯度経度CSVファイルを読み込みます。

	// Load base point
	fmt.Println("Load base point")
	bases := []int64{}
	titles,records := gomc.ReadCSV("./base.csv")
	for _,v:=range records{
		lat,_ := strconv.ParseFloat(v[titles["lat"]], 64)
		lon,_ := strconv.ParseFloat(v[titles["lon"]], 64)
		bases = append(bases,h3.Find(g,h3indexes,goraph.LatLon{lat,lon},9))
	}

ボロノイ図の作成

github.com/takoyaki-3/goraph/searchモジュールのVoronoi関数にグラフと基地ノードリストを投げると、各ノードがどの基地（今回は県庁）に属するかを返してくれます。

	// Create voronoi diagram
	fmt.Println("Create voronoi diagram")
	nodes := search.Voronoi(g,bases)

結果の出力

返ってきた結果をCSVとして出力します。
全ノードがどの基地に属するかを出力すると、ファイルサイズが大きくなってしまうので、100ノードに1つの割合で出力しています。

	// Output
	fmt.Println("Output")
	wf, err := os.Create("./output.csv")
	if err != nil {
		log.Println(err)
	}
	defer wf.Close()

	w := csv.NewWriter(wf) // utf8
	w.Write([]string{"base","lat","lon"})
	for k,v := range nodes{
		if k%100!=0{
			continue
		}
		lat := strconv.FormatFloat(g.LatLons[k].Lat, 'f', -1, 64)
		lon := strconv.FormatFloat(g.LatLons[k].Lon, 'f', -1, 64)
		w.Write([]string{strconv.Itoa(int(v)),lat,lon})
	}
	w.Flush()
}

可視化

最後に、kepler.glにより、プログラムの出力を可視化します。

kepler.glは、大変便利なことに、列名にlat,lon,geometryなどが含まれるCSVファイルをドラッグアンドドロップするだけで可視化してくれます。

先ほど作成したプログラムにより出力されたoutput.csvをkepler.glにドラッグアンドドロップすると、ノードがすべてプロットされます。

kepler.glのレイヤー設定より、表示する色をbaseによって分けるよう設定すると、各県庁毎に塗分けられます。

完成です。

各県庁からの距離も可視化してみる

ちなみに、上記プログラムのボロノイ図の作成から先を次のように変更することで、次のような各県庁からの距離を可視化することもできます。
プログラムの全体像はGithubにもアップしています。

最寄りのベースノードからの距離をすべて算出

	// Calculate the cost to all nodes
	cost := search.AllDistance(g,bases)

結果の出力

こちらは、1kmごとに色を分け、km表記の1の位のみを出力することで、同一距離でいける範囲を分かりやすくしています。
すべてのノードを出力すると結果が膨大となるので、10ノードに1つ出力しています。

	// Output
	wf, err := os.Create("./allnode_distance.csv")
	if err != nil {
		log.Println(err)
	}
	defer wf.Close()

	w := csv.NewWriter(wf) // utf8
	w.Write([]string{"d","lat","lon"})
	for k,v := range cost{
		if k%10!=0{
			continue
		}
		lat := strconv.FormatFloat(g.LatLons[k].Lat, 'f', -1, 64)
		lon := strconv.FormatFloat(g.LatLons[k].Lon, 'f', -1, 64)
		w.Write([]string{strconv.Itoa(int(v/1000.0)%10),lat,lon})
	}
	w.Flush()
}