地図タイルの計算まとめ #GIS

地図タイルとは？

地図の配信に使われる手法
画像やGeoJSON/Protocol Buffersとかで配信される
これの計算メモ。

拡大率(Level)Z/とXYで表される。

|
© OpenStreetMap contributors

左がGoogle Maps/OSM/地理院地図などインターネットで使われる地図の表示方法(EPSG:3785)
※実際地球は丸く(楕円形)正方形で表すと歪むので距離単位mはあくまで投影している図面上の距離や面積なので注意

EPSG:3785
Pseudo-Mercator/Spherical Mercator
900913 3587 54004 41001 102113 102100
http://epsg.io/3785

タイルのレベルの計算

タイルのレベルは拡大率 L0は全世界でL1,L2とレベルが上がると表示領域が拡大される。

世界を正方形で表す = 球体
座標値:メートル単位

半径:6378137の球
2×r×π= [円周長:] 40075016.6855784861531768177614

つまり正方形の世界の一辺の長さ=球の円周長でこれを256pxで表したのがレベル0。
レベルは1pxが何mかで表され、レベル0を4タイル(1タイルは256pxの正方形)で表現したのがレベル1(Quadtree)。

[円周] / 256pxで約15,654メートル/px

L08で1px=611.49メートルくらい
L17で1px= 1.19メートルくらい

これを求める計算

const GEO_R:number = 6378137;
let level:number = 15

window.document.write(
  (2 * GEO_R * Math.PI / 256 / Math.pow(2, level)).toString()
);

緯度経度からタイルXYを求める

グリニッジ子午線-赤道がxy:00なので円周長を半分に割った値が北西端(左上端)
[円周長:] 40075016.6855784861531768177614 / 2

-20037508.3427....
20037508.3427....

レベルが指定された場合、
左上端を原点としてタイルは始まるので256px=何メートルがあれば原点XYとのm差で割るだけで行番号(y)列番号(x)が求まる。

緯度経度とm単位の換算は探せばいくらでも出てくる。
degrees2meters.js
https://gist.github.com/onderaltintas/6649521
緯度方向が高緯度になるほどm単位が加算される

let degrees2meters = function(lat:number, lon:number) {
  var x = lon * 20037508.34 / 180.0;
  var y = Math.log(Math.tan((90.0 + lat) * Math.PI / 360.0)) / (Math.PI / 180.0);
  y = y * 20037508.34 / 180.0;
  return [x, y]
}

const GEO_R:number = 6378137;
const orgX = -1 * (2 * GEO_R * Math.PI / 2);
const orgY = (2 * GEO_R * Math.PI / 2);

//let xy = degrees2meters(36.104600,140.085871);
//let level =17;
let xy = degrees2meters(35, 135);
let level =15;

let unit = 2 * GEO_R * Math.PI / Math.pow(2, level)

let xtile = Math.floor((xy[0] - orgX) / unit);
let ytile = Math.floor((orgY - xy[1]) / unit);

//https://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png

window.document.write( `https://cyberjapandata.gsi.go.jp/xyz/std/${level}/${xtile}/${ytile}.png` )

都合こんな感じで35,135辺りのタイルを参照できる。

タイルXYZから緯度経度を求める

タイルXYZから緯度経度(北西/左上起点)を求める。(unit足し引きで4隅)

let meters2degress = function(x,y) {
  let lon = x *  180 / 20037508.34 ;
  //thanks magichim @ github for the correction
  let lat = Math.atan(Math.exp(y * Math.PI / 20037508.34)) * 360 / Math.PI - 90; 
  return [lon, lat]
}


let level = 15;
let xtile = 28671;
let ytile = 12979;

const GEO_R:number = 6378137;
const orgX = -1 * (2 * GEO_R * Math.PI / 2);
const orgY = (2 * GEO_R * Math.PI / 2);

let unit = 2 * GEO_R * Math.PI / Math.pow(2, level)

let x = orgX + xtile * unit;
let y = orgY - ytile * unit;

let latlon = meters2degress(x,y);

window.document.write( `https://maps.gsi.go.jp/#${level}/${latlon[1]}/${latlon[0]}/` )