LoginSignup
3
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@nokonoko_1203(のこのこ)inMIERUNE

4次元時空間情報利活用のための「空間ID」をRustで触ってみる！

Last updated at Posted at 2025-05-13

image.png

みなさんこんにちわ！！！

「空間ID」ってご存知でしょうか？？まぁ普通は知らないっすよね！

空間IDというのはいわゆる空間インデックスの一種であり、地理空間をボクセルで8分木分割し、各ボクセルに一意のIDを振ったものの総称になります。

仮想空間と現実空間をシームレスに融合させるいわゆる「デジタルツイン」やGISの文脈で近年利用され始めています。

類似の概念としては、Uberの「H3」やGoogleの「S2」などがありますが、空間IDはその中でも特に3次元を意識しています。

8分木のため、Zオーダーやヒルベルト曲線などの空間インデックスのアルゴリズムを用いることができるので「空間インデックス」という単語を用いましたが、単にDBで高速に検索するだけなら他のインデックスで十分なため、どちらかというと複数のシステムを連携させるための「共通キー」といった認識の方が良いとは思います。

詳しいことは、すべて4次元時空間情報利活用のための空間IDガイドライン（1.0版）に記載されており、基本的にはここを見れば全て解決するのですが、今回はその中から少しだけ掻い摘んで説明してみます。
（前半はガイドラインを眺めていく感じになります。）

空間IDの概要

空間IDは特定の領域を一意に識別するために生まれました。

3次元空間を8分木で格子状に再帰的に分割していくことで、さまざまな粒度のボクセルを定義できます。

分割された各種ボクセルは{z}/{f}/{x}/{y}という形式の「空間ID」が割り当てられます。

image.png
（4次元時空間情報利活用のための空間IDガイドライン（1.0版） 図2-1）

空間ボクセルは以下のような構造を持ちます。

  • 最上位の階層をズームレベル0とし、ズームレベルが1つ増えるごとに空間ボクセルの8分割を繰り返す階層構造を持つ
  • ズームレベル間で親子関係を持つ
  • 同一ズームレベルにおいて重複する空間ボクセルは存在しない

image.png
（4次元時空間情報利活用のための空間IDガイドライン（1.0版） 図2-3）

空間IDの適用範囲

上の図を見てなんとなく察した方もいるかもしれませんが、水平面の空間分割範囲はいわゆる「XYZタイル」と同一で、緯度±85.0511度外の範囲は切り捨てています。

鉛直方向を考慮しなければXYZタイルと全く同一の範囲になるため、既存のWebGISとは相性が良いです。

image.png
（4次元時空間情報利活用のための空間IDガイドライン（1.0版） 図2-10）

鉛直方向については±33,554,432mの範囲となっており、ズームレベルが上がるたびに2分（水平方向と掛け合わせて8分）されていきます。
この分割数に応じて指し示すインデックスが、{z}/{f}/{x}/{y}の中のfにあたりfインデックスと呼ばれます。

image.png
（4次元時空間情報利活用のための空間IDガイドライン（1.0版） 図2-7）

高さの基準は「標高0m」となっており、東京湾平均海面に一致します。
また、測地系はJGD2024またはWGS84となっています。

image.png
（4次元時空間情報利活用のための空間IDガイドライン（1.0版） 図2-5）

経緯度と標高・ズームレベルから空間IDを算出する方法はこちらに記載がありました。
https://github.com/unvt/zfxy-spec

Three-dimentional spatial identifier candidate ZFXY of a point with longitude lng [decimal degrees], latitude lat [decimal degrees] (lat_rad in radians), and elevation h [m] shall be encoded {z}/{f}/{x}/{y} when a zoom level z is given as an integer and f, x, y are given in the following formulae:

f = floor(n * h / H)
x = floor(n * ((lng + 180) / 360))
y = floor(n * (1 - log(tan(lat_rad) + (1 / cos(lat_rad))) / PI) / 2) where n = 2 ^ z, Z = 25 and H = 2 ^ Z [m].
Definitions of z, x, y are the same as ones in the industry standard Slippy Map Tilenames.

このほかにも時系列（4次元）を管理し、工事に利用してみたり…

image.png
（4次元時空間情報利活用のための空間IDガイドライン（1.0版） 図3-5）

ローカル空間IDを定義して、自由なIDを振る方法なんかも書いていました！

image.png
（4次元時空間情報利活用のための空間IDガイドライン（1.0版） 図2-19）

ざっくり概要がわかったところで、実装してみましょう！

経緯度・標高・ズームレベルから空間IDを算出する

すでに利用可能なJavaScriptの実装が存在するので、今回はこちらを参考にさせていただきました。

こんな感じのプロジェクト構成になりました。

❯ tree     
.
├── Cargo.toml
├── README.md
└── src
    ├── cell
    │   ├── mod.rs
    │   ├── spatial_cell.rs
    │   ├── xyz.rs
    │   └── zfxy.rs
    ├── coords
    │   ├── mod.rs
    │   └── types.rs
    ├── lib.rs
    └── projection
        ├── mod.rs
        └── web_mercator.rs

ソースコードはこちらにあります。
（調整中なので、あくまで参考程度にしてください）

重要なファイルはspatial-id-core/src/cell/spatial_cell.rsとspatial-id-core/src/cell/zfxy.rsです。

  • spatial_cell.rs
use std::fmt::Debug;
use std::hash::Hash;

pub type Vertices2D<T> = [T; 4];
pub type Vertices3D<T> = [T; 8];

pub trait SpatialCell: Debug + Clone + Copy + PartialEq + Eq + Hash {
    type Coord: Copy + Default + Debug;
    type Vertices: Copy + Default + Debug + AsRef<[Self::Coord]>;

    fn centroid(&self) -> Self::Coord;
    fn vertices(&self) -> Self::Vertices;
    fn bbox(&self) -> (Self::Coord, Self::Coord);
}

現状、2次元と3次元しかないですが、将来的には4次元のセルに拡張することも考えられます。
また、ローカルの空間IDは定義済みのVertices3Dをそのまま利用できると思っていますが、その場合には以下のようなリファクタリングが必要になりそうです。

  • 原点の位置を保存
  • 全体の各辺の長さを計算するロジック
  • 回転角

次はこのVertices3Dを実装したロジックです。

zfxy.rsでは色々書いていますが、メインはこのあたりになります。

  • zfxy.rs
const MAX_ALT: f64 = 33_554_432.0;

#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub struct ZFXY {
    pub zoom: u8,
    pub floor: i32,
    pub x: u32,
    pub y: u32,
}

impl SpatialCell for ZFXY {
    type Coord = LatLonAlt;
    type Vertices = Vertices3D<Self::Coord>;

    ...
    
    pub fn from_lat_lon_alt(lat: f64, lon: f64, alt: f64, zoom: u8) -> Self {
        let Some((x, y)) =
            crate::projection::web_mercator::latlon_to_tile_xy(lat, lon, zoom.into())
        else {
            panic!("Invalid lat/lon for tile conversion in from_lat_lon_alt");
        };
        let floor = alt.floor() as i32;
        ZFXY { zoom, floor, x, y }
    }

    pub fn to_lat_lon_alt(&self) -> LatLonAlt {
        let n = 2.0f64.powi(self.zoom as i32);
        let vz = MAX_ALT / n;
        let (lat, lon) =
            crate::projection::web_mercator::tile_xy_to_latlon(self.x, self.y, self.zoom);
        let alt = (self.floor as f64 + 0.5) * vz;
        LatLonAlt { lat, lon, alt }
    }

    pub fn to_spatial_id_str(&self) -> String {
        format!("/{}/{}/{}/{}", self.zoom, self.floor, self.x, self.y)
    }
    
    ...

}

実際の座標変換はspatial-id-core/src/projection/web_mercator.rsにおいています。

use std::f64::consts::PI;

pub const MAX_LATITUDE: f64 = 85.05112877980659;

pub fn latlon_to_tile_xy(lat: f64, lon: f64, zoom: u32) -> Option<(u32, u32)> {
    let lat_clamped = lat.clamp(-MAX_LATITUDE, MAX_LATITUDE);
    let lat_rad = lat_clamped.to_radians();
    let n = (2.0_f64).powi(zoom as i32);

    let tile_x_f = (lon + 180.0) / 360.0 * n;

    let tile_y_f = (1.0 - lat_rad.tan().asinh() / PI) / 2.0 * n;

    let tile_x = tile_x_f.floor().max(0.0) as u32;
    let tile_y = tile_y_f.floor().max(0.0) as u32;

    let max_tile_index = (1u32 << zoom).saturating_sub(1u32);
    Some((tile_x.min(max_tile_index), tile_y.min(max_tile_index)))
}

pub fn tile_xy_to_latlon(x: u32, y: u32, zoom: u8) -> (f64, f64) {
    let n = 2.0f64.powi(zoom as i32);
    let lon_deg = (x as f64 / n) * 360.0 - 180.0;

    let lat_rad = (PI * (1.0 - 2.0 * (y as f64) / n)).sinh().atan();
    let lat_deg = lat_rad.to_degrees();

    (lat_deg, lon_deg)
}

経緯度標高からインスタンスの生成することができ、生成されたZFXY構造体が空間IDを保持しまs

このファイルの中で以下のようなテストを書いています。

    #[test]
    fn test_from_lat_lon_alt() {
        let cell = ZFXY::from_lat_lon_alt(0.0, 0.0, 10.5, 25);
        assert_eq!(cell.zoom, 25);
        assert_eq!(cell.floor, 10);
        assert_eq!(cell.x, 16777216);
        assert_eq!(cell.y, 16777216);
    }

    #[test]
    fn test_to_lat_lon_alt() {
        let cell = ZFXY {
            zoom: 25,
            floor: 10,
            x: 16777216,
            y: 16777216,
        };
        let expected_lat_lon_alt = LatLonAlt {
            lat: 0.0,
            lon: 0.0,
            alt: 10.5,
        };
        assert_eq!(cell.to_lat_lon_alt(), expected_lat_lon_alt);

        let cell = ZFXY {
            zoom: 25,
            floor: 0,
            x: 16777216,
            y: 16777216,
        };
        let expected_lat_lon_alt = LatLonAlt {
            lat: 0.0,
            lon: 0.0,
            alt: 0.5,
        };
        assert_eq!(cell.to_lat_lon_alt(), expected_lat_lon_alt);

        let cell = ZFXY {
            zoom: 25,
            floor: 1,
            x: 16777216,
            y: 16777216,
        };
        let expected_lat_lon_alt = LatLonAlt {
            lat: 0.0,
            lon: 0.0,
            alt: 1.5,
        };
        assert_eq!(cell.to_lat_lon_alt(), expected_lat_lon_alt);

        let cell = ZFXY {
            zoom: 20,
            floor: 0,
            x: 524288,
            y: 524288,
        };
        let expected_lat_lon_alt = LatLonAlt {
            lat: 0.0,
            lon: 0.0,
            alt: 16.0,
        };
        assert_eq!(cell.to_lat_lon_alt(), expected_lat_lon_alt);
    }

こちらは実際に問題なく動作するんですが、先ほど紹介したJavaScriptのコードのテストを参考にしています。

  • 参考
import * as zfxy from "../src/zfxy";

describe('zfxy', () => {
  describe('getCenterLngLatAlt', () => {
    it('works', () => {
      const center1 = zfxy.getCenterLngLatAlt({z: 25, f: 0, x: 16777216, y: 16777216});
      expect(center1.alt).toStrictEqual(0.5);

      const center2 = zfxy.getCenterLngLatAlt({z: 25, f: 1, x: 16777216, y: 16777216});
      expect(center2.alt).toStrictEqual(1.5);

      const center3 = zfxy.getCenterLngLatAlt({z: 20, f: 0, x: 524288, y: 524288});
      expect(center3.alt).toStrictEqual(16);

      const center4 = zfxy.getCenterLngLatAlt({z: 20, f: 1, x: 524288, y: 524288});
      expect(center4.alt).toStrictEqual(32 + 16);

      const center5 = zfxy.getCenterLngLatAlt({z: 20, f: 10, x: 524288, y: 524288});
      expect(center5.alt).toStrictEqual((32 * 10) + 16);
    });
  });
});

cargo testなどのコマンドを実行するとテストが正常に実行されますね！

まとめ

今回は空間IDの概要とRustでの実装について説明してみました。

概念としてはほぼXYZタイルなのでWebGISに慣れている人であれば全く違和感なく理解できると思います。

鉛直方向の分割については、「標高0m」が基準のため日本独自路線感が否めませんが…日本での活用がちゃんと進めば（実装がたくさん出てきて、DB周りに取り込まれるなど）少なくとも日本国内ではデファクトになっていくんだろうと思います。

何度か説明したJavaScript実装の方ではタイルハッシュの出力や近傍のボクセルを見つける、親ズームレベルのボクセルを返す、などいろんな便利メソッドが実装されてそうだったので、そちらも今後実装できたら良さそうだなーと思っています！

3
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?