都道府県県庁間の距離行列（単位：m）を作成してみました。
dist_Pref <- read.csv("https://raw.githubusercontent.com/yoshida-gisc/stockroom/master/distmat_Pref.csv")
dist_Pref[1:6,1:6]
# PrefName Hokkaido Aomori Iwate Miyagi Akita Yamagata
# 1 Hokkaido 0.0 253809.6 373582.28 534013.96 385851.16 542058.09
# 2 Aomori 253809.6 0.0 129308.02 283958.87 134229.07 288698.14
# 3 Iwate 373582.3 129308.0 0.00 161119.79 90055.14 176229.35
# 4 Miyagi 534014.0 283958.9 161119.79 0.00 174197.75 44628.38
# 5 Akita 385851.2 134229.1 90055.14 174197.75 0.00 165634.86
# 6 Yamagata 542058.1 288698.1 176229.35 44628.38 165634.86 0.00
都道府県県庁間の距離は、国土地理院 がGRS80楕円体に基づく最短距離（測地線長）を公開しています。
空間計量経済モデルの空間重み行列を作成するときなどにも使ったことがあったのですが、単位がkmで小数第一位までになっていて、m だったらなと思ったことがありました（なぜだっけ...）。
次のような疑問や希望がありえそうということで下記コードより作成しました。
- GISソフトを用いれば2点間距離は簡単に計算できるが投影座標系の区域境界をまたぐような長距離はどう計算したらよいか。
- 全国の市区町村（役所）間距離や町丁目間距離を計算しようとするとなると単位は m がよい。
library(geosphere)
library(data.table)
# PrefCapital.csv は https://www.benricho.org/chimei/latlng_data.html から
# 都道府県庁所在地の緯度経度をDLし、簡単に整形したもの
# このリンクのソースは https://www.gsi.go.jp/KOKUJYOHO/center.htm
dat <- fread("C:/Users/hoge/PrefCapital.csv")
head(dat)
# PrefCapital Lon Lat
# 1: 北海道 141.3469 43.06417
# 2: 青森県 140.7400 40.82444
# 3: 岩手県 141.1525 39.70361
# 4: 宮城県 140.8719 38.26889
# 5: 秋田県 140.1025 39.71861
# 6: 山形県 140.3633 38.24056
## geosphere::distGeo() 関数のデフォルトはWGS84楕円体のため、GRS80楕円体のパラメータに設定し直す
## WGS84楕円体とGRS0楕円体は扁平率 f がごく僅かに異なる（長半径 a は同じ）
formals(distGeo)$f <- 1/298.257222101
dist_Pref <- geosphere::distm(dat[,c("Lon","Lat")], fun = distGeo)
colnames(dist_Pref) <- dat$PrefCapital
rownames(dist_Pref) <- dat$PrefCapital
# m 単位の都道府県庁間の距離
# km 単位の https://www.gsi.go.jp/KOKUJYOHO/kenchokan.html とも整合
dist_Pref[1:6,1:6]
# 北海道 青森県 岩手県 宮城県 秋田県 山形県
# 北海道 0.0 253809.6 373582.28 534013.96 385851.16 542058.09
# 青森県 253809.6 0.0 129308.02 283958.87 134229.07 288698.14
# 岩手県 373582.3 129308.0 0.00 161119.79 90055.14 176229.35
# 宮城県 534014.0 283958.9 161119.79 0.00 174197.75 44628.38
# 秋田県 385851.2 134229.1 90055.14 174197.75 0.00 165634.86
# 山形県 542058.1 288698.1 176229.35 44628.38 165634.86 0.00
市区町村や町丁目のポリゴンデータを持っていればその重心点の緯度経度と上のコードから、測地線距離に基づいた市町村間、町丁目間の距離行列が生成できる。市区町村の役所の緯度経度は 国土地理院 がまとめているので、PDFを整形し、60進法表記を @uri さんのkuniezuパッケージ や pazerパッケージを活用して10進法に変換後、上のコードから市区町村役所間の距離行列が作成できそう。
次は、OSMの道路ネットワークを使った最短経路距離について、rgrass7パッケージでやってみよう...。