mottie0911’s blog

技術や北海道について

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)が求められた

OSM の建物のデータをMySQL8系にImportしてQGISで表示してみた

こちらはQiitaに投稿した記事になります。


この記事は2019年12月14日に行われたFoss4Gの"GeoSaturday"にて私mottie0911がTryした内容になります

中堅WebエンジニアでQGIS初心者ですが、MySQL8 でGIS機能が使えると聞いてOSM(OpenStreetMap)のデータをQGISで表示することに挑戦してみました

ogr2ogr の変換や全体的な手法については以下を参考にしました

https://github.com/YoshiakiYamasaki/MySQL-GIS-Data-Japan-eStat

環境

QGISWindowsで使用しているので、WindowsにMySQL8系をInstallして実行した

OSM データのDL

OSMのデータをDLできるサイトHOTから、今回はShapeFileを選択しDLした(※要OSMの会員登録 & ログイン) 日本全国すべてのデータだと大きいので今回は札幌周辺の建物データとする。


1. DLしたい地域を選択(今回は札幌周辺をBOX選択) "1.Describe" タブの中のName等を適当に記載(今回のNameはsapporo_test)   

  1. "2 Formats" タブの中からDLしたいフォーマットを選択(今回はShapefile)

  2. "3 Data" タブからBuildingsを選択   ※このサイトはyaml等でDataの抽出方法を設定できるのでここの辺を使いこなせたら便利そう

  3. "4 Summary" タブより内容を確認し "Create Export"を押下。しばらく経つとダウンロードできるようになる。   (今回は2分程度、データ作成処理が完了すると完了通知メールが登録メールアドレス宛に届く)

ogr2ogr 変換

DLしzipを展開するとShapefile一式があるはずなので、.shpファイルに対しogr2ogrで変換&MySQL取り込みを行う
下準備として接続したいMySQLサーバにデータを格納したいDBを作成しておく。今回は以下の通り。

> create database gis_test;


データベースができたら、DL&回答したディレクトリに移動し、WindowsのPoworShell等より以下のコマンドを実行

> ogr2ogr -f "ESRI Shapefile" -lco ENCODING=UTF-8 -oo  .\sapporo_test_planet_osm_polygon_polygons_utf8.shp .\sapporo_test_planet_osm_polygon_polygons.shp


  • MySQL取り込み(今回は同じPC内にあるMySQLサーバに対し接続を行っている、うまく接続できない場合はHost, Portなどのパラメータを増やすこと)
> ogr2ogr -f "MySQL" MySQL:"gis_test,user=root,password=【パスワード】" .\sapporo_test_planet_osm_polygon_polygons_utf8.shp


mysql> use gis_test
Database changed
mysql> show tables;
+-----------------------------------------------+
| Tables_in_gis_test                            |
+-----------------------------------------------+
| geometry_columns                              |
| sapporo_test_planet_osm_polygon_polygons_utf8 |
+-----------------------------------------------+
2 rows in set (0.01 sec)

mysql> desc sapporo_test_planet_osm_polygon_polygons_utf8;
+------------+--------------+------+-----+---------+----------------+
| Field      | Type         | Null | Key | Default | Extra          |
+------------+--------------+------+-----+---------+----------------+
| OGR_FID    | int(11)      | NO   | PRI | NULL    | auto_increment |
| SHAPE      | geometry     | NO   | MUL | NULL    |                |
| osm_id     | double       | YES  |     | NULL    |                |
| building   | varchar(80)  | YES  |     | NULL    |                |
| addrhousen | varchar(137) | YES  |     | NULL    |                |
| addrstreet | varchar(133) | YES  |     | NULL    |                |
| buildingma | varchar(80)  | YES  |     | NULL    |                |
| name       | varchar(183) | YES  |     | NULL    |                |
| accessroof | varchar(80)  | YES  |     | NULL    |                |
| roofmateri | varchar(80)  | YES  |     | NULL    |                |
+------------+--------------+------+-----+---------+----------------+
10 rows in set (0.01 sec)

どうやら入っている模様

QGIS側のMySQLプラグインのInstall

QGISからMySQLにつなぐには"MySQL Importer"というプラグインが必要 "QGISプラグインプラグインの管理とインストール→すべてのプラグイン"から "MySQL Importer"を検索しインストールすればよいのだが、そのままだと起動に失敗する。


回避方法は以下のサイトにのっているので実施し、Uninstall&Install し、有効化しておく How to Install and Enable MySQL Importer Plugin on QGIS 3.4

今回DLしたPythonのMysqlClientは "mysqlclient‑1.4.6‑cp37‑cp37m‑win_amd64.whl"

whlで圧縮されたファイルから対象フォルダを抽出しなければならないが、7zipで場所を指定して展開、等ではなく".whlファイルを右クリック→一番上の展開"より7zipのブラウザを起動し、フォルダに対してコピーをかけると上手く"C:\Program Files\QGIS 3.4\apps\Python37\Lib\site-packages"配下に対象ディレクトリ(mysqlclient-1.4.6.dist-info, MySQLdb)を持っていくことができる

QGISで表示

QGISの "レイヤ→データソースマネージャ" のベクタよりデータベースを選択

"MySQL Importer" がインストールされれば、データベースのタイプをプルダウンからMySQLが選択できるはずなのでMySQLを選択

接続の部分より新規の接続を作成(今回認証はベーシック認証を使用しました)

  • DBのユーザとPasswordを入れて構成に変換

  • 変換された認証設定を使用し接続テストを行い接続に成功すればOK

  • 接続設定ができたら"追加"を押下、追加するテーブル名を聞かれるので "sapporo_test_planet_osm_polygon_polygons_utf8" を選択

  • ちゃんと表示できました

余談

Qiitaにメモ程度以外の記事を書くのは初めてなのでこの記事はQiitaのお作法にのっとってないかもしれません。

Web系エンジニアはPosgreではなくMySQLの方が馴染み深いと思うので、そういう人にとってGISを始めるハードルが下がれば幸いです。 (MySQLで接続できちゃえば、例えば建物のnameをとってきたりUpdateするのは簡単ですね)

追記 (MySQL workbenchでの内容確認)

yymazaki1さんより、「そのデータMySQL workbenchでもみれるよ」との情報を頂いたので今回のデータをMySQL 8.0にシェープファイルを取り込んで、MySQL Workbenchでポリゴンの形を確認するの記事にそって試してみました。もしQGISをインストールしていなくても、こちらで建物の図形や情報を確認することができます。

尚、MySQL workbenchはMySQL8に同封されてインストールされていました
メニューの"Database→Connect to Database"より接続 mysqlb1.png


左側のTabよりSchemasを選択 mysqlb2.png


テーブルを右クリックで選択し、"Select Rows -- Limit 1000" を実行 mysqlb3.png


実行され、画面中央に結果が表示されました mysqlb4.png


select文を実行したり、FormEditorを使用すればgeometry型(POLYGON)の値や図を確認することができる mysqlb5.png

yymazaki1さん、ありがとうございました