地図XMLの仕様書における測地系記載変更に対するメモ

[koswatana]

地図XML(2024/07/09時点)の仕様書において、「測地系は日本測地系2000 (JGD2000) 又は日本測地系2011 (JGD2011)」と追記があり、データ上はこれらの区別がつかないことが判明。

• https://github.com/user-attachments/files/16144453/000116464.pdf (2024/07/09時点、P.2)

現状、mojxml2geojsonでは、
(1) 平面直角座標系 (JGD2000) →緯度経度 (JGD2000)
(2) 緯度経度 (JGD2000) → 緯度経度 (JGD2011)

• mojxml2geojson/src/package/MojXMLtoGeoJSON.py

  Line 73 in aa854b1

  # 測地系の変換 to 6668(JGD2011)

  
  の座標変換を行なっているため、入力がJGD2011の場合、不要な (2) の変換をかけてしまうことになる。

---

方針

• （検証中）調査の結果、現状の実装はJGD2010の平面直角座標を、JGD2011の緯度経度に変換することを想定していたが、pyprojでパラメータファイル（touhokutaiheiyouoki2011.gsb）を指定していないため、適用範囲内であってもJGD2000の緯度経度が返っている
• 地図XMLの仕様上、入力がJGD2000かJGD2011かは判定できない
  • JGD2011が入力された場合、パラメータファイルを適用してしまうと二重で2011の変換処理をかけてしまうことになる
• そのため地図XMLで判別できるようになるまでは、mojxml2geojsonでは入力の測地系に関わらず平面直角座標から緯度経度への変換のみを行うものとする
  • 現状の実装には影響がない
    • 入力がどちらかわからない以上fromは変更のしようがない
    • toのEPSGが6668（JGD2011）になってしまっているので、修正するか？
• ドキュメントは修正する

---

地籍調査

参考資料1

• https://www.sangiin.go.jp/japanese/annai/chousa/rippou_chousa/backnumber/2015pdf/20151001118.pdf
• 地籍調査は1951(昭和26)年に開始
• 地籍調査の進捗率は2014(平成26)年度末で51%
  • 2011年以前の測量成果に基づく登記所備付地図が多い割合で存在する可能性。

参考資料2

• 「測地成果2000」と「測地成果2011」
  • https://www.gsi.go.jp/sokuchikijun/jgd2000-2011.html
• 平成23年（2011年）東北地方太平洋沖地震に伴う基準点測量成果の改定の詳細
  • https://www.gsi.go.jp/common/000064457.pdf
• 地殻変動により旧成果と新成果で東北地方を中心に最大で5mを超える差
  • 登記所備付地図の利用シーンによっては看過できない影響である可能性。

---

mojxml2geojson/src/package/ConvertDatum.py

Line 7 in aa854b1

def xy2bl(coordinates, from_epsg, to_epsg):

mojxml2geojson/src/package/MojXmlDef.py

Line 4 in aa854b1

def GetCrs(crsText: str):

[koswatana]

pyprojが2011のパラメータファイルを使っているかどうか
使っている場合inputが2011に対して再度2011の処理をかけて、存在しない測地系になってしまう

[koswatana]

• pyprojはEPSGを指定するだけではパッチが提供されていなさそう
  石巻市役所
  JGD2000 10系 EPSG:2452
  y, x | 40978.932, -173713.629

