1点からの半径検索は、Google 投影に変換してから検索するのがベストプラクティス？ -- 第1回: Google 投影にふれる --

geotools.rb

1点からの半径検索は、データを単純なメルカトル投影 (例えばいわゆる Google 投影)に投影しておいて、検索点を検索距離だけ高緯度にずらした所の縮尺係数を掛けた矩形を使って検索をするのも手。

というような主張を、何回かのエントリに分けて行っていきたいと思っています。今回は、「いわゆる Google 投影」を導入します。

きっかけ

ここギコさんのエントリ：

latlongをUTM（ユニバーサル横メルカトル図法。m単位の直角座標系）に変更し、その上でdistance関数を使って範囲検索を行っていました。

http://kokogiko.net/m/archives/002135.html

問題があるとすれば、UTMは経緯度範囲で適用する座標系が異なるので、複数の座標系をまたがるような範囲に対してサービス提供するアプリケーションの場合、その範囲の座標系毎に関数インデックスを張っておかないといけないので、データ挿入・更新の多いアプリケーションだと遅くなるのかも。
また座標系を大きくまたがるような範囲の距離計算は不正確になるので、検索半径も見合った範囲に抑えておく必要がありそうです...まあ、経度6°を上回るような検索は誰もしないと思うけど。

http://kokogiko.net/m/archives/002135.html

に触発されて、

1点からの半径検索は、データを単純なメルカトル投影 (例えばいわゆる Google 投影)に投影しておいて、検索点を検索距離だけ高緯度にずらした所の縮尺係数を掛けた矩形を使って、空間インデクスを使った検索をするのも手。

という主張をしたいと思いました。

Google 投影の基本的性質を可視化してみる

そこで、geotools.rb で Google 投影を体験してみることにしました。

# geotools.rb を使う
require 'geotools'

# 経緯度座標系
wgs84 = Geo::import_epsg_crs(4326)
# 地球を球とするメルカトル投影。いわゆる Google 座標系
google = Geo::import_wkt_crs(<<-EOS
PROJCS["Google Mercator", 
  GEOGCS["WGS 84", 
    DATUM["World Geodetic System 1984", 
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], 
      AUTHORITY["EPSG","6326"]], 
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], 
    UNIT["degree", 0.017453292519943295], 
    AXIS["Geodetic latitude", NORTH], 
    AXIS["Geodetic longitude", EAST], 
    AUTHORITY["EPSG","4326"]], 
  PROJECTION["Mercator (1SP)", AUTHORITY["EPSG","9804"]], 
  PARAMETER["semi_major", 6378137.0], 
  PARAMETER["semi_minor", 6378137.0], 
  PARAMETER["latitude_of_origin", 0.0], 
  PARAMETER["central_meridian", 0.0], 
  PARAMETER["scale_factor", 1.0], 
  PARAMETER["false_easting", 0.0], 
  PARAMETER["false_northing", 0.0], 
  UNIT["m", 1.0], 
  AXIS["Easting", EAST], 
  AXIS["Northing", NORTH], 
  AUTHORITY["EPSG","900913"]]
EOS
)
# 経緯度座標系から Google 座標系への変換
tr = Geo::Transform.new(wgs84, google)
# Google 座標系から経緯度座標系への変換
itr = Geo::Transform.new(google, wgs84)

# 長さの歪みを検証する"基線"距離。単位はメートル
L = 250000
# Geotools 提供の測地計算器。大圏距離を計算してくれる。
gc = Geo::Tools::GeodeticCalculator.new
# JTS の幾何ファクトリ
gf = Geo::Tools::GeometryFactory.new
# 座標配列。JRuby 向け。
coords = com.vividsolutions.jts.geom.Coordinate[2].new

# "基線"の可視化のための Shapefile
distance = Geo::FeatureList.open('distances.shp')
# 縮尺係数を表っぽく保存するための Shapefile
scale_factor = Geo::FeatureList.open('scale_factor.shp')
# 計算をする経度についての繰り返し。メルカトル投影は経度に依存しないのでひとつで十分
135.step(135, 10) do |lng|
  # 計算をする緯度についての繰り返し。メルカトル投影は高緯度では発散。80度あたりで止める。
  -80.step(80, 10) do |lat|
    # 計算をする位置を経緯度座標系で
    g_wgs84 = Geo::import_wkt_geometry("POINT (#{lng} #{lat})")
    # 計算をする位置を Google 座標系で
    g_google = tr.transform(g_wgs84)
    # "基線" の原点は計算をする位置
    coords[0] = g_google.getCoordinate
    # "基線" の原点は計算をする位置
    gc.setStartingGeographicPoint(g_wgs84.getX, g_wgs84.getY)
    # scale_factor.shp に記憶するための地物ハッシュを作成
    f = {:the_geom => g_google}
    # 30 度おきの繰り返し。ここでの角度は通常の数学の通りに定義される
    0.step(2 * Math::PI, Math::PI / 6) do |rad|
      # Google 座標系における "基線" のもう一方の端 h_google
      h_google = g_google.clone
      # h_google は、"基線" の原点を 角度 rad の方向に L だけ translate したもの
      h_google.apply(Geo::Tools::AffineTransformation.translationInstance(
        L * Math.cos(rad), L * Math.sin(rad)))
      # h_google を経緯度座標系に変換した h_wgs84
      h_wgs84 = itr.transform(h_google)
      # "基線" g_wgs84 - h_wgs84 の間の大圏距離を計算させる
      gc.setDestinationGeographicPoint(h_wgs84.getX, h_wgs84.getY)
      # 縮尺係数「座標系上の距離 / 大圏距離」を地物ハッシュに格納する
      f["sf#{(rad / (2 * Math::PI) * 360).round}".to_sym] = 
        L / gc.getOrthodromicDistance
      # 今度は大圏距離で L 離れた点が投影座標上ではどこになるかを調べる。
      # 測地的な意味での方位角 (azimuth) と距離を測地計算器に与える。
      gc.setDirection(((rad - Math::PI) / (2 * Math::PI) * 360).round, L)
      # 大圏距離で測った "大圏距離基線" の端点 pt_wgs84
      pt_wgs84 = Geo::import_wkt_geometry("POINT (#{gc.getDestinationGeographicPoint.getX} #{gc.getDestinationGeographicPoint.getY})")
      # pt_wgs84 を Google 座標系に投影した pt_google
      pt_google = tr.transform(pt_wgs84)
      # pt_google を "大圏距離基線" の端点にセット
      coords[1] = pt_google.getCoordinate
      # "大圏距離基線" を distance.shp に書き込み
      distance.write({:the_geom => gf.createLineString(coords),
                     :sf => L / gc.getOrthodromicDistance})
    end
    # 縮尺係数表を scale_factor.shp に書き込み
    scale_factor.write(f)
  end
end
# distance.shp を閉じる
distance.close
# scale_factor.shp を閉じる
scale_factor.close

distances.shp の内容はこのような感じになります：

続きは、あとで書くことにします。