Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
18
Help us understand the problem. What is going on with this article?
@MALORGIS

地図タイルの計算まとめ

More than 3 years have passed since last update.

地図タイルとは?

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

タイル座標.jpg
拡大率(Level)Z/とXYで表される。

osm00.png|平射.png
© 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]}/` )
18
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
MALORGIS

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
18
Help us understand the problem. What is going on with this article?