root@68ea51e6c429:/# python
Python 3.10.14 (main, Jul  2 2024, 22:12:36) [GCC 12.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from pyproj import Transformer

① 平面直角座標（JGD2000）→ 緯度経度（JGD2000）
>>> Transformer.from_proj(2452, 4612).transform(-173713.629, 40978.932) # JGD2000
(38.434192786760825, 141.3027177136399)

② 平面直角座標（JGD2000）→ 緯度経度（JGD2011）
>>> Transformer.from_proj(2452, 6668).transform(-173713.629, 40978.932) #JGD2011
(38.434192786760825, 141.3027177136399)

③ 平面直角座標（JGD2000）→ 緯度経度（JGD2000）→ 緯度経度（JGD2011）
>>> Transformer.from_proj(4612, 6668).transform(38.434192786760825, 141.3027177136399)
(38.434192786760825, 141.3027177136399)

---

参考

https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/xy2blf.html

https://vldb.gsi.go.jp/sokuchi/surveycalc/patchjgd/index.html

[koswatana]

>>> Transformer.from_proj(6668, 6678).transform(38.434192786760825, 141.3027177136399)
(-173713.62900000004, 40978.93200000026)
>>> Transformer.from_proj(6678,2452).transform(-173713.62900000004, 40978.93200000026)
(-173713.62900000004, 40978.93200000026)

[koswatana]

https://pyproj4.github.io/pyproj/stable/transformation_grids.html

[koswatana]

https://github.com/tohka/JapanGridShift/blob/master/README.ja.md

wget https://github.com/tohka/JapanGridShift/raw/master/gsb_files/touhokutaiheiyouoki2011.gsb

NTv2 (*.gsb) の仕様について
https://qiita.com/tohka383/items/e73c7235efc15efe2c1b#ntv2-gsb-%E3%81%AE%E4%BB%95%E6%A7%98%E3%81%AB%E3%81%A4%E3%81%84%E3%81%A6

EPSG の towgs84 問題について
https://qiita.com/tohka383/items/43434990e30ad76f6651

[koswatana]

pyprojでtouhokutaiheiyouoki2011を適用する方法

root@68ea51e6c429:/# wget https://github.com/tohka/JapanGridShift/raw/master/gsb_files/touhokutaiheiyouoki2011.gsb
root@68ea51e6c429:/# python
Python 3.10.14 (main, Jul  2 2024, 22:12:36) [GCC 12.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.

from pyproj import CRS, Transformer

jgd2000 = CRS.from_epsg(4612)  # EPSGコードでJGD2000を定義
jgd2011 = CRS.from_epsg(6668)  # EPSGコードでJGD2011を定義

grid_shift_file = "/touhokutaiheiyouoki2011.gsb"

transformer = Transformer.from_pipeline(
    f"+proj=pipeline "
    f"+step +proj=latlong +ellps=GRS80 "
    f"+step +proj=hgridshift +grids={grid_shift_file} "
)

longitude = 141.302717714
latitude = 38.434192787
jgd2011_longitude, jgd2011_latitude = transformer.transform(longitude, latitude)

jgd2011_longitude_formatted = f"{jgd2011_longitude:.9f}"
jgd2011_latitude_formatted = f"{jgd2011_latitude:.9f}"

print(f"JGD2011座標: 緯度 = {jgd2011_latitude_formatted}, 経度 = {jgd2011_longitude_formatted}")

> JGD2011座標: 緯度 = 38.434178585, 経度 = 141.302768728

• 地理院のツールと比較して一致することを確認
  • 緯度の値が合わない
    https://vldb.gsi.go.jp/sokuchi/surveycalc/patchjgd/index.html
    

参考：パラメータファイルを指定しないでEPSGのみ指定した場合

from pyproj import CRS, Transformer

# JGD2000とJGD2011のCRSを定義
jgd2000 = CRS.from_epsg(4612)  # EPSGコードでJGD2000を指定
jgd2011 = CRS.from_epsg(6668)  # EPSGコードでJGD2011を指定
transformer = Transformer.from_crs(jgd2000, jgd2011)
longitude = 141.302717714
latitude = 38.434192787

# 座標変換
jgd2011_x, jgd2011_y = transformer.transform(longitude, latitude)

print(f"JGD2011座標: X = {jgd2011_x}, Y = {jgd2011_y}")

> JGD2011座標: X = 141.302717714, Y = 38.434192787

• EPSGを指定しているが入力値と変化がないことを確認