Ruby on GeoTools で距離計算 2

Ruby on GeoTools geotools.rb を使って、2点間の直線距離を計算してみる第2回です。

地点名 経度 緯度
渋谷のハチ公 139.700638888889 35.6590111111111
つくば市のエキスポセンターの H-II ロケット 140.111305555556 36.0862305555556
http://d.hatena.ne.jp/hfu/20071002/1191299399

ということが分かりましたので、geotools.rb を使って距離を計算したいと思います。

計算する距離

を算出して、比べてみたいと思います。

渋谷とつくばは平面直角座標系の何系か

http://vldb.gsi.go.jp/sokuchi/patchjgd/download/Help/jpc/jpc.htmによると、東京都の島しょ部を除く全ての市町村と茨城県の全市町村は、平面直角座標系の IX 系に所属することがわかります。
この座標系の EPSG コードを調べることにより、平面直角座標系 IX 系という座標系を GeoTools に分からせることができます。例えば、「平面直角座標系 IX EPSG」を Google で検索することにより、例えば id:hogeman さんのエントリから平面直角座標系 IX 系の EPSG コードは 2451 であることが分かります。
geotools.rb に、座標系の名前から EPSG コードを検索する機能を、あとで追加するかもしれません。*1

系をまたぐ場合には投影後の距離計算は難しい

逆に考えると、系をまたぐ場合、平面直角座標系に投影したあとで距離を計算するのは難しいことになります。
大圏距離にはこのような制約はありませんので、大圏距離が平面直角座標系投影後の距離とあまり違わないようであれば、大圏距離を使った方が余計な心配をしなくていいことになります。

渋谷とつくばは UTM 座標系の何系(ゾーン)か

UTM (Universal Transverse Marcator) という座標系のオリジナルの仕様は、とくにオーソライズされたものはないようです。UTM は、もともと NATOワルシャワ条約機構が統一規格として普及をはじめたもののようですが、現在では世界中で使われていて、日本の 1:25,000 地形図にも、UTM 座標系が使われています。
UTM に関する情報としては、http://www.uwgb.edu/DutchS/FieldMethods/UTMSystem.htm にあるものが良い情報であるように見えます。
ここによると、UTM のゾーンを求めるには、

To find the grid zone for any longitude:
Treat west longitude as negative and east as positive.
Add 180 degrees; this converts the longitude to a number between zero and 360 degrees.
Divide by 6 and round up to the next higher number.

http://www.uwgb.edu/DutchS/FieldMethods/UTMSystem.htm

とせよとあります。 zone という表現だとか、round ではなくて round up (ゼロを嫌う) させるところとかが、mil 規格っぽい感じがしないでもありません。
実際に計算してみると、

$ irb
irb(main):001:0> [139.700638888889, 140.111305555556].map {|lng| ((lng + 180) / 6).ceil}
=> [54, 54]

ですので、渋谷もつくばもともに UTM Zone 54 にあることが分かりました。

UTM Zone 54 の EPSG コードを求める

ここで、geotools.rb を使って UTM Zone 54 の EPSG コードを探してみます。ここでは、GeoTools の EPSGCRSAuthorityFactory を裸で使います。あとでここの処理はラップする予定です。逆に言うと、geotools.rb の Geo モジュールのクラスやメソッドは以下のようなコードを足がかりに作っています。

require 'geotools'

f = Geo::Tools::EPSGCRSAuthorityFactory.new
list = f.getAuthorityCodes(Geo::Tools::ProjectedCRS)
iter = list.iterator
while(iter.hasNext)
  c = iter.next
  begin
    crs = f.createCoordinateReferenceSystem(c)
    name = crs.getName.toString
    print "code: #{c} => name: #{name}\n" if /UTM.*54/.match name
  rescue
    # GeoTools レベルで例外が発生するコードがあるので無視する
  end
end

このスクリプトを実行すると、

$ ruby epsg_search.rb 
DEBUG: rjb primitive_conversion mode
epsg_search.rb:10: warning: already initialized constant VERSION_KEY
code: EPSG:23894 => name: EPSG:ID74 / UTM zone 54S
code: EPSG:32254 => name: EPSG:WGS 72 / UTM zone 54N
code: EPSG:32354 => name: EPSG:WGS 72 / UTM zone 54S
code: EPSG:32454 => name: EPSG:WGS 72BE / UTM zone 54N
code: EPSG:32554 => name: EPSG:WGS 72BE / UTM zone 54S
code: EPSG:32654 => name: EPSG:WGS 84 / UTM zone 54N
code: EPSG:32754 => name: EPSG:WGS 84 / UTM zone 54S

となります。日本は北半球ですので、WGS 84 / UTM zone 54N がもっとも良くあてはまりますから、EPSG:32654 を使えばよいことになります。
このことは、http://bubble.atnifty.com/modules/bwiki/index.php?ProjParamList のちずろぐさんのところでも裏をとれます。

緯度経度の EPSG コードは?

日本における世界測地系は、厳密には WGS84 ではなく、JGD2000 ですから、以下のように検索してみたところ、EPSG:4612 を使えば良いことが分かります。

require 'geotools'

f = Geo::Tools::EPSGCRSAuthorityFactory.new
list = f.getAuthorityCodes(Geo::Tools::GeographicCRS)
iter = list.iterator
while(iter.hasNext)
  c = iter.next
  begin
    crs = f.createCoordinateReferenceSystem(c)
    name = crs.getName.toString
    print "code: #{c} => name: #{name}\n" if /JGD2000/.match name
  rescue
    # GeoTools レベルで例外が発生するコードがあるので無視する
  end
end
$ ruby epsg_search.rb 
DEBUG: rjb primitive_conversion mode
epsg_search.rb:10: warning: already initialized constant VERSION_KEY
code: EPSG:4612 => name: EPSG:JGD2000

このエントリの結論

  • 渋谷とつくばが含まれる座標系は、平面直角座標系 IX 系で、その EPSG コードは EPSG:2451
  • 渋谷とつくばが含まれる座標系は、UTM 54N で、その EPSG コードは EPSG:32654
  • 日本のいわゆる世界測地系は、JGD2000 で、その EPSG コードは EPSG:4612

ということが分かりました。肝心の距離計算は、次のエントリで。

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

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

geotools.rb

http://svgmapdata.sakura.ne.jp/geotools/ に公開されています。今回、org.opengis.referencing の座標系クラスをとりこむための変更を加えています。

*1:その場合、GeoTools では座標系の情報を英語でしか持っていないので、検索文字列は英語になってしまいます。Geo::search_epsg_from_string('Japan Planar IX') などといったインタフェースを考えることになると思います