RDBMS-GIS界隈でははじめまして（あるいはただいま？）、@kochizufanと申します。
20年ほど前はRDBMS-GIS界隈の第一人者だった¹のですが、ここ6～7年はサーバサイドを使わないフロントエンドだけで解決する技術を追求していたこともあり、RDBMS…？何それお高い²んでしょう…？と距離を置いていました。
そのため完全に浦島太郎状態ですが、リハビリかねてちょっといろいろ使ってみました。
どうぞよろしくお願いいたします。

今回の成果物

今回の成果はこちらのgithubに公開しています。

MySQLのJSON対応機能

WebGISエンジニア、特にフロントエンドしか使わねえぜと心に決めてきたものにとって、一番の友達はGeoJSON仕様³なので、MySQLに最近ついた機能でやっぱり興味津々なのはJSON対応機能です。

静的Webストレージに置いておくだけで使えるのと、あとはテーブル型ではない構造的データも扱えるので重宝したGeoJSONですが、それがテーブル型データベースで扱えるとは面白い、ということで早速使ってみました。

Node.JSからのmysqlへの接続

まずは接続。
フロントエンド開発者なので基本は開発に使うのはNode.JSです。
Node.JSからmysqlに繋げるライブラリを探していたところ、今はmysql2というのを使うようです⁴。

あと、async/awaitを使いたければ、promiseインタフェースに対応する必要があるので、呼び出すライブラリはmysql2/promiseになります。
まとめると接続はこんな感じ。

import mysql from "mysql2/promise";
import dotenv from "dotenv";
const res = dotenv.config();
const env = res.parsed;

const main = async () => {
  // Connection
  const connection = await mysql.createConnection({
    host: 'localhost',
    user: env.USER,
    password: env.PASS
  });

  // Drop Database
  await connection.query("DROP DATABASE geojson;");

  // Create & Use Database
  await connection.query("CREATE DATABASE geojson;");
  await connection.query("USE geojson;");

まとめてgeojsonを叩き込んでみた...がそれだと使いようがない

続いて、さっそくJSONカラムを持ったテーブルを作って、GeoJSONを入れてみました。
サンプルに使ったデータは、私の活動である館林の石造物調査データを用いました。

まずは全部読み込んで、叩き込んでみましたが...。

import fs from "fs-extra";
const geojson = fs.readJsonSync("./tatebayashi_stones.geojson");

export async function whole_geojson(connection) {
  console.log("## Result of geojson case (Whole geojson)");
  // Create Table
  // For insert whole GeoJSON
  await connection.query(`CREATE TABLE whole_geojson(
    id INT NOT NULL AUTO_INCREMENT,
    geojson JSON,
    PRIMARY KEY (id)
  );`);

  // Insert whole GeoJson
  console.log(`Total number of features: ${geojson.features.length}`);
  console.log("Installing...");
  await connection.query('INSERT INTO whole_geojson (geojson) VALUES (?);', JSON.stringify(geojson));
  console.log("Installed");

  console.log("No idea...")
}

GeoJSONの全件データが入った行が1つだけできましたが、これ、これ以上どうしようもないですね。
検索しようにも検索条件が設定できないし（DBとしては行1つしかないから検索もクソもない）、なるほどそういう使い方ではないと。
最低限、FeatureCollectionから各featureに分けるところはアプリ側で行って、DBに流し込まないといけないようです。

FeatureCollectionを分解してDBに登録、検索

行を分解して登録し、検索まで行ってみます。

import fs from "fs-extra";
const geojson = fs.readJsonSync("./tatebayashi_stones.geojson");

export async function each_row_geojson(connection) {
  console.log("## Result of geojson case (Each rows)");
  // Create Table
  // For insert GeoJSON records
  await connection.query(`CREATE TABLE each_row_geojson(
    id INT NOT NULL,
    geojson JSON,
    latlong GEOMETRY AS (ST_GeomFromGeoJSON(geojson, 1, 4326)) STORED NOT NULL SRID 4326,
    PRIMARY KEY (id),
    SPATIAL INDEX sp_index1 (latlong)
  );`);

  // Insert GeoJson's each record
  console.log(`Total number of features: ${geojson.features.length}`);
  console.log("Installing...");
  for (let i = 0; i < geojson.features.length; i++) {
    const record = geojson.features[i];
    await connection.query('INSERT INTO each_row_geojson (id, geojson) VALUES (?, ?)', [record.properties.fid, JSON.stringify(record)]);
  }
  console.log("Installed");

  // Select by MBR: (139.52246270948356, 36.251589025204765) To (139.54034174583242, 36.24221711199563)
  console.log("Select by MBR: (139.52246270948356, 36.251589025204765) To (139.54034174583242, 36.24221711199563)");
  await connection.query('SET @poly1 = ST_SRID(ST_MakeEnvelope(Point(?, ?), Point(?, ?)), 4326)', [139.52246270948356, 36.251589025204765, 139.54034174583242, 36.24221711199563]);
  const [buf] = await connection.query('SELECT ST_AsGeoJSON(@poly1) AS geom');
  // Save MBR as GeoJSON
  fs.writeJsonSync("./mbrGeojson.geojson", {type:"Feature", geometry: buf[0]["geom"]});
  const [res1] = await connection.query('SELECT geojson FROM each_row_geojson WHERE ST_Within(latlong, @poly1)');
  console.log(`Number of selected features: ${res1.length}`);
  // Save result as GeoJSON
  geojson.features = res1.map((row) => { return row.geojson });
  fs.writeJsonSync("./selectedGeojson.geojson", geojson);

  // Select by MBR with index
  console.log("EXPLAIN ANALYZE: With index");
  const [res2] = await connection.query('EXPLAIN ANALYZE SELECT geojson FROM each_row_geojson WHERE ST_Within(latlong, @poly1)');
  console.log(`EXPLAIN: ${res2[0]['EXPLAIN']}`);

  // Select by MBR without index
  // Invisible Index
  await connection.query('ALTER TABLE each_row_geojson ALTER INDEX sp_index1 INVISIBLE');
  console.log("EXPLAIN ANALYZE: Without index");
  const [res3] = await connection.query('EXPLAIN ANALYZE SELECT geojson FROM each_row_geojson WHERE ST_Within(latlong, @poly1)');
  console.log(`EXPLAIN: ${res3[0]['EXPLAIN']}`);
}

これのテーブル生成SQLだけ取り出してみると、

CREATE TABLE each_row_geojson(
  id INT NOT NULL,
  geojson JSON,
  latlong GEOMETRY AS (ST_GeomFromGeoJSON(geojson, 1, 4326)) STORED NOT NULL SRID 4326,
  PRIMARY KEY (id),
  SPATIAL INDEX sp_index1 (latlong)
);

となりますが、ここでポイントは3つくらいありそうです。

geojsonカラムに含まれる位置情報から検索インデックスを作るには、Generated Columnを使う
既存カラムから検索可能なインデックスを作るには、Generated Columnを使う方法と、関数インデックスを使う方法が考えられますが、今回はGenerated Columnを使う方法で試しました。
試してないのでわかりませんが、空間検索は関数出力に演算子で判定するような形ではなく、ST_Withinだのの関数の中で利用されるので、関数インデックスアプローチだと動作しない気がします。
Generated ColumnのタイプはSTOREDを選ぶ
Generated Columnは生成したカラムをメモリ内にしか持たないか、ストレージにも蓄積するかでVIRTUALとSTOREDという設定があるようですが、STOREDでないと空間インデックスが張れないようです。ストレージ容量は食いますがインデックスなしとかDBに入れる意味がないので、STOREDにしましょう。
SRIDを与える
SRIDというのは、GISで扱う座標系を見分ける番号です。
専門でない人には意外かもしれませんが、位置を表すのに経緯度というのは一意な尺度ではなく、世界には様々な経緯度の尺度があります。
GNSS⁵などない時代から人類は経緯度を使ってきているので、地球というでこぼこの天体を近似するのにどういう回転楕円体を選ぶかというところから始まって、様々な定義が世界中に存在します。
さらに経緯度だけではなく、それをメルカトル図法だのユニバーサル横メルカトル図法だの、平面直角座標系だの様々な地図図法に変換した後の座標系も存在します。
そのような、どの経緯度定義を使うか、さらに地図図法として何を使うか、というセットに識別IDを与えたものがSRIDです。
このSRIDを指定しないと、異なるSRIDの座標の間で様々な幾何演算を行っても、それは電圧計算で0ボルト電位を合わせていないようなものなので、意味のある結果は出ないので、利用するSRIDを指定する必要があります⁶。
今回指定しているSRID 4326は、GNSSで採用されて、世界標準となっているWGS84測地系を、地図座標返還せず経緯度のまま用いる場合に指定するSRIDです。

この点に気を付けてテーブルを作成し、GeoJSONのFeatureCollectionを分解してテーブルに流し込み、範囲で検索すると、ちゃんと結果が返ってきます。
示している図は結果をQGISで表示したもので、緑の点が全体のGeoJSONデータ、赤い枠が検索範囲、赤の点が検索結果です。
ちゃんと選択結果がGeoJSONとしての属性も維持して結果が返ってきているのがわかります。

EXPLAIN ANALYZEの結果も、インデックスを使わない場合との比較含め出力してみました。
ちゃんとINDEXが使われているのがわかります。

Select by MBR: (139.52246270948356, 36.251589025204765) To (139.54034174583242, 36.24221711199563)
Number of selected features: 151
EXPLAIN ANALYZE: With index
EXPLAIN: -> Filter: st_within(each_row_geojson.latlong,<cache>((@poly1)))  (cost=37.61 rows=83) (actual time=0.277..8.141 rows=151 loops=1)
    -> Index range scan on each_row_geojson using sp_index1 over (latlong unprintable_geometry_value)  (cost=37.61 rows=83) (actual time=0.196..1.428 rows=151 loops=1)

EXPLAIN ANALYZE: Without index
EXPLAIN: -> Filter: st_within(each_row_geojson.latlong,<cache>((@poly1)))  (cost=192.25 rows=1360) (actual time=0.294..102.266 rows=151 loops=1)
    -> Table scan on each_row_geojson  (cost=192.25 rows=1360) (actual time=0.060..5.924 rows=1997 loops=1)

FlatGeoBufを使った例

面白かったので、GeoJSONだけでなく、FlatGeoBufを使った例も試してみました。
FlatGeoBufとは、GeoJSONと相互互換でバイナリ化することによりデータサイズをはるかに小さくできると同時に、HTTPのByte Rangeリクエストでバイトレンジ部分リクエストすることにより、単独ファイルでタイル配信的機能を実現できるクラウドフレンドリーなすごい仕様です⁷。

FlatGeoBufのサイトを見ると、MySQLのライバル？のオープンソースGIS RDBMS、PostGISはこのFlatGeoBufに対応しているようなので、MySQLも将来対応してほしいですね⁸。

とりあえずサンプルデータとしては、FlatGeoBufのLeafletサンプルで使われていた、USの郡ポリゴン情報を使ってみました。

import { geojson as fgb } from 'flatgeobuf';
import { readFileSync }  from 'node:fs';
import fs from "fs-extra";
const buffer = readFileSync('./UScounties.fgb');
const geojson = fgb.deserialize(new Uint8Array(buffer));

export async function us_flatgeobuf(connection) {
  console.log("## Result of flatgeobuf case");
  // Create Table
  // For insert FlatGeoBuf
  await connection.query(`CREATE TABLE us_flatgeobuf(
    id INT NOT NULL AUTO_INCREMENT,
    geojson JSON,
    latlong GEOMETRY AS (ST_GeomFromGeoJSON(geojson, 1, 4326)) STORED NOT NULL SRID 4326,
    PRIMARY KEY (id),
    SPATIAL INDEX sp_index2 (latlong)
  );`);

  // Insert FlatGeoBuf's each record
  console.log(`Total number of features: ${geojson.features.length}`);
  console.log("Installing...");
  for (let i = 0; i < geojson.features.length; i++) {
    const record = geojson.features[i];
    await connection.query('INSERT INTO us_flatgeobuf (geojson) VALUES (?)', [JSON.stringify(record)]);
  }
  console.log("Installed");

  // Select by MBR: (-124.42278337796137, 39.3058971532753) To (-114.7085706097631, 35.73909681729993)
  console.log("Select by MBR: (-124.42278337796137, 39.3058971532753) To (-114.7085706097631, 35.73909681729993)");
  await connection.query('SET @poly2 = ST_SRID(ST_MakeEnvelope(Point(?, ?), Point(?, ?)), 4326)', [-124.42278337796137, 39.3058971532753, -114.7085706097631, 35.73909681729993]);
  const [buf] = await connection.query('SELECT ST_AsGeoJSON(@poly2) AS geom');
  // Save MBR as GeoJSON
  fs.writeJsonSync("./mbrFlatgeobuf.geojson", {type:"Feature", geometry: buf[0]["geom"]});
  const [res1] = await connection.query('SELECT geojson FROM us_flatgeobuf WHERE ST_Within(latlong, @poly2)');
  console.log(`Number of selected features: ${res1.length}`);
  // Save result as GeoJSON
  geojson.features = res1.map((row) => { return row.geojson });
  fs.writeJsonSync("./selectedFlatgeobuf.geojson", geojson);

  // Select by MBR with index
  console.log("EXPLAIN ANALYZE: With index");
  const [res2] = await connection.query('EXPLAIN ANALYZE SELECT geojson FROM us_flatgeobuf WHERE ST_Within(latlong, @poly2)');
  console.log(`EXPLAIN: ${res2[0]['EXPLAIN']}`);

  // Select by MBR without index
  // Invisible Index
  await connection.query('ALTER TABLE us_flatgeobuf ALTER INDEX sp_index2 INVISIBLE');
  console.log("EXPLAIN ANALYZE: Without index");
  const [res3] = await connection.query('EXPLAIN ANALYZE SELECT geojson FROM us_flatgeobuf WHERE ST_Within(latlong, @poly2)');
  console.log(`EXPLAIN: ${res3[0]['EXPLAIN']}`);
}

QGISで見た結果はこんな感じ。

EXPLAIN ANALYZEの結果は以下の通り。

Select by MBR: (-124.42278337796137, 39.3058971532753) To (-114.7085706097631, 35.73909681729993)
Number of selected features: 32
EXPLAIN ANALYZE: With index
EXPLAIN: -> Filter: st_within(us_flatgeobuf.latlong,<cache>((@poly2)))  (cost=28.57 rows=32) (actual time=0.410..78.698 rows=32 loops=1)
    -> Index range scan on us_flatgeobuf using sp_index2 over (latlong unprintable_geometry_value)  (cost=28.57 rows=32) (actual time=0.205..28.761 rows=36 loops=1)

EXPLAIN ANALYZE: Without index
EXPLAIN: -> Filter: st_within(us_flatgeobuf.latlong,<cache>((@poly2)))  (cost=1859.35 rows=1608) (actual time=262.368..2130.596 rows=32 loops=1)
    -> Table scan on us_flatgeobuf  (cost=1859.35 rows=1608) (actual time=0.046..55.046 rows=3221 loops=1)

MySQLの空間実装：回転楕円体面 or 球面として計算されている？

GeoJSONやFlatGeoBufをMySQLのJSON処理機能に流し込んで、GIS機能とも連携してみた実験については、前章までの内容になります。
最後に、今回試してみた中で見つけた面白い結果を紹介します。

今回、USの郡の検索で、SRID4326座標系（世界測地系）で格納された地理情報に対し、矩形でST_Withinでの検索を試みたのですが、その結果をQGISで表示すると面白いことがわかりました。
面白い結果はここに表れています。

ST_Withinは「内部に含む」結果を計算するはずなのに、QGISで見ると、この選択された郡は、検索元の矩形（上の赤線）をはみ出ています。
これ、バグでなく正しい結果が出ているとして考えると、おそらくMySQLの空間演算が、平面演算ではなくちゃんと回転楕円体面演算、あるいは球体近似演算されていることが原因じゃないかと思います。
検索対象の矩形をQGIS上で描画すると、今回の描画は球面メルカトル座標（SRID:3857）を用いているため、同じ緯度の2点を結ぶと、その間の線は緯度の線上に載ります。
しかし、実際の最短距離を結ぶ線（大圏コース）は、北半球では緯度の線上よりだいぶん北を通ります。
そのため、検索結果に含まれたのではないかと推察されます。

私のように古いGIS RDBMSユーザの思い込みからだと、各SRIDの指定の差は座標返還処理の基準としてのみ使われて、空間演算は各座標系をデカルト座標系として計算されるものと思い込んでいましたが、MySQLのリファレンスの5.6から8.0への記述の違いを見ると

ただし、すべての計算が実際の SRID 値には関係なく、デカルト (平面) 座標を表す SRID 0 を想定して実行されます。

SRID 0 は、軸に単位が割り当てられていない無限平坦なデカルト平面を表します。 SRID 0 の動作を保証するには、SRID 0 を使用してジオメトリ値を作成します。 SRID 0 は、SRID が指定されていない場合の新しいジオメトリ値のデフォルトです。

となっているので、SRID:0以外は回転楕円体面(or 球面)として計算するようになったのかもしれません。

多くのユースケースに合うように改善されたともいえると思いますが、一方で、

たとえば南北回帰線内に存在する地物を検索するなどのように、意図的に球面検索したくない場合はどうなるのか？
地球の上のデータだけでなく、例えば月の上の地物を管理して空間演算などをする場合、計算できるのか？⁹

などが疑問に感じるので、また詳細を調査、実験してみたうえで、結果を報告してみたいと思います。

技術的第一人者という意味ではなく、情報発信してたので検索などすれば必ず行き当たるという意味での第一人者です。技術的には昔から変わらず、全方面圧倒的ど素人です。戦線は当時よりえらい拡大しましたが。 ↩
金銭的コストだけでなく、運用コスト含め、余暇での開発ではフロントエンドとバックエンドの双方を網羅したフルスタックの技術力を維持するのは困難と判断したので、フロントエンドのみの技術を維持して開発することに専念していました。 ↩
ただし、最近はFlatGeoBuf仕様が後継として現れてきています。後述。 ↩
mysql というライブラリもありますが、こちらはmysqlの最新のデフォルト認証プラグインであるcaching_sha2_passwordに対応していないので、今は使えないようです。 ↩
米国のGPSだけでなく、EUのGalileo、ロシアのGLONASS、中国の北斗など、全世界のGPSライクなシステム全体を包含する一般語彙。 ↩
というか、MySQL8.0以降では、SRIDが指定されないと空間インデックスが効かない仕様になっているようですね。 ↩
個人的にサーバサイドの不要なフロントエンドだけで成立する技術を追求してきた関係上、2018年頃から、このFlatGeoBufのように単独ファイルとして扱えつつHTTR Byte Rangeリクエストでクラウド対応できる仕様が欲しい（というか作れるなら作りたい）と言い続けてきたのですが、ここ1～2年でベクトルデータはFlatGeoBuf、ラスタデータはCloud Optimized GeoTiffとまさにやりたかった事が実現できる仕様が揃ってきたのでうれしい限りです。 ↩
もっとも、FlatGeoBufに対応して嬉しいのは大容量のFeatureCollectionのようなケースだと思うので、個々のFeature単位で管理するRDBMSでFlatGeoBufを使えてどのくらい嬉しいかは疑問符も点灯しますが...。ユースケースを考えた実装が求められそうです。 ↩
これについては、計算の基準となる楕円体を、与えられた楕円体を用いて計算していたら問題なし、WGS84の回転楕円体に変換したうえで計算していたら対応不可と思われます。しかし与えられた楕円体を用いて空間計算した場合は、今度は基準楕円体が違う空間参照系を用いたデータセットでは、空間計算の結果がずれるという問題が発生しそうな気がします。 ↩