PROJ.4 - Windows (VisualStudio) 向けのインストール方法と、C++での緯度経度→Webメルカトルのピクセル座標変換

最近PROJ.4を触る機会があり、日本語の情報があまりない(特にWindows用、Python/jsで使ってみた以外)のでまとめてみる

PROJ.4とは

PROJ is a generic coordinate transformation software that transforms geospatial coordinates from one coordinate reference system (CRS) to another. This includes cartographic projections as well as geodetic transformations. PROJ is released under the X/MIT open source license

公式サイトより

つまりFromとToの座標参照システム(※)を指定すると、変換してくれるオープンソースのソフトウェアライブラリのこと

ライブラリ自体はC/C++で書かれていて、インストールにはC/C++11のコンパイラが必要になる。またコンパイルされたものはコマンドラインインターフェースで処理を受け付けるAPIも含まれているらしい…が今回はAPIは使用せずコンパイルして生成されたヘッダファイルを使用していきます

※ 座標参照システム(CRS)とは、測地系・座標系等の定義のことで、その定義にはOGPが管理しているEPSGコードが割り当てられている(GPSで使われているWGS84はEPSG:4326, Google Map API等で使われているWebメルカトルはEPSG:3857)。PROJ.4ではこのEPSGコードごとの定義をコードに組み込み使用していく

EPSGコードの説明は以下のサイトが詳しい

GeoDjangoではじめる地理空間情報: 測地系と座標系

DL & Windows10, Visual Studio 2017 を用いてコンパイル

公式サイトより、今回は"proj-5.2.0.tar.gz"をDLし、7zip等で解凍する

(5系と6系でコンパイル方法が変わっているので注意)

解凍した"proj-5.2.0"フォルダ直下にある"makefile.vc"を使っていく

ここで脳死でREADMEに書いてあるとおりnmakeを打っていくとstdio.hがない等エラーがでてしまうので、VS2017インストール時に付随している NativeTools コマンドプロンプトを使用する

(私はこれがわからず、脳死でPowerShellでコマンドを叩いたり、VSのPathをシステム環境設定に追加する等無駄足を踏んでしまった)

あとはREADMEにある通り、以下のコマンドを実行していく

C:\>cd C:\tools\proj-5.2.0

C:\tools\proj-5.2.0>nmake -f makefile.vc

C:\tools\proj-5.2.0>nmake -f makefile.vc install-all

実行すると"C:\tools\"で実行した場合 "C:\tools\PROJ\"が生成されているはず

C/C++のプログラムで使うには"C:\tools\PROJ\lib\"の"proj.lib", "proj_i.lib"と"C:\tools\PROJ\include\"中にある "proj.h" や "proj_api.h"を使用すればOK

このときVS2017上のプロジェクトのプロパティより"[C/C++]→[追加のインクルードディレクトリ]"に"C:\tools\PROJ\include\"、”[リンカー]→[全般]→[追加のライブラリディレクトリ]”に"C:\tools\PROJ\lib\、"[リンカー]→[入力]→[追加の依存ファイル]"に"proj.lib;proj_i.lib;"の追加を忘れないこと

あとは"proj.h"や"proj_api.h"を使用して、コードを書けば良い

C++での緯度経度→メルカトル座標の投影値変換

公式サイトのPROJ v. 5のサンプルより

※入力値受け入れの部分だけ改変

#include <iostream>
#include <proj.h>


int main(int argc, char **argv) {
    PJ *P;
    PJ_COORD c;

    P = proj_create(PJ_DEFAULT_CTX, "+proj=merc +ellps=clrk66 +lat_ts=33");
    if (P == 0)
        return 1;

    c.lp.lam = 141.242;
    c.lp.phi =  45.1785;
    c.lp.lam = proj_torad(c.lp.lam);
    c.lp.phi = proj_torad(c.lp.phi);
    c = proj_trans(P, PJ_FWD, c);
    printf("%.2f\t%.2f\n", c.xy.x, c.xy.y);
    

    proj_destroy(P);
}

この

"+proj=merc +ellps=clrk66 +lat_ts=33"

の部分がPROJ.4に与える、出力したい座標参照システムの定義

この場合はメルカトル図法で投影座標に変換している

これはこちらのページのProj4のリンクより定義式を確認できる

+proj の部分は投影方法(メルカトルならmerc、緯度経度ならlonglat、UTM座標系ならutm)、+ellpsの部分は楕円体の扱い方、それ以降の+lat_ts等で投影座標の中心を指定できる

冒頭に"FromとToの座標参照システムを指定"と書いたが、ver5以上では入力はデフォルトで緯度経度で扱われている模様

詳細はwikibooksのPROJ.4のページで確認できる

C++での緯度経度→WEBメルカトルのピクセル座標への変換

まずはWebメルカトル図法(EPSG:3857)のProj4の定義を確認

SR-ORG:7483より

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs

であることがわかる

だが、こちらの定義式を使った、zoomレベルを考慮したWebメルカトルのピクセル座標の出し方がわからなかった

検索した結果

proj.4で経緯度からいきなりWebメルカトルのピクセル座標やタイル座標を出す方法

こちらより、以下のproj4定義を使用した

cpp +proj=merc +a={Math.pow(2,z+7) / Math.PI} +b={Math.pow(2,z+7) / Math.PI} +lat_ts=0.0 +lon_0=0.0 +x_0={Math.pow(2,z+7)} +y_0=-{Math.pow(2,z+7)} +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

{}の中をZoomレベルによって定数化したのは以下の通り

