オンライン地図のタイル座標

　最近ベースマップに使っているOpenStreetMapが読み込めなくなっていた。

　タイル利用規約には
　『データは誰でも自由に使えます。我々のタイルサーバは違います。』
とある。

　多用しているものの単発的な利用でヘビーユースには該当しないはずだが……。

　元に戻ったようだから、まあいいか。

目　次

タイル座標について
1. 緯度、経度　⇒　タイル座標
Link

タイル座標について

　タイルレイヤーの使い方については下のページ。

Leafletのレイヤー

レイヤーのコントロール。Leafletのサイト参照　<div id = "map"></div><script>　の下　var map = L.map('map').setView([35.30, 138.30] , 5);

　これまでLeafletで使ってきたタイルはpng、jpgの画像タイル。

　一方、txt、geojson、pbf（pdfではない）など画像ではないタイルもある。
　こちらはプログラミングできないと使いこなせないっぽい。

　追）否、何とかなるかも

ベクトルタイルをいじってみる

梅雨明けはまだ。本来、今頃五輪の最中。懸念されていた猛暑も来年に延期か。さて、前々回少し触れた「画像ではないタイル」のうちベクトルタイルについて、いじってみた。拡張子がgeojson、pbf（pdfではない）のタイル

　広範でシームレスに利用するのでなければ、タイル座標を入力して必要な範囲のデータを使うこともできる。

　タイル座標とは、https://tile.openstreetmap.jp/{z}/{x}/{y}.png

の{z}/{x}/{y}の部分。

　x、yは各々経度、緯度に該当するが、変換されたピクセル値　pixelから割り出した値。
　zは高さには違いないが、ズーム値。小さい値ほど小縮尺（分母数が大の縮尺）の地図。

　国土地理院　タイル座標確認ページ（maps.gsi.go.jp/development/tileCoordCheck.html）
　産業技術総合研究所　タイル座標確認ツール（gsj-seamless.jp/labs/tools/tileCoord.html）

などから確認することができる。
　地理院地図はズーム18まで。

　ズーム０は世界全図。
　正方形のタイル１つにおさめられている。

　タイルの一辺を赤道一周約４万kmとして、地点をメートル　mで表した座標系が以前ちょっととりあげたEPSG3857：WGS84　球面メルカトル図法。

GIS　CRSとEPSG

GADMデータから各国の行政界を描画している途中、サン・マリノの地図がベースマップ（OpenStreetMap）とずれているのが気になった。簡素化しすぎて合わないというのではなく、全体的に西にずれている

　一辺を256ピクセル　pixelとして、地点をピクセル値で表した座標がタイル座標。１ピクセルあたり約４万／256＝約156km。
　ただし、x＝0は経度0°ではなく経度180°、y＝0は北緯およそ85°。

　ズーム１は256×256のタイルが２倍×２倍で計４タイル。
　ズーム２は（２×２倍）×（２×２倍）で計16タイル。
　ズーム３は（２×２×２（＝２³）倍）×２³倍で計64タイル。
　見る機会が多いのはズーム15の25000分の１地形図だが、
　ズーム15は２¹⁵×２¹⁵で計２³⁰（＝1073741824）タイル。

緯度、経度　⇒　タイル座標

　経度、緯度からタイル座標への変換は、
　経度の場合、
　（経度＋180）×256／360
でズーム０におけるピクセル値（ただし０～256の実数）が求まり、
　２のズーム乗（z乗）を掛けて256で割った値の整数部分が求めるタイルx座標。
　つまり、ズーム値　zのピクセル値は、
　（経度＋180）／360 × ２^z

　緯度の場合、等分でないことに加えて、南北85.0511°までの範囲なので、すんなりいかない。

　逆双曲線関数との遭遇。球面幾何学の世界。

　頭がねじ切れてしまうので、Slippy map tilenamesのページ（wiki.openstreetmap.org/wiki/Slippy_map_tilenames）を参考に計算。

　Ｒ言語で、経度　lon、緯度　latからタイル座標　tilex、tileyへ変換。
　tilexyzという関数（関数名は任意）を作って、値を入力。

tilexyz <- function(lon,lat,zoom){
tilex <- (lon + 180) / 360 * 2^zoom
tiley <- (1 - log(tan(lat * pi / 180)+ (1/cos(lat * pi / 180))) / pi) * 2^(zoom-1)
return(list(tilex,tiley,zoom))
}

　例えばズーム15の箱根・強羅の25000分の１地形図であれば