Future Tech Blog
フューチャー技術ブログ
Programmingカテゴリ

go-projを用いて日本測地系/世界測地系の変換

はじめに

TIG DXユニットの真野です。

ある緯度経度の座標を日本測地系から世界測地系へ変換する際に、everystreet/go-projというパッケージを用いました。私にとって初めてのcgoを用いたライブラリ利用であり、環境構築に少し悩んだので手順をまとめておきます。

なお、類似にgo-spatial/projもありますが、そちらは日本測地系がサポートされていないようでした。測地系変換の知識が無いのでノータッチですがPure Go実装にできると嬉しいですね。

日本測地系 / 世界測地系 とは

今泉さんの 郵便番号・住所・緯度経度の体系について に記載がありますので参考ください。東京付近では日本測地系と世界測地系との誤差が450m程度あるらしく、大抵のユースケースではコードを統一した方が良さそうです。

他にも以下のサイトが参考になります。

日本測地系 / 世界測地系 の変換について

日本測地系への対応 - プログラマー’sペイジ さんによると、概ね3つの変換パターンがあるようです。

  1. 1次式による近似
    • ググるとすぐ出てくる↓の式
      • lonJ = lonW + latW * 0.000046047 + lonW * 0.000083049 – 0.010041
      • latJ = latW + latW * 0.00010696 – lonW * 0.000017467 – 0.0046020
    • シンプルで素晴らしい。約2.5~2.8メートルズレる
  2. 3次元直交座標系に変換して平行移動
    • PROJ が利用している方式
    • 3に比べて、約0.2~0.4メートルズレる
  3. 地域ごとの変換パラメータで変換する
    • パラメータファイルを用いた最も精度が高い方式
    • Web版 TKY2JGD といったツールがある

今回は1は許容できないけど、2のレベルであれば問題ないユースケースだったので、PROJを用いることにしました。

※ズレの計算は、3のWeb版TKY2JGDで 日本測地系→世界測地系に変換したものを正とし、1,2の変換後の緯度経度を、測量計算(距離と方位角の計算) のWebツールで計測した。場所は五反田駅など関東の複数地点で試しました

この変換を覚えると嬉しいの？

日本測地系は、旧日本測地系（Tokyo Datum）と呼ばれるだけあり、明治時代に作られ、2002年4月1日まで利用されていた規格です。それ以降は日本測地系2000（JGD2000）あるいは、日本測地系2011（JGD2011）が使われています。JGD2000とJGD2011のどちらも世界測地系と言える経度・緯度の体系です。GD2011へ移行したきっかけは、東日本大震災による大規模な地殻変動の発生への対応だそうで、基本的にJGD2000と、JGD2011はTokyo Datumに比べると、スマホアプリのマップにピンを立ててナビゲートするようなケースでは同一と見なして良いレベルです。具体的には日本近辺で5cm程度の差があるとのことです。

何が言いたいかというと、2002年以降に公開されたデータは大概、世界測地系であるJGD2000またはJGD2011です。そのため、日本測地系→世界測地系にしたいユースケースってそんなに多くないと思います（古いデータを移行したいとか、その古いデータを元にGoogle Mapsなどに表示させたいなど？）。この記事に書かれたサンプルコードは、古い資産を活用したい場合など、限定された場面で参考になるんだと認識してもらえればです。

例えば、国土交通省の街区レベル位置参照情報はJGD2000形式です。同じく国土交通省の国有林野データはJGD2011でした。少なくても公共機関から取得できるようなデータは世界測地系になっているのかなと思います。

PROJとは

PROJ はこの界隈では有名な、地図投影や測地変換が行えるツールで、コマンドラインから呼び出しと、ライブラリのようにも使えるAPIも提供しています。

コマンドラインから日本測地系から世界測地系に変換する例です。

# 日本測地系（EPSG 4301）から世界測地系（EPSG 4612）に変換
$ echo "139.477799479 35.4891015625" | cs2cs +init=epsg:4301 +to +init=epsg:4612 -f "%.10f"
139.4745977186  35.4923543046 0.0000000000

参考: http://pen.envr.tsukuba.ac.jp/~torarimon/?%C2%AC%C3%CF%B7%CF%3A+datum

各言語からPROJを用いる

Java（proj4j）, JavaScript（proj4js）, Python（pyproj）など様々な言語によるポーティングが公開されています。

Goの場合は冒頭でも紹介したeverystreet/go-projがそれにあたります。

PROJはC++で実装されており、cgo経由で扱うため、PROJ環境を構築した後に、go-projから呼び出す必要があります。

go-projの利用開始

WindowsとWSL2での構築手順をまとめます。Macは手元に無かったので割愛します。

全体の注意ですが、2022.11.20時点でPROJの最新バージョンは 9.1.0 ですが、go-projは 8.1.0 に対応しています。8.1系だと最新が 8.1.1 が存在するので、この記事ではそのバージョンを利用します。

PROJのインストールは公式ドキュメント に細かく手順が載っていますが、conda 経由で行います（WSL2はapt-getも追記しておきます）。

各環境からは次のようなgo-projを用いたコードが動作することを確かめます。

main.go
package main

import (
	"fmt"
	"github.com/everystreet/go-proj/v8/proj"
	"github.com/golang/geo/s1"
)

func main() {
	wgsLng, wgsLat := Tky2Wgs(128542740/float64(60*60*256), 32706756/float64(60*60*256))
	fmt.Printf("%f %f\n", wgsLng, wgsLat) // 139.474598 35.492354
}

func Tky2Wgs(lng, lat float64) (float64, float64) {
	coord := proj.LP{
		Lng: s1.Angle(lng),
		Lat: s1.Angle(lat),
	}

	err := proj.CRSToCRS(
		"EPSG:4301", // 日本測地系(TOKYO) 緯度経度
		"EPSG:4326", // 世界測地系(WGS84) 緯度経度
		func(pj proj.Projection) {
			proj.TransformForward(pj, &coord)
		})
	if err != nil {
		panic(err) // サンプルコードなのでpanicにしていますが、errorを戻り値にした方が良いです
	}
	return float64(coord.Lng), float64(coord.Lat)
}

Windows環境構築手順

conda は miniconda を利用します。環境がない方は以下からインストールください。

もし、プロキシ環境で構築する場合は %USERPROFILE% 直下に .condarc ファイルを作成し、プロキシ情報を追記ください。

.condarc
proxy_servers:
    http: http://{user}:{password}@proxy.example.co.jp:8000
    https: http://{user}:{password}@proxy.example.co.jp:8000

※ {password} に記号が入っている人で、上記で認証が通らない方は、URLエンコードして設定してみてください

conda環境ができたら、インストールします。

> conda install -c conda-forge proj=8.2.1 proj-data=1.11.0

成功すれば以下のようなバージョンが表示されると思います。インストールが成功しているけどprojコマンドが動かない人は、 %USERPROFILE%\Miniconda3\Library\bin をPATHに追加してください。

>proj
Rel. 8.2.1, January 1st, 2022
usage: proj [-bdeEfiIlmorsStTvVwW [args]] [+opt[=arg] ...] [file ...]

次に環境変数を設定し、go get でパッケージをインストールします。

# 環境変数
> set CGO_CFLAGS=-I%USERPROFILE%\Miniconda3\Library\include
> set CGO_LDFLAGS=-L%USERPROFILE%\Miniconda3\Library\lib

# go-projの取得
> go get -u github.com/everystreet/go-proj/v8

これで動くと思います。

>go run main.go
139.474598 35.492354

もし、以下のようなエラーが出てのであれば、Cのヘッダーファイルが見つかっていないということですので、 CGO_CFLAGS で指定している、includeフォルダのパスを確認し再設定ください。

失敗例1
# github.com/everystreet/go-proj/v8/cproj
~\go\pkg\mod\github.com\everystreet\go-proj\v8@v8.0.0\cproj\cgo_helpers.go:8:10: fatal error: proj.h: No such file or directory
 #include "proj.h"
          ^~~~~~~~
compilation terminated.

あるいは、以下のエラーの場合は、 CGO_LDFLAGS で設定したライブラリのパスを確認し再設定ください。

失敗例2
~\go\go1.19.1\pkg\tool\windows_amd64\link.exe: running gcc failed: exit status 1
C:/Program Files/mingw64/bin/../lib/gcc/x86_64-w64-mingw32/8.1.0/../../../../x86_64-w64-mingw32/bin/ld.exe: cannot find -lproj
...
collect2.exe: error: ld returned 1 exit status

WSL2インストール手順

condaを使う場合

ほぼWindows側と同じです。miniconda3 をインストールし、PROJをインストール＆環境変数を設定しておしまいです。

# minicondaインストール（プロキシ設定は割愛）
$ wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
$ bash Miniconda3-latest-Linux-x86_64.sh

# PATHを通す（この例では .bashrcに export PATH=$PATH:~/miniconda3/bin を追記）
$  ~/.bashrc
$ source ~/.bashrc

# PROJインストール
$ conda install -c conda-forge proj=8.2.1 proj-data=1.11.0

# CGOのパラメータ設定
$ export CGO_CFLAGS=-I${HOME}/miniconda3/include
$ export CGO_LDFLAGS=-L${HOME}/miniconda3/lib

# go-projの取得
> go get -u github.com/everystreet/go-proj/v8

これでgo-projが実行できるようになると思います。

$ go run main.go
139.474598 35.492354

apt を使う場合

conda を使いたくない方向けに手順を残しておきます。まず、公式のパッケージリポジトリでインストールできるprojが6系と古いため、3rdパーティのパッケージリポジトリを参照する必要があります。以下から利用するリポジトリを選択し、運用していただいていることに感謝します。

私は以下のように、 proj.list というファイルを作って追加しました。

$ sudo vim /etc/apt/sources.list.d/proj.list
deb http://kr.archive.ubuntu.com/ubuntu jammy main universe

これでPROJをインストールします。

# アップデート
$ sudo apt update

# PROJ 8.2.1-1 が存在することを確認します
# $ sudo apt-cache show libproj-dev

# インストール
$ sudo apt install proj-bin=8.2.1-1 libproj-dev=8.2.1-1

これでgo-projが実行できるようになると思います。

$ go run main.go
139.474598 35.492354

apt経由の場合は、libproj-dev でインストールすると /usr/include にヘッダファイルなどがインストールされるため、CGO系のフラグは設定しなくても動くと思います。

高速化したい場合

proj.CRSToCRS() ですが、処理件数が多い処理に利用すると物足りない性能でした（1変換に10msほどかかりました）。理由は、C++側のライブラリに測地系のロードを都度行うためだと思われます。そのため、proj.CRSToCRS() をそのまま利用するのでなく、その內部で利用している cproj を直接利用すると良いかもしれません。

例えば以下のような実装です（※ちゃんとメモリリークしないか確かめていないです）。これで1000倍くらい早くなりました。

改良版
package main

import (
	"fmt"
	"github.com/everystreet/go-proj/v8/cproj"
	"github.com/everystreet/go-proj/v8/proj"
	"github.com/golang/geo/s1"
)

func main() {
	pj := NewProjection("EPSG:4301", "EPSG:4326")
	defer pj.Close()

	for i := 0; i < 1000*1000; i++ {
		wgsLng, wgsLat := pj.CRSToCRS(128542740/float64(60*60*256), 32706756/float64(60*60*256))
		fmt.Printf("%f %f\n", wgsLng, wgsLat)
	}
}

type Projection struct {
	ctx        *cproj.PJ_CONTEXT
	src        *cproj.PJ
	dst        *cproj.PJ
	pj         *cproj.PJ
	normalized *cproj.PJ
}

func NewProjection(source, target string) Projection {
	ctx := cproj.Context_create()
	src := cproj.Create(ctx, source)
	dst := cproj.Create(ctx, target)
	pj := cproj.Create_crs_to_crs_from_pj(ctx, src, dst, nil, nil)
	normalized := cproj.Normalize_for_visualization(ctx, pj)
	return Projection{
		ctx:        ctx,
		src:        src,
		dst:        dst,
		pj:         pj,
		normalized: normalized,
	}
}

func (p Projection) CRSToCRS(lng, lat float64) (float64, float64) {
	coord := proj.LP{
		Lng: s1.Angle(lng),
		Lat: s1.Angle(lat),
	}
	proj.TransformForward(p.normalized, &coord)
	return float64(coord.Lng), float64(coord.Lat)
}

func (p Projection) Close() {
	cproj.Context_destroy(p.ctx)
	cproj.Destroy(p.src)
	cproj.Destroy(p.dst)
	cproj.Destroy(p.pj)
	cproj.Destroy(p.normalized)
}

Close() 内でじゃかじゃか Destroy() しているのはC側のオブジェクトはGo側でGCされないためで、自分でリソースを開放する必要があるためです。

まとめ

GoからPROJを利用する方法をまとめました。 構築は cgo を使ったことがある方であればハマることは少ないと思いますが、conda経由だとCGOオプションが必要です。

測地系の変換を行うケースは一般的にはあまり多くないと思いますが、何かの参考になれば幸いです。最後まで読んでいただき、ありがとうございました！