zoom_level 0:      +proj=merc +a=40 +b=40 +lat_ts=0.0 +lon_0=0.0 +x_0=128 +y_0=-128 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 1:      +proj=merc +a=81 +b=81 +lat_ts=0.0 +lon_0=0.0 +x_0=256 +y_0=-256 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 2:      +proj=merc +a=162 +b=162 +lat_ts=0.0 +lon_0=0.0 +x_0=512 +y_0=-512 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 3:      +proj=merc +a=325 +b=325 +lat_ts=0.0 +lon_0=0.0 +x_0=1024 +y_0=-1024 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 4:      +proj=merc +a=651 +b=651 +lat_ts=0.0 +lon_0=0.0 +x_0=2048 +y_0=-2048 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 5:      +proj=merc +a=1303 +b=1303 +lat_ts=0.0 +lon_0=0.0 +x_0=4096 +y_0=-4096 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 6:      +proj=merc +a=2607 +b=2607 +lat_ts=0.0 +lon_0=0.0 +x_0=8192 +y_0=-8192 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 7:      +proj=merc +a=5215 +b=5215 +lat_ts=0.0 +lon_0=0.0 +x_0=16384 +y_0=-16384 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs
z
oom_level 8:      +proj=merc +a=10430 +b=10430 +lat_ts=0.0 +lon_0=0.0 +x_0=32768 +y_0=-32768 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 9:      +proj=merc +a=20860 +b=20860 +lat_ts=0.0 +lon_0=0.0 +x_0=65536 +y_0=-65536 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 10:     +proj=merc +a=41721 +b=41721 +lat_ts=0.0 +lon_0=0.0 +x_0=131072 +y_0=-131072 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 11:     +proj=merc +a=83443 +b=83443 +lat_ts=0.0 +lon_0=0.0 +x_0=262144 +y_0=-262144 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 12:     +proj=merc +a=166886 +b=166886 +lat_ts=0.0 +lon_0=0.0 +x_0=524288 +y_0=-524288 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 13:     +proj=merc +a=333772 +b=333772 +lat_ts=0.0 +lon_0=0.0 +x_0=1048576 +y_0=-1048576 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 14:     +proj=merc +a=667544 +b=667544 +lat_ts=0.0 +lon_0=0.0 +x_0=2097152 +y_0=-2097152 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 15:     +proj=merc +a=1335088 +b=1335088 +lat_ts=0.0 +lon_0=0.0 +x_0=4194304 +y_0=-4194304 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 16:     +proj=merc +a=2670176 +b=2670176 +lat_ts=0.0 +lon_0=0.0 +x_0=8388608 +y_0=-8388608 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 17:     +proj=merc +a=5340353 +b=5340353 +lat_ts=0.0 +lon_0=0.0 +x_0=16777216 +y_0=-16777216 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 18:     +proj=merc +a=10680707 +b=10680707 +lat_ts=0.0 +lon_0=0.0 +x_0=33554432 +y_0=-33554432 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

zoom_level 19:     +proj=merc +a=21361414 +b=21361414 +lat_ts=0.0 +lon_0=0.0 +x_0=67108864 +y_0=-67108864 +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs

検算

trail-noteさんの座標の変換（世界座標、ピクセル座標、タイル座標、緯度・経度）の"緯度、経度からピクセル座標への変換"より以下のコードに代入して結果を比較する

#define _USE_MATH_DEFINES

#include "pch.h"
#include <iostream>
#include <proj.h>
#include <string>
#include <tuple>

using namespace std;

tuple<double, double > get_pixel(int, double, double);


int main(int argc, char** argv)
{

    double lat = atof(argv[1]);
    double lng = atof(argv[2]);
    int zoom = atoi(argv[3]);



    double result_pix_x, result_pix_y;
    std::tie(result_pix_x, result_pix_y) = get_pixel(zoom, lat, lng);

    
    printf("result_pix_x:%.2f\tresult_pix_y:%.2f\n", result_pix_x, result_pix_y);

    return 0;
}

std::tuple<double, double> get_pixel(int zoom_level, double x, double y) {

    PJ *P;
    PJ_COORD c;

    long double value1 = pow(2, zoom_level + 7) / M_PI;
    long double value2 = pow(2, zoom_level + 7);

    std::string value1_s = std::to_string((value1));
    std::string value2_s = std::to_string((value2));

    std::string zoom_calc = "+proj=merc +a=" + value1_s + " +b=" + value1_s + " +lat_ts=0.0 +lon_0=0.0 +x_0=" + value2_s + " +y_0=-" + value2_s + " +k=1.0 +axis=esu +units=m +nadgrids=@null +no_defs";


    P = proj_create(PJ_DEFAULT_CTX, zoom_calc.c_str());

    c.lp.lam = y;
    c.lp.phi = x;

    c.lp.lam = proj_torad(c.lp.lam);
    c.lp.phi = proj_torad(c.lp.phi);
    c = proj_trans(P, PJ_FWD, c);

    return forward_as_tuple(floor(c.xy.x), floor(c.xy.y));

}

結果、記事に掲載されているZoomレベルが17のときのピクセル座標(x,y)が求められた

mottie0911’s blog

技術や北海道について

PROJ.4 - Windows (VisualStudio) 向けのインストール方法と、C++での緯度経度→Webメルカトルのピクセル座標変換

PROJ.4とは

DL & Windows10, Visual Studio 2017 を用いてコンパイル

C++での緯度経度→メルカトル座標の投影値変換

C++での緯度経度→WEBメルカトルのピクセル座標への変換

検算