LineString 近辺の点を拾い出すには、 distance を使う方が buffer + within よりも 10 倍速い
distance と within、どちらが速い？
ある LineString から一定距離 d の中にある点を拾い出す処理を JTS / OpenGIS / ISO 19100 の Geometry を使って実現するには、例えば以下のような方法が考えられます。
- [distance 法]LineString と対象点の距離を計測し、その距離が d 未満であればその点を拾う。
- [buffer 法]LineString に距離 d のバッファを発生させ、対象点がそのバッファへ包含されればその点を拾う。
どちらが速いのでしょうか。また、速さにはどのくらいの違いが生じるのでしょうか。
distance 法が速い？
常識的に考えると、distance 法のほうが速いように予想されます。その理由として、以下が考えられます：
- 次元が違う。distance 法は 0次元対1次元の演算であるが、buffer 法は 1次元対2次元の演算になる。次元があがることにより、計算は飛躍的に複雑になるのではないか。
- バッファ作成により、計算に使う図形の点数も増加してしまう。バッファの円を多角形で近似するためである。点数が増加すれば、自ずと計算時間も増大してしまう。
buffer 法が速い？
buffer 法が速いとすれば、その理由としては以下が考えられます。
- within は面が点を内包しているかどうか、真偽値を返せばよいが、distance は距離という数値を返さなければならない。距離を計算するには、例えば開平などの演算が必要と想像されるが、内包の演算ではこのあたりの演算を省略でき、それによる速度向上があるかもしれない。
- LineString が複雑である場合には、バッファ図形の方が LineString よりも単純である場合があるかもしれない。その場合には内包演算の方が速いかもしれない。
実験
地球地図日本 Version 1.1 の transl_1_1.shp と hydrol_1_1.shp を入力データとして、各 LineString データの bounding box 内に 10000 点のランダムな点を発生し、その点が LineString から 10m 以内*1にあるかどうかを判定する時間を計算させました。
距離が関わる計算なので、LineString の重心の位置に基づいて UTM 投影した上で諸計算をしています。投影の時間は計測時間から除外しています。
distance 法を採用することによる速度の向上（buffer 法の処理時間を distance 法の処理時間で割ったもの）を、LineString の長さに対してプロットすると、以下のようになりました。
どの長さにおいても、10倍程度の速度の違いが発生しているようです。
考察
ArcGIS などの GUI ベースの GIS では、ここで buffer 法と呼んだ、
近傍点の検索をするのに、バッファ付きポリゴンを計算してからバッファ付きポリゴンへの内包演算をする
という方便がよく使われています。これは GUI という UI によって計算の遅さがそれほど気にならないためであると思われます。しかし、GUI の GIS でよく使われている方法が必ずしも計算上効率的であるとは限らないようであることが今回分かりました。
ArcGIS で良くやっているからといって、buffer 法でプログラムを書くと、計算時間を 10 倍のオーダーで損する可能性があります。
ソースコード
実験に用いたスクリプトは、以下のようなものです。
require 'geotools' src = '(...){transl_1_1.shp|hydrol_1_1.shp}' n_trials = 10000 buf_width = 10 wgs84 = Geo::import_epsg_crs(4326) utms = {} trs = {} 51.upto(56) do |i| utms[i] = Geo::import_epsg_crs(32600 + i) trs[i] = Geo::Transform.new(wgs84, utms[i]) end wf = File.open('buffer_or_distance.csv', 'w') Geo::Reader.foreach(src, {:whitelist => []}) do |geom, attrs| # projection geom = geom.getGeometryN(0) next if geom.getLength == 0 # because there are such strange data... c = geom.getCentroid tr = trs[((c.getX + 180) / 6).ceil] geom = tr.transform(geom) n_points = geom.getNumPoints env = geom.getEnvelopeInternal points = [] n_trials.times do |i| points << Geo::import_array_geometry([ rand(nil) * env.getWidth + env.getMinX, rand(nil) * env.getHeight + env.getMinY]) end start_time = Time.now buf = geom.buffer(buf_width) buffering_time = Time.now - start_time count_buffer = 0 start_time = Time.now points.each do |point| count_buffer += 1 if point.within(buf) end process_time_buffer = Time.now - start_time print "#{count_buffer} points are in by within (#{process_time_buffer}s)\n" count_distance = 0 start_time = Time.now points.each do |point| count_distance += 1 if point.distance(geom) < buf_width end process_time_distance = Time.now - start_time print "#{count_distance} points are in by distance (#{process_time_distance}s)\n" print "n_points = #{n_points}\n\n" wf.print "#{geom.getLength},#{env.getWidth},#{env.getHeight},#{env.getWidth / env.getHeight},#{process_time_buffer},#{process_time_distance},#{process_time_buffer / process_time_distance},#{count_distance.to_f / n_trials},#{n_points},#{count_buffer},#{count_distance},#{count_buffer == count_distance}\n" end wf.close
