位置と空間アドレス

MIRGRID では、地球上の場所をセルという立方体の単位で扱います。この文書では、セルとレベルの考え方と、緯度・経度や平面直角座標系の座標をセルや空間アドレスに変換する方法を説明します。

以下の例は、最初に次のように読み込んだものとして書いています。Node.js で試すときは、init() の呼び方をはじめにのとおりに変えてください。

import init, {
  spatial_hilbert_cellid_decode,
  spatial_voxel_to_ecef,
  spatial_voxel_bounds,
  spatial_cellid_to_ecef,
  spatial_ecef_to_llh_rad,
  spatial_ecef_height,
  spatial_ecef_to_enu,
  spatial_enu_to_ecef,
  spatial_ecef_to_enu_many,
  spatial_llh_to_plane_rect,
  spatial_plane_rect_to_llh_deg,
  spatial_plane_rect_to_ecef,
  spatial_ecef_to_plane_rect,
  spatial_plane_rect_scale_factor,
} from '@mirgrid/web-sdk';

await init();

セルとレベル

セルの位置は ECEF をもとに決まっています。ECEF は、地球の中心を原点にした x・y・z の直交座標で、単位は m です。地球をすっぽり包む大きな立方体を考え、それを各辺 2 つずつに分けていったものがセルです。分けた回数をレベルと呼びます。レベルは 1〜30 で、1 つ上がるごとにセルの一辺が半分になります。

セルの一辺は spatial_voxel_size_m(level) で求められます。レベルごとの値は次のとおりです（単位は m）。

レベル	一辺	レベル	一辺
1	8388608	16	256
2	4194304	17	128
3	2097152	18	64
4	1048576	19	32
5	524288	20	16
6	262144	21	8
7	131072	22	4
8	65536	23	2
9	32768	24	1
10	16384	25	0.5
11	8192	26	0.25
12	4096	27	0.125
13	2048	28	0.0625
14	1024	29	0.03125
15	512	30	0.015625

レベル 24 でちょうど一辺 1 m になります。

それぞれのセルは、x・y・z の 3 つの整数の番号（index）で区別します。番号は 0 から、2 を「レベル」回掛けた数より 1 小さい数までです。たとえばレベル 24 なら 0〜16777215 です。

セルの面は ECEF の軸に沿っています。ECEF の軸は地面の東西南北とは揃っていないので、地上から見るとセルは少し傾いた立方体になります。地面に水平な格子が欲しいときは、ファイルをセルにして表示するで説明する、現場に揃えた格子を使ってください。

空間アドレス

空間アドレスは、セルの番号とレベルを 1 つの文字列にまとめた ID です。32 桁の 16 進で書きます。番号はヒルベルト曲線の順に振っているので、近くにあるセルどうしは、空間アドレスの前のほうの桁が揃いやすくなります。

空間アドレスから番号とレベルに戻すには、spatial_hilbert_cellid_decode を使います。

const cell = spatial_hilbert_cellid_decode('0000000002e9380374bc9e7c5e140000');
console.log(cell);

{ vx: 4428917, vy: 11738702, vz: 12088151, level: 24 }

逆に、番号とレベルから空間アドレスを作るには spatial_hilbert_cellid_hex(vx, vy, vz, level) を使います。緯度・経度から空間アドレスを求める手順は、はじめにの例を見てください。

セルの中心と境界

セルの中心は spatial_voxel_to_ecef、境界は spatial_voxel_bounds で求めます。どちらも ECEF の値を返します。境界は [最小の x, 最小の y, 最小の z, 最大の x, 最大の y, 最大の z] の順で、最小の側はセルに含み、最大の側は含みません。

const center = spatial_voxel_to_ecef(cell.vx, cell.vy, cell.vz, cell.level);
const bounds = spatial_voxel_bounds(cell.vx, cell.vy, cell.vz, cell.level);
console.log('中心:', ...center);
console.log('最小:', ...bounds.subarray(0, 3));
console.log('最大:', ...bounds.subarray(3, 6));

中心: -3959690.5078125 3350094.4921875 3699543.4921875
最小: -3959691.0078125 3350093.9921875 3699542.9921875
最大: -3959690.0078125 3350094.9921875 3699543.9921875

空間アドレスから直接中心を求めるときは、spatial_cellid_to_ecef が使えます。

console.log(...spatial_cellid_to_ecef('0000000002e9380374bc9e7c5e140000'));

-3959690.5078125 3350094.4921875 3699543.4921875

ECEF から緯度・経度に戻す

ECEF の点を緯度・経度に戻すには spatial_ecef_to_llh_rad を使います。戻り値の lon と lat はラジアンなので、度にするときは 180 / Math.PI を掛けます。高さ（楕円体高）は spatial_ecef_height で別に求めます。

const { lon, lat } = spatial_ecef_to_llh_rad(center[0], center[1], center[2]);
const height = spatial_ecef_height(center[0], center[1], center[2]);
console.log(lon * 180 / Math.PI, lat * 180 / Math.PI, height);

139.7671232664105 35.68123605033351 40.203132109716535

はじめにで変換した点は、経度 139.767125・緯度 35.681236・高さ 40 m でした。セルの中心はその点を含むセルの真ん中なので、元の点とは少しずれます。ずれは最大でセルの対角線の半分です。

ENU に変換する

ENU は、ある地点を原点にして、そこから見た東・北・上を軸にした座標です。単位は m です。地図や 3D の画面にセルを描くときは、ECEF のままではなく ENU に変換してから使います。

ECEF の値は数百万 m の大きさがあります。そのまま Float32Array に入れて GPU に渡すと、数十 cm 刻みに丸められてしまいます。表示したい場所の近くに原点を置いた ENU に変換すれば、値が小さくなるので、この丸めを避けられます。

spatial_ecef_to_enu には、ECEF の点と、原点の経度・緯度・楕円体高を渡します。

const origin = [139.767125, 35.681236, 40];
const enu = spatial_ecef_to_enu(center[0], center[1], center[2], ...origin);
console.log(...enu);

-0.15693441620920887 0.00558469607271514 0.20313210846190108

セルの中心は、原点から西に約 0.16 m、北に約 0.006 m、上に約 0.20 m の位置にあることがわかります。

ENU から ECEF に戻すには spatial_enu_to_ecef を使います。

const east100 = spatial_enu_to_ecef(100, 50, 10, ...origin);
console.log(...east100);
console.log(...spatial_ecef_to_enu(east100[0], east100[1], east100[2], ...origin));

-3959739.0120321712 3350004.33505192 3699589.815651013
100.00000000017756 50.00000000030757 10.00000000015702

たくさんの点をまとめて変換するときは、spatial_ecef_to_enu_many と spatial_enu_to_ecef_many を使います。点を [x, y, z, x, y, z, …] のように並べた Float64Array を渡すと、同じ並びの結果が返ります。何十万点あっても、1 回の呼び出しで済みます。

const points = new Float64Array([...center, ...east100]);
console.log(...spatial_ecef_to_enu_many(points, ...origin));

-0.15693441620920887 0.00558469607271514 0.20313210846190108 100.00000000017756 50.00000000030757 10.00000000015702

平面直角座標系と変換する

平面直角座標系は、日本の測量や図面でよく使われる座標系です。日本を 19 の系に分け、それぞれの原点からの距離（m）で位置を表します。この SDK は、JGD2011 の平面直角座標系 I〜XIX 系（EPSG:6669〜6687）に対応しています。

関数の引数と戻り値では、東向きの座標を easting、北向きの座標を northing と呼びます。測量の慣習では北向きを X、東向きを Y と書くので、順番に気をつけてください。

緯度・経度と変換する

spatial_llh_to_plane_rect は、経度・緯度を [easting, northing] に変換します。spatial_plane_rect_to_llh_deg はその逆で、度の [経度, 緯度] を返します。東京駅は IX 系（EPSG:6677）に入ります。

const [easting, northing] = spatial_llh_to_plane_rect(139.767125, 35.681236, 6677);
console.log(easting, northing);
const [lonDeg, latDeg] = spatial_plane_rect_to_llh_deg(easting, northing, 6677);
console.log(lonDeg, latDeg);

-5992.919570423696 -35363.23774492909
139.76712500000005 35.68123599994409

複数の点をまとめて変換するときは、[経度, 緯度, 経度, 緯度, …] を渡す spatial_llh_to_plane_rect_many を使います。

ECEF と変換する

高さも含めて ECEF と行き来するときは、spatial_plane_rect_to_ecef と spatial_ecef_to_plane_rect を使います。高さは楕円体高です。

const fromPlane = spatial_plane_rect_to_ecef(easting, northing, 40, 6677);
console.log(...fromPlane);
console.log(...spatial_ecef_to_plane_rect(fromPlane[0], fromPlane[1], fromPlane[2], 6677));

-3959690.4857008546 3350094.26791384 3699543.369164228
-5992.919570420158 -35363.23775113177 40.00000000093132

対応していない EPSG を渡すと Error になります。

try {
  spatial_llh_to_plane_rect(139.767125, 35.681236, 4326);
} catch (error) {
  console.log(error.message);
}

unsupported EPSG (expected 6669..=6687)

縮尺係数

平面直角座標系の上で測った距離は、実際の距離と少し違います。その比が縮尺係数 k で、「座標の上の距離 = 実際の距離 × k」の関係があります。spatial_plane_rect_scale_factor で、その地点の k を求められます。

const k = spatial_plane_rect_scale_factor(easting, northing, 6677);
console.log(k, 100 / k);

0.9999004424246567 100.0099567488041

東京駅のあたりでは、座標の上で 100 m 離れた 2 点は、実際には約 100.01 m 離れています。

高さとジオイド高

この SDK の関数が受け取る高さと返す高さは、すべて楕円体高です。楕円体高は、地球の形を近似した回転楕円体の面から測った高さです。一方、図面や測量で使う標高は、ジオイド（平均海面にあたる面）から測った高さです。

2 つの高さの差をジオイド高 N と呼び、「楕円体高 = 標高 + N」の関係があります。手元の高さが標高のときは、N を足してから SDK に渡してください。

SDK には、ジオイド高を求める関数はありません。国土地理院の「ジオイド高計算」などで、その場所の N を調べてください。次の例は、静岡市の駿府城東御門橋のあたり（VIII 系、EPSG:6676）で、標高 22.45 m の点を ECEF にしています。この場所の N は 40.53 m です。

const elevation = 22.45; // 標高（m）
const geoidHeight = 40.53; // ジオイド高 N（m）
const bridge = spatial_plane_rect_to_ecef(-10492.5, -113402.1, elevation + geoidHeight, 6676);
console.log(spatial_ecef_height(bridge[0], bridge[1], bridge[2]));

62.980000000447035

点群や 3D モデルをセルにするときは、VoxelizerJs の入力の指定にジオイド高を書けます（入力の座標系を指定するを参照）。