Ruby on GeoTools で距離計算３ - Relevant, Timely, and Accurate

geotools.rb を使って二点間の距離を求める記事の第3回です。第1回、第2回で、

地点名経度緯度

渋谷のハチ公 139.700638888889 35.6590111111111

つくば市のエキスポセンターの H-II ロケット 140.111305555556 36.0862305555556

http://d.hatena.ne.jp/hfu/20071002/1191299399

渋谷とつくばが含まれる座標系は、平面直角座標系 IX 系で、その EPSG コードは EPSG:2451

渋谷とつくばが含まれる座標系は、UTM 54N で、その EPSG コードは EPSG:32654

日本のいわゆる世界測地系は、JGD2000 で、その EPSG コードは EPSG:4612

http://d.hatena.ne.jp/hfu/20071003/1191473160

ということが分かりましたので、これを使って実際に距離を計算してみようと思います。

平面直角座標系 IX 系でのユークリッド距離

require 'geotools'

jgd2000 = Geo::import_epsg_crs(4612)
ix = Geo::import_epsg_crs(2451)

shibuya = Geo::import_array_geometry([139.700638888889, 35.6590111111111])
tsukuba = Geo::import_array_geometry([140.111305555556, 36.0862305555556])

tr = Geo::Transform.new(jgd2000, ix)
print "#{tr.transform(shibuya).distance(tr.transform(tsukuba))} m\n"

ユークリッド距離の計算に、JTS の関数 Geometry#distance を使っています。

$ ruby ix.rb 
DEBUG: rjb primitive_conversion mode
60180.7390166327 m

UTM 54N でのユークリッド距離

require 'geotools'

jgd2000 = Geo::import_epsg_crs(4612)
utm54n = Geo::import_epsg_crs(32654)

shibuya = Geo::import_array_geometry([139.700638888889, 35.6590111111111])
tsukuba = Geo::import_array_geometry([140.111305555556, 36.0862305555556])

tr = Geo::Transform.new(jgd2000, utm54n)
print "#{tr.transform(shibuya).distance(tr.transform(tsukuba))} m\n"

$ ruby utm.rb 
DEBUG: rjb primitive_conversion mode
60169.8952557551 m

大圏距離 1/2

まず、http://www.whitehat.net.nz/articles/2007/01/24/great-circle-distance-in-ruby にある、大圏距離算出の B アプローチのコードを使って計算:

# from: 
# http://www.whitehat.net.nz/articles/2007/01/24/great-circle-distance-in-ruby
def distance(lat1, lon1, lat2, lon2, radius = 6_372_795.0)
  include Math

  # Convert coordinates to radians
  lat1, lon1, lat2, lon2 = [lat1, lon1, lat2, lon2].map do |v|
      v * PI / 180.0
  end

  # Perform distance calculation
  distance = acos(
      cos(lat1) * cos(lat2) * cos(lon1 - lon2) + sin(lat1) * sin(lat2)
  ) * radius

  return distance
end

print "#{distance(139.700638888889, 35.6590111111111, 140.111305555556, 36.0862305555556)} m\n"

このコードは、地球を回転楕円体ではなくて球と近似して算出した大圏距離を求めているように見えます。

$ ruby gcd.rb 
58375.8226658971 m

平面直角座標系 IX 系、UTM　54N 上の「地図上の距離」とずいぶん違ってきています。

大圏距離 2/2

それから、GeoTools にも、GeodeticCalculator というクラスがあるので、これを使って計算してみます。このクラスは、楕円体を考慮して距離を計算してくれるようです。スレッドセーフでないクラスだそうです。

require 'geotools'

gc = Geo::Tools::GeodeticCalculator.new
gc.setStartingGeographicPoint(139.700638888889, 35.6590111111111)
gc.setDestinationGeographicPoint(140.111305555556, 36.0862305555556)
print gc.toString
print "#{gc.getOrthodromicDistance}m\n"

ここで、orthodromic distance という言葉が出てきますが、これは、

The orthodromic distance is the shortest distance between two points on a sphere's surface. The orthodromic path is always on a great circle. An other possible distance measurement is the loxodromic distance, which is a longer distance on a path with a constant direction on the compas.
http://deegree.sourceforge.net/deegree1.x.x_javadoc/org/deegree_impl/model/cs/Ellipsoid.html

ということですので、大圏距離と同じ意味であると思われます。

$ ruby geodetic.rb 
DEBUG: rjb primitive_conversion mode
Ellipsoid:     WGS84               
Source point:  139°42.0'E 35°39.5'N
Target point:  140°06.7'E 36°01.1'N
60186.64147204m

まとめ

渋谷のハチ公とつくばの H-II ロケットの距離は、以下のように計算されました。

距離の定義	距離	WGS84 大圏距離との比
平面直角座標系 IX 系での地図上の距離	60180.7390166327 m	0.99990193080619
UTM 54N 系での地図上の距離	60169.8952557551 m	0.999721761908036
地球を球と近似した大圏距離	58375.8226658971 m	0.969913277068565
地球を WGS84 で近似した大圏距離	60186.64147204 m	-

ほとんど同じ結果が出ています。特に、平面直角座標系については60km動いて5m台の違いしかなく、やはり投影法として優秀なものだと思います。

考察

大圏距離の計算には単純なユークリッド距離の計算と比べ、計算上のコストがかかります。距離の算出を頻繁に行う場合には、データを投影しておいてユークリッド距離の計算ですますということが有効かもしれません。

このエントリの信頼性について

このエントリの内容には誤りがあるかもしれません。

ChangeLog

Tue Oct 9 05:51:40 CEST 2007: 大圏距離の計算の入力の値に、誤った数字が混入していることに気がつきました。これを修正しました。

geotools.rb

このエントリで使用した geotools.rb は、http://svgmapdata.sakura.ne.jp/geotools/ に公開されています。

地点名	経度	緯度
渋谷のハチ公	139.700638888889	35.6590111111111
つくば市のエキスポセンターの H-II ロケット	140.111305555556	36.0862305555556