昨日(22/1/23)から、法務局の地図データがG空間情報センターで無償ダウンロードできるようになりました。
早速アクセスしてみましたが、昨日はずっとビジーで今日ようやくデータ(zip)を2つダウンロード(25MBと50MB)。
ファイル名は市区町村コードで始まっているので、整理はしやすいですね。
中には
(1)索引ファイル(csv)
(2)地図データ(xmlをzipにそれぞれ圧縮)
が入っています。索引ファイルで見たい地番に対応する地図データを探し、解凍してxmlを開く感じでしょうか。それぞれ圧縮しているのは、圧縮率が高く、解凍すると30~40倍の容量を食うためでしょう。試しに10個解凍したところ、35倍でしたから、まとめて1回圧縮しただけだと、単純計算で25MB・50MB→875MB・1750MBになります。ただ、数が多い(25MBの方は46個、50MBの方は585個)ので、手作業で扱うのは大変そうです。
地図データは、
- 地図名、市区町村コード、市区町村名、座標系など。
- 点データ:XY座標
IDは、「P+数字9桁」
- 線データ:土地と土地の境界線データか。起点と終点が 2 の点データで設定されている。
IDは、「C+数字9桁」
- 面データ:1筆の土地のデータ。土地を囲む線が 3 の線データで指定されている。
IDは、「F+数字9桁」
- 属性データ:2~4の各点・線・面の属性(三角点・筆界点、行政区界線・筆界線、地番など)が設定されている。
面の属性データのみ「筆」として「H+数字9桁」のIDが割り振られている。
理由は分からないが、飛び地も表現できるようにするためかも。
- 図郭データ:元の地図番号や地図の種類、各地図に含まれる土地のリストなど。
各地図に含まれる土地は、「筆」データが指定されている。
という感じになっています。地番から目的土地を探すなら、
- 索引ファイルで必要なxmlファイルを探して開く。
- 属性データでほしい地番を探す。見つかったら、面データのIDを拾う。
- 面データのIDでほしい面データを探す。見つかったら、線データのIDを拾う。
- 線データのIDでほしい線データを探す。見つかったら、点データのIDを拾う。
- 点データのIDでほしい点データを探す。見つかったら、XY座標を拾い、直線で結ぶ。
- 4・5を 3 で拾ったすべての線データについて行う。
となりましょうか。検索しやすいように独自の形式を変換したほうが早いかもしれません。
さて、このデータを活用するのに必須と思われる「点データのXY座標」⇒「経緯度」変換プログラムを作りました。
ただし、地図データの最初の方にある「座標系」データが「公共座標●系」となっているものにしか使えません。「任意座標系」はその地図の中でしか使えない独自の座標系なので、経緯度に換算できません。
また、平面直角座標系は、原点から遠くなるほど誤差が大きくなるため、日本国内に19の原点が設けられ、1~19までの「系番号」を割り振っています。どの原点を参照するかは、地域ごとに決まっています(詳しくは国土地理院ウェブサイトの解説参照)。「公共座標●系」の「●」の部分が系番号です。
X座標・Y座標・系番号から、北緯・東経を返す関数です。系番号は K に投入します。度(1度未満は小数)で返しますが、ReturnType を 1 にすると、秒単位(度は1万、分は千と百の位に表示)で返します。
計算方法は、国土地理院ウェブサイトの計算式を参照してください。
<エクセルVBA>
Public Function SurveyCalc_xy2bl(ByVal x#, ByVal y#, ByVal k#, Optional ByVal ReturnType% = 0) As Double()
Dim ret#(1)
If k < 0 Or k > 19 Then GoTo Err '存在しない座標系の指定
'定数の設定
Const small_a& = 6378137, F# = 298.257222101 '楕円体の長半径(メートル)及び逆扁平率(測量法施行令3条)
Const m0# = 0.9999 '平面直角座標系のX軸上における縮尺係数(定義)
Dim r#: r = 180 / WorksheetFunction.Pi 'ρ''。ここではすべて度単位とする。
'平面直角座標系原点:https://www.gsi.go.jp/LAW/heimencho.html
Dim Lats As Variant, Lons As Variant, bLat#, bLon#
Lats = Array(0, 33, 33, 36, 33, 36, 36, 36, 36, 36, 40, 44, 44, 44, 26, 26, 26, 26, 20, 26)
Lons = Array(0, 129.5, 131, 132 + 1 / 6, 133.5, 134 + 2 / 6, 136, 137 + 1 / 6, 138.5, 139 + 5 / 6, 140 + 5 / 6, 140.25, 142.25, 144.25, 142, 127.5, 124, 131, 136, 154)
bLat = Lats(k) / r: bLon = Lons(k) / r 'bLat(baseLatitude=緯度)はφ0、bLon(baseLongtitude=経度)はλ0。三角関数はラジアンで処理するため、先に変換している。
'最初にnの計算をする。n2~n6は累乗表記の簡便のため
Dim n#, n2#, n3#, n4#, n5#, n6#, a#(0 To 5), b#(1 To 5), d#(1 To 6)
n = 1 / (2 * F - 1): n2 = n ^ 2: n3 = n ^ 3: n4 = n ^ 4: n5 = n ^ 5: n6# = n ^ 6
'A0~5、β1~5、δ1~5の計算
a(0) = 1 + n2 / 4 + n4 / 64: a(1) = -3 / 2 * (n - n3 / 8 - n5 / 64): a(2) = 15 / 16 * (n2 - n4 / 4): a(3) = -35 / 48 * (n3 - n5 * 5 / 16): a(4) = 315 / 512 * n4: a(5) = -693 / 1280 * n5
b(1) = n / 2 - n2 * 2 / 3 + n3 * 37 / 96 - n4 / 360 - n5 * 81 / 512: b(2) = n2 / 48 + n3 / 15 - n4 * 437 / 1440 + n5 * 46 / 105: b(3) = n3 * 17 / 480 - n4 * 37 / 840 - n5 * 209 / 4480: b(4) = n4 * 4397 / 161280 - n5 * 11 / 504: b(5) = n5 * 4583 / 161280
d(1) = n * 2 - n2 * 2 / 3 - n3 * 2 + n4 * 116 / 45 + n5 * 26 / 45 - n6 * 2854 / 675: d(2) = n2 * 7 / 3 - n3 * 8 / 5 - n4 * 227 / 45 + n5 * 2704 / 315 + n6 * 2323 / 945
d(3) = n3 * 56 / 15 - n4 * 136 / 35 - n5 * 1262 / 105 + n6 * 73814 / 2835: d(4) = n4 * 4279 / 630 - n5 * 332 / 35 - n6 * 399572 / 14175: d(5) = n5 * 4174 / 315 - n6 * 144838 / 6237: d(6) = n6 * 601676 / 22275
'計算の展開
With WorksheetFunction
Dim bar_a#, bar_S#, eps#, eta#, eps2#, eta2#, khi#, tLat#, tLon#
Dim s#, j% 'シグマ部分の計算に使う
bar_a = (m0 * small_a) / (1 + n) * a(0) 'Aバー
s = 0: For j = 1 To 5: s = s + a(j) * Math.Sin(2 * j * bLat): Next 'Sバーφ0のシグマ部分
bar_S = (m0 * small_a) / (1 + n) * (a(0) * bLat + s) 'Sバーφ0
eps = (x + bar_S) / bar_a: eta = y / bar_a 'εとη
s = 0: For j = 1 To 5: s = s + b(j) * Math.Sin(2 * j * eps) * .Cosh(2 * j * eta): Next
eps2 = eps - s 'ε'
s = 0: For j = 1 To 5: s = s + b(j) * Math.Cos(2 * j * eps) * .Sinh(2 * j * eta): Next
eta2 = eta - s 'η'
khi = .Asin(Math.Sin(eps2) / .Cosh(eta2)) 'χ
s = 0: For j = 1 To 6: s = s + d(j) * Math.Sin(2 * j * khi): Next
tLat = (khi + s) * r 'φ
tLon = (bLon + Math.Atn(.Sinh(eta2) / Math.Cos(eps2))) * r 'λ
End With
If ReturnType = 1 Then '度分秒形式への変換
ret(0) = Int(tLat) * 10000 + Int(tLat * 60 - Int(tLat) * 60) * 100 + (tLat * 3600 - Int(tLat * 60) * 60)
ret(1) = Int(tLon) * 10000 + Int(tLon * 60 - Int(tLon) * 60) * 100 + (tLon * 3600 - Int(tLon * 60) * 60)
Else
ret(0) = tLat
ret(1) = tLon
End If
Err:
SurveyCalc_xy2bl = ret
End Function
<VB2019>
Public Function SurveyCalc_xy2bl(x#, y#, k#, Optional ReturnType% = 0) As Double()
If k < 0 Or k > 19 Then GoTo Err '存在しない座標系の指定
'定数の設定
Dim small_a% = 6378137, F# = 298.257222101 '楕円体の長半径(メートル)及び逆扁平率(測量法施行令3条)
Dim m0# = 0.9999 '平面直角座標系のX軸上における縮尺係数(定義)
Dim r# = 180 / Math.PI 'ρ''。ここではすべて度単位とする。
'平面直角座標系原点:https://www.gsi.go.jp/LAW/heimencho.html
Dim Lats#() = {0, 33, 33, 36, 33, 36, 36, 36, 36, 36, 40, 44, 44, 44, 26, 26, 26, 26, 20, 26}
Dim Lons#() = {0, 129.5, 131, 132 + 1 / 6, 133.5, 134 + 2 / 6, 136, 137 + 1 / 6, 138.5, 139 + 5 / 6, 140 + 5 / 6, 140.25, 142.25, 144.25, 142, 127.5, 124, 131, 136, 154}
Dim bLat# = Lats(k) / r, bLon# = Lons(k) / r 'bLat(baseLatitude=緯度)はφ0、bLon(baseLongtitude=経度)はλ0。三角関数はラジアンで処理するため、先に変換している。
'最初にnの計算をする。n2~n6は累乗表記の簡便のため
Dim n# = 1 / (2 * F - 1), n2# = n ^ 2, n3# = n ^ 3, n4# = n ^ 4, n5# = n ^ 5, n6# = n ^ 6
'A0~5、β1~5、δ1~5の計算=ここから原則下から上にたどる。
Dim a#() = {1 + n2 / 4 + n4 / 64, -3 / 2 * (n - n3 / 8 - n5 / 64), 15 / 16 * (n2 - n4 / 4), -35 / 48 * (n3 - n5 * 5 / 16), 315 / 512 * n4, -693 / 1280 * n5}
Dim b#() = {0, n / 2 - n2 * 2 / 3 + n3 * 37 / 96 - n4 / 360 - n5 * 81 / 512, n2 / 48 + n3 / 15 - n4 * 437 / 1440 + n5 * 46 / 105, n3 * 17 / 480 - n4 * 37 / 840 - n5 * 209 / 4480, n4 * 4397 / 161280 - n5 * 11 / 504, n5 * 4583 / 161280}
Dim d#() = {0, n * 2 - n2 * 2 / 3 - n3 * 2 + n4 * 116 / 45 + n5 * 26 / 45 - n6 * 2854 / 675, n2 * 7 / 3 - n3 * 8 / 5 - n4 * 227 / 45 + n5 * 2704 / 315 + n6 * 2323 / 945, n3 * 56 / 15 - n4 * 136 / 35 - n5 * 1262 / 105 + n6 * 73814 / 2835, n4 * 4279 / 630 - n5 * 332 / 35 - n6 * 399572 / 14175, n5 * 4174 / 315 - n6 * 144838 / 6237, n6 * 601676 / 22275}
'計算の展開
Dim bar_a# = (m0 * small_a) / (1 + n) * a(0) 'Aバー
Dim s# = 0, j% 'シグマ部分の計算に使う
For j = 1 To 5 : s += a(j) * Math.Sin(2 * j * bLat) : Next 'Sバーφ0のシグマ部分
Dim bar_S# = (m0 * small_a) / (1 + n) * (a(0) * bLat + s) 'Sバーφ0
Dim eps# = (x + bar_S) / bar_a, eta# = y / bar_a 'εとη
s = 0 : For j = 1 To 5 : s += b(j) * Math.Sin(2 * j * eps) * Math.Cosh(2 * j * eta) : Next
Dim eps2# = eps - s 'ε'
s = 0 : For j = 1 To 5 : s += b(j) * Math.Cos(2 * j * eps) * Math.Sinh(2 * j * eta) : Next
Dim eta2# = eta - s 'η'
Dim khi# = Math.Asin(Math.Sin(eps2) / Math.Cosh(eta2)) 'χ
s = 0 : For j = 1 To 6 : s += d(j) * Math.Sin(2 * j * khi) : Next
Dim tLat# = (khi + s) * r 'φ
Dim tLon# = (bLon + Math.Atan(Math.Sinh(eta2) / Math.Cos(eps2))) * r 'λ
If ReturnType = 1 Then
Dim Lat# = Int(tLat) * 10000 + Int(tLat * 60 - Int(tLat) * 60) * 100 + (tLat * 3600 - Int(tLat * 60) * 60)
Dim Lon# = Int(tLon) * 10000 + Int(tLon * 60 - Int(tLon) * 60) * 100 + (tLon * 3600 - Int(tLon * 60) * 60)
SurveyCalc_xy2bl = {Lat, Lon}
Else
SurveyCalc_xy2bl = {tLat, tLon}
End If
Exit Function
Err:
SurveyCalc_xy2bl = {0, 0}
End Function
例えば、次のように設定すると、自動的にChromeが立ち上がり、グーグルマップで該当する経緯度の地図を表示してくれます。
下のサンプルで表示されるのは当職の事務所(神奈川総合法律事務所)です。
<VB2019>
Private Sub Test()
Dim i#() = SurveyCalc_xy2bl(-61160, -17675, 9)
Process.Start(fileName:="""C:\Program Files\Google\Chrome\Application\chrome.exe""", "https://www.google.co.jp/maps/search/" & i(0).ToString & "," & i(1).ToString)
<エクセルVBA>
Private Sub Test()
Dim i#()
i = SurveyCalc_xy2bl(-61160, -17675, 9)
Shell """C:\Program Files\Google\Chrome\Application\chrome.exe"" https://www.google.co.jp/maps/search/" & i(0) & "," & i(1)
End Sub
ちなみに、国土地理院のウェブサイトは地図や測量について丁寧に解説されており、経緯度など地図を扱うための数値の換算を計算サイトで簡単に行えるうえ、換算方法の解説までついています。
とても勉強になりますし、高校程度の知識は必要かもしれませんが、読めば(取り扱える程度には)理解できます。