SQL Server

2015-06-17

緯度経度データ用にカスタムインデックスを設定する

今回は、緯度と経度から成る位置情報を含む地理データに、カスタムインデックスを設定する方法を説明します。カスタムインデックスを設定すると、データをリアルタイムにすばやく取得できます。たとえば、次のようなシナリオを考えてみます。以下の形式のデータが大量に (数百万項目) あるとします。

0001 47.610 -122.330
0002 36.175 -115.136
0003 47.613 -122.332
...

1 列目はユーザー ID (または位置と関連付けられる任意のオブジェクトの ID)、2 列目はユーザーが位置する緯度、3 列目は経度です。緯度の範囲は -90.0 ～ +90.0 度で、上下方向の座標を表します。経度の範囲は -180.0 ～ +180.0 度で、左右方向の座標を表します。1 人目と 3 人目のユーザーの位置はシアトル近辺です。赤道の緯度がゼロ、イギリスのグリニッジを通る本初子午線の経度がゼロです。データが SQL テーブルに格納されており、緯度と経度で 1 度ずつ区切ってメッシュを作成し、同じメッシュ内にいるユーザーをすべて見つけてユーザー 0001 と識別することにします。これは、次の簡単な SQL ステートメントで実現できます。

    SELECT userID
    FROM tblUsers
    WHERE latitude >= 47.0 AND latitude < 48.0
    AND longitude >= -123.0 AND longitude < -122.0

多くの場合、この方法でうまくいきます。ただし、リアルタイムのパフォーマンス (おそらく 3 秒以内の応答時間) が必要な場合、データのサイズが (システムのリソースに応じて決まる) しきい値を超えると、応答時間が著しく低下し始めます。

こうした状況でリアルタイムのパフォーマンスを実現する場合、各位置にセクター ID を割り当てる方法があります。この方法で、もともとのデータから次のような補助テーブルを生成できます。

ここでは、1 列目がユーザー ID、2 列目はセクター ID です。セクター ID にインデックスを設定すると、データセットに項目が数億個あっても、次の SQL ステートメントで非常に迅速に処理できます。

    SELECT userID
    FROM tblSectors
    WHERE sector = 49377

セクター ID は、主なデータセットのカスタムインデックスとして機能します。この場合、本当に難しいのは、緯度と経度を受け取ってセクターを返す関数を作成することです。

カスタムインデックスという考え方をよく理解するため、Windows Presentation Foundation (WPF) のサンプルアプリケーションを用意しました (図 2 参照)。このアプリケーションには、Bing マップコントロールが埋め込まれています。ワシントン州のレドモンドを地図の中心にして拡大したうえで、地図上でダブルクリックしました。ダブルクリックイベントは、クリックした地点の緯度と経度を割り出して赤い目印を表示するイベントハンドラーに関連付けられています。次に、アプリケーションはクリックした位置のセクター ID を算出したうえで、1 億項目のランダムのユーザー位置を格納するバックエンド SQL データベースを検索し、同じセクターにいるユーザーを 8 人見つけました。検索時間は 0.36 秒、非常に優れたパフォーマンスです。

図 1 カスタムインデックスの設計

メッシュの大きさ = 0.5	[-180.0, -179.5) 0	[-179.5, -170.0) 1	[-170.0, -169.5) 2		[179.5, 180.0] 719
[-90.0, -89.5) -> 0	0	1	2	...	719
[-89.5, -80.0) -> 1	720	721	722	...	1439
[-80.0, -79.5) -> 2	1440	1441	1442	...
	...

[89.5, 90.0] -> 359	258,480			...	259,199

Microsoft SQL Server 2008 には、このような場合に使用できる洗練された空間インデックス機能があります。さらに、多くの場合、SQL Server の空間インデックスを使用するのが最善の選択肢です。ただし、少なくとも 2 つのシナリオでは、カスタムインデックスを設定する方が良いでしょう。第 1 に、空間インデックスを作成できるのは、SQL 型のジオメトリか地理の列に限られます。レガシデータベースやレガシデータソースで緯度と経度が単なる数値として格納されている場合、その値を地理型に変換するのは現実的ではありません。第 2 に、SQL Server 2008 の空間インデックスは非常に強力でカスタマイズしやすいとはいえ、シナリオによってはインデックスの設定方法を完全に管理する必要が生じます。

カスタムインデックスの設定は新しい考え方ではありませんが、具体例はあまりありません。今回は、カスタムインデックスを設定して使用する方法の完全な例を紹介します。説明を完全に理解するには、少なくとも C# と T-SQL の中級レベルのコーディング能力が必要です。ここで紹介するサンプルの重要な C# コードは、archive.msdn.microsoft.com/mag201206Indexing (英語) から入手できます。

緯度と経度をセクター ID に対応付ける

緯度と経度の組み合わせをセクター ID に対応付けるにはさまざまな方法がありますが、ここでは 1 番簡単な方法を紹介します。これは、図 1 のダイアグラムに示す方法です。カスタムインデックスの設定方法を決めるには、まず地図をメッシュに分割する方法を選択します。図 1 の表の左上に「メッシュの大きさ = 0.5」と示しているように、それぞれの緯度 (行) と経度 (列) の 2 つの部分に分けます。その結果、緯度は [-90.0, -89.5) から [+89.5, +90.0] までの 360 の範囲に、経度は [-180.0, -179.5) から [+179.5, +180.0] までの 720 の範囲に分けられます。角かっこはその値を含むことを示し、丸かっこは含まないことを示します。このようにメッシュの大きさを 0.5 にすると、360 * 720 = 259,200 のセクターに分かれるため、0 ～ 259,199 の数を各メッシュに割り当てます。図 1 に示すように、セクターは左から右、上から下の順に並びます。

Custom Indexing Design
図 2 動作中のカスタム空間インデックス設定

メッシュの大きさを示すパラメーター (fraction) の値を小さくすると、より細かく分割することになるため、セクター数が増加します。この例では、緯度と経度を分割するのに同じ大きさを使用していますが、データの分布によっては異なる値を使用してもかまいません。ただし、すべてのセクターの面積は同じにはなりません。このことは後で説明します。行 (緯度) のインデックスを ri、列 (経度) のインデックスを ci とすると、セクター ID sid は以下のように求めます。

sid = (ri * 360 * 1/fraction) + ci

たとえば図 1 で、セクター 1441 は行インデックス 2 と列インデックス 1 に関連付けられていることから、(2 * 360 * 1/0.5) + 1 = (2 * 720) + 1 = 1441 と算出されています。360 * 1/fraction という項は、それぞれの行がいくつの列に分割されるかを計算し、これに ri を掛けて該当する行の先頭セクター ID を求めます。これに項 ci を追加することで、行の先頭から右へ ci 移動したセクター ID を求めています。

このパズルの最後は、緯度の値を行インデックスに、経度の値を列インデックスに対応付ける部分です。単にシンプルな線形検索を行う方法もありますが、大量のデータセットを扱った苦い経験から、この場合はバイナリサーチの方が適していることを学習しました。考え方としては、まず一番小さいインデックスをゼロとし、中央のインデックスと最も大きなインデックスを使用します。ここで、中央インデックスの緯度 (または経度) のメッシュを計算します。目的の緯度 (または経度) がそのメッシュに当てはまるれば、その中央インデックスが正しい戻り値です。目的の緯度 (または経度) が現在のメッシュよりも小さければ、最も小さいインデックスとそれまでの中央インデックスとの中央に、新しい中央インデックスを移動します。目的の緯度 (または経度) が現在のメッシュよりも大きければ、それまでの中央インデックスと最も大きなインデックスとの中央に、新しい中央インデックスを移動します。正しい中央インデックスが見つかるまで、このプロセスを繰り返します。

緯度の値とメッシュの大きさを受け取り一致する行インデックスを返す LatIndex メソッドの C# の実装を図 3 に示します。これと対になる経度の列インデックスを返す LonIndex メソッドは、定数 180 が 360 に、また定数 -90.0 と 90.0 が -180.0 と 180.0 に置き換わること以外はまったく同じです。

図 3 緯度の行/緯度インデックス

static int LatIndex(double latitude, double fraction)
{
  latitude = Math.Round(latitude, 8);
  int lastIndex = (int)(180 * (1.0 / fraction) - 1);
  double firstLat = -90.0;
  if (latitude == -90.0) return 0;
  if (latitude == 90.0) return lastIndex;
  int lo = 0;
  int hi = lastIndex;
  int mid = (lo + hi) / 2;
  int ct = 0; // To prevent infinite loop
  while (ct < 10000)
  {
    ++ct;
    double left = firstLat + (fraction * mid); // Left part interval
    left = Math.Round(left, 8);
    double right = left + fraction;         // Right part interval
    right = Math.Round(right, 8);
    if (latitude >= left && latitude < right)
      return mid;
    else if (latitude < left)
    {
      hi = mid - 1; mid = (lo + hi) / 2;
    }
    else
    {
      lo = mid + 1; mid = (lo + hi) / 2;
    }
  }
  throw new Exception("LatIndex no return for input = " + latitude);
}

図 3 のメソッドは緯度に double 型を使用しているため、入力緯度パラメーターを四捨五入して小数点第 8 位までに丸め、2 つの double 型を比較する場合の厄介なエラーを回避します。無限ループに陥るのを防ぐため、ct 変数をちょっとずるい方法で使用しています。コードを短く保つため、入力パラメーターの検証などの一般的なエラーチェックはほとんど省略しましたが、実際の運用シナリオでは省略しないでください。メインループはメッシュを [a,b) の形式でチェックするため、最後の 90.0 値を明示的にチェックする必要があることに注意します。

セクターの設計とヘルパーメソッドを配置し、次のようなメソッドで緯度経度の組み合わせをセクターに対応付けられます。

static int LatLonToSector(double latitude, double longitude,
  double fraction)
{
  int latIndex = LatIndex(latitude, fraction); // row
  int lonIndex = LonIndex(longitude, fraction);  // col
  return (latIndex * 360 * (int)(1.0 / fraction)) + lonIndex;
}

セクターサイズ

今回説明しているカスタムインデックスの設定方法では、メッシュの大きさのパラメーター値に基づいて地図を数値的に等間隔なメッシュに分割しています。たとえば、メッシュの大きさが 0.5 であれば、緯度と経度を分割する間隔は 0.5 度です。ただし、この場合、すべてのセクターの地理上の面積が同じなるわけではありません。図 4 を見てください。緯度の線は互いに平行に引かれており、同じ距離 (約 111 キロメートル) 離れています。そのため、図 4 のラベル A の距離は、ラベル B の距離と同じです。しかし、経度の線は極に近づくほど互いに近くなります。そのため、ラベル C の距離は、ラベル D の距離より短くなります。

Sectors Have Different Geographical Areas
図 4 セクターの地理上の面積の違い

地図上の任意の 2 点間のキロメートル単位の距離は、次のメソッドで測定できます。

static double Distance(double lat1, double lon1, 
  double lat2, double lon2)
{
  double r = 6371.0; // approx. radius of earth in km
  double lat1Radians = (lat1 * Math.PI) / 180.0;
  double lon1Radians = (lon1 * Math.PI) / 180.0;
  double lat2Radians = (lat2 * Math.PI) / 180.0;
  double lon2Radians = (lon2 * Math.PI) / 180.0;
  double d = r * Math.Acos((Math.Cos(lat1Radians) *
    Math.Cos(lat2Radians) *
    Math.Cos(lon2Radians - lon1Radians) +
   (Math.Sin(lat1Radians) * Math.Sin(lat2Radians)));
  return d;
}

そのため、あるセクターに対応する緯度と経度がわかれば、おおよその幅と高さを計算できるため、おおよその面積を計算できます。セクターを含むメッシュの左端の点を求めるには、次のメソッドを使用します。

static double SectorToLat(int sector,
    double fraction)
  {
    int divisor = 360 * (int)(1.0 / fraction);
    int row = sector / divisor;
    return -90.0 + (fraction * row);
  }

このメソッドは、本質的には LatLonToSector メソッドと逆のプロセスを実行します。たとえば、図 1 の入力セクターが 1441 でメッシュの大きさのパラメーター値が 0.5 の場合は、そのセクターのローカル divisor 変数は 360 * 1.0 / 0.5 = 720 で、経度のメッシュの数になります。そのため、変数 row の値は 1441 / 720 = 2 となり、これは行インデックスです (整数除算を行っていることに注意してください)。最後に、-90.0 + 0.5 * 2 = -90.0 + 1.0 = -89.0 と計算し、これがセクター 1441 に関連付けられた [-89.0, -79.5) のメッシュの左端になります。

セクターの経度のメッシュの左端を算出するメソッドもこれと似ていますが、まったく同じではありません。

static double SectorToLon(int sector, double fraction)
{
  int divisor = 360 * (int)(1.0 / fraction);
  int col = sector % divisor;
  return -180.0 + (fraction * col);
}

これら 2 つのヘルパーメソッドを使って、メッシュの大きさのパラメーターで定義したセクターのおおよその面積を平方キロメートル単位で算出できます。

static double Area(int sector, double fraction)
{
  double lat1 = SectorToLat(sector, fraction);
  double lon1 = SectorToLon(sector, fraction);
  double lat2 = lat1 + fraction;
  double lon2 = lon1 + fraction;
  double width = Distance(lat1, lon1, lat1, lon2);
  double height = Distance(lat1, lon1, lat2, lon1);
  return width * height;
}

隣接セクター

開発シナリオによっては、あるセクターに隣接するセクターを求めることがあります。たとえば図 1 で、ユーザーがセクター 721 にいる場合、隣接するセクター 0、1、2、720、722、1440、1441、および 1442 にいるユーザーを求めるとします。左右のセクターは、セクターが左端や右端でなければ、ID がプラス 1 またはマイナス 1 異なります。また、上下のセクターは、1 番上と 1 番下の行にあるセクターを除き、列のメッシュ数、つまり 360 * (1 / fraction) 分だけプラスまたはマイナスします。セクター (とセクターを生成するメッシュの大きさ) を受け取る AdjacentSectors メソッドを作成するには、セクターが地図の左端か右端、または 1 番上か 1 番下に位置するかどうかを把握すると役立ちます。そのための 4 つのメソッドを、図 5 に示します。

図 5 AdjacentSectors メソッド用のヘルパーメソッド

static bool IsLeftEdge(int sector, double fraction)
{
  int numColumns = (int)(1.0 / fraction) * 360;
  if (sector % numColumns == 0) return true;
  else return false;
}
static bool IsRightEdge(int sector, double fraction)
{
  if (IsLeftEdge(sector + 1, fraction) == true) return true;
  else return false;
}
static bool IsTopRow(int sector, double fraction)
{
  int numColumns = (int)(1.0 / fraction) * 360;
  if (sector >= 0 && sector <= numColumns - 1) return true;
  else return false;
}
static bool IsBottomRow(int sector, double fraction)
{
  int numColumns = (int)(1.0 / fraction) * 360;
  int numRows = (int)(1.0 / fraction) * 180;
  int firstValueInLastRow = numColumns * (numRows - 1);
  int lastValueInLastRow = numColumns * numRows - 1;
  if (sector >= firstValueInLastRow && sector <= lastValueInLastRow)
    return true;
  else
    return false;
}

AdjacentSectors メソッドの作成する方法はたくさんあるため、わかりやすさとコードサイズのバランスが重要です。1 つは、セクター値の配列を返す方法です。つまり、返す配列の [0] は入力セクターの左上に位置する隣接セクター、[1] は真上、[2] は右上、[3] は左、[4] は右、[5] は左下、[6] は真下、および [7] は右下のセクターとします。図 6 のコードは、AdjacentSectors を実装する 1 つの方法を大まかに示しています。入力セクターが地図の端ではない一般的な場合はわかりやすいですが、斜め方向の 4 つのセクターを含めて、いくつかの特殊な場合はわかりにくいでしょう。完全な実装を確認するには、この記事に付属のコードダウンロードを参照してください。

図 6 AdjacentSectors メソッドの実装

static int[] AdjacentSectors(int sector, double fraction)
{
  int[] result = new int[8]; // Clockwise from upper-left
  int numCols = (int)(1.0 / fraction) * 360;
  int numRows = (int)(1.0 / fraction) * 180;
  int firstValueInLastRow = numCols * (numRows - 1);
  int lastValueInLastRow = numCols * numRows - 1;
  // General case
  if (IsLeftEdge(sector, fraction) == false &&
    IsRightEdge(sector, fraction) == false &&
    IsTopRow(sector, fraction) == false &&
    IsBottomRow(sector, fraction) == false)
  {
    result[0] = sector - numCols - 1;
    result[1] = sector - numCols;
    result[2] = sector - numCols + 1; 
    result[3] = sector - 1;
    result[4] = sector + 1;
    result[5] = sector + numCols - 1;
    result[6] = sector + numCols;
    result[7] = sector + numCols + 1;
    return result;
  }
  // Deal with special cases here. See code download.
  throw new Exception("Unexpected logic path in AdjacentSectors");
}

地理的インデックス設定の構想を完成させる

今回説明したカスタム地理インデックスの設定方法を理解する優れた方法は、簡単な例を最初から最後まで調べることです。まず、次のような T-SQL コードを使用して、ダミーの SQL データベースを作成します。

    use master
    if exists(select * from sys.sysdatabases where name='dbGeoData')
      drop database dbGeoData
    create database dbGeoData on
    (
    name=dbGeoData,
    filename='D:\SomePlace\dbGeoData.mdf'
    )

次に、ダミーデータを格納するテーブルを作成します。

    use dbGeoData
    create table tblUsers
    (
    userID int primary key,
    latitude real not null,
    longitude real not null
    )

続いて、次のような C# コードを作成してダミーデータを生成します。

Random r = new Random(0);
string initialDataFile = "..\\..\\UserIDLatLon.txt";
FileStream ofs = new FileStream(initialDataFile, FileMode.Create);
StreamWriter sw = new StreamWriter(ofs);
for (int i = 0; i < 1000000; ++i)
{
  double lat = (90.0 - (-90.0)) * r.NextDouble() + (-90.0);
  double lon = (180.0 - (-180.0)) * r.NextDouble() + (-180.0);
  sw.WriteLine(i.ToString().PadLeft(6,'0') + "," +
    lat.ToString("F8") + "," + lon.ToString("F8"));
}
sw.Close(); ofs.Close();

ここでは、コンマ区切りのデータを百万行作成します。範囲 [hi,lo) の緯度と経度をランダムに生成するため、(hi - lo) * random.NextDouble() + lo という標準パターンを使用しています。ここで、コマンドラインで bcp.exe プログラムを使用して、ダミーデータを SQL テーブルにコピーします。

> bcp.exe dbGeoData..tblUsers in UserIDLatLon.txt -c -t , -r \n -S(local) -T

このコマンドは、「UserIDLatLon.txt ファイルのデータを、dbGeoData データベースにある tblUsers テーブルにコピーします。これは、各列がコンマによって区切られており、各行が改行記号で終わっている文字データ (つまりテキストファイル) です。データベースサーバーはローカルのコンピューターにあり、信頼関係 (Windows) 認証を使用して接続します」という意味です。

さらに、次のような C# コードで、補助のカスタムインデックスを作成するセクターデータを作成します。

string initialDataFile = "..\\..\\UserIDLatLon.txt";
string sectorDataFile = "..\\..\\ UserIDSector.txt";
FileStream ifs = new FileStream(initialDataFile, FileMode.Open);
StreamReader sr = new StreamReader(ifs);
FileStream ofs = new FileStream(sectorDataFile, FileMode.Create);
StreamWriter sw = new StreamWriter(ofs);
string line = "";
string[] tokens = null;
while ((line = sr.ReadLine()) != null)
{
  tokens = line.Split(',');
  int userID = int.Parse(tokens[0]);
  double lat = double.Parse(tokens[1]);
  double lon = double.Parse(tokens[2]);
  int sector = LatLonToSector(lat, lon, 0.5);
  sw.WriteLine(userID.ToString().PadLeft(6, '0') + "," + sector);
}
sw.Close(); ofs.Close(); sr.Close(); ifs.Close();

コードはダミーデータファイルを 1 行ずつループし、ユーザー ID、緯度、および経度を解析し、メッシュの大きさのパラメーターを 0.5 として、関連付けられるセクターを算出し、ユーザー ID とセクターをコンマ区切り形式でテキストファイルに書き込みます。

次に、セクターデータを格納する SQL テーブルを作成します。

    create table tblSectors
    (
    userID int primary key,
    sector int
    )

続いて、セクターデータをテーブルにコピーします。

> bcp.exe dbGeoData..tblSectors in UserIDSector.txt -c -t , -r \n -S(local) -T

最後に、セクターデータにインデックスを付けます。

    if exists(select name from sys.sysindexes where name='ndxSector')
      drop index ndxSector on tblSectors
    create nonclustered index ndxSector on tblSectors(sector)

これは、複数の方法で試せます。たとえば、次のような T-SQL コードで、同じセクターにいるすべてのユーザーをユーザー 3 として選択できます。

declare @userID int
set @userID = 3
declare @sector int
select @sector = sector from tblSectors where userID=@userID
print @sector
select userID from tblSectors where sector = @sector

当然、アプリケーションから SQL データにアクセスしている場合は、この種の SQL コードに相当する ADO.NET や LINQ を使用できます。

リアルタイムのパフォーマンスを実現する

カスタムインデックスの設定は、多くのアプリケーションのシナリオで使用できます。たとえば、Bing マップコントロールを備えた ASP.NET Web アプリケーションまたは Silverlight アプリケーションがあるとします。ここで、ユーザーが地図をクリックできるようにするシステムを作成しました。クリックイベントは、クリックした緯度と経度を計算し、対応するセクターを算出したうえ、バックエンドの SQL データベースからセクター内のすべてのユーザーを取得するコードと関連付けられています。図 1 に示すような行が数億行ある場合でも、リアルタイムのパフォーマンスを実現できます。

複数のカスタムインデックスを設定することがあります。たとえば、メッシュの大きさ = 1.0 (緯度経度 1 度ずつの範囲を 1 セクターとする) という粒度の低いインデックス、メッシュの大きさ = 0.25 の粒度が中程度のインデックス、およびメッシュの大きさ = 0.05 の粒度の高いインデックスが必要な場合を考えてみます。複数のインデックスを使用すれば、SQL データを段階的にフィルター処理できます。まず、粒度の低いセクターに何人ユーザーがいるかを算出します。取得したユーザー数が多すぎる場合は、次に中程度の粒度のセクターにいるユーザー数を算出します。ユーザーが少なすぎれば、7 つの隣接セクターとそれらのセクター内のユーザー数を算出する、といった具合です。

繰り返しになりますが、SQL Server 2008 には強力な空間インデックス機能が組み込まれており、カスタムインデックスを設定する前に使用すべきかどうかを検討する必要があります。カスタムインデックスは、SQL テーブルを明示的に追加する必要があり、システム管理の負担が増します。ただし、シナリオによっては、カスタムインデックスを設定するとパフォーマンスが非常に速くなります。GPS 機能の付いたモバイルデバイスの使用が増え、大量の地理的データを収集する機能は拡大するばかりです。カスタムインデックスの設定を設計する能力は、こうしたデータの分析に非常に役立ちます。

Dr. James McCaffrey は Volt Information Sciences に勤務し、ワシントン州レドモンドにあるマイクロソフト本社で働くソフトウェアエンジニアの技術トレーニングを管理しています。これまでに、Internet Explorer、MSN サーチなどの複数のマイクロソフト製品にも携わってきました。また、『.NET Test Automation Recipes』(Apress、2006 年) の著者でもあります。連絡先は、jammc@microsoft.com (英語のみ) です。

この記事のレビューに協力してくれた技術スタッフの Isaac Kunen と Dan Liebling に心より感謝いたします。

次の方法で共有