2016年7月3日日曜日

OpenLayers 3 を使ってみよう(その17:OpenLayers v3 で 経路の距離を測る)

 これは,OpenLayers 3 を使ってみよう(その16:OpenLayers v3 で Google Map と地理院地図を重ねる)からの続きになる。 OpenLayers 3 を使ってみよう(その0:はじめに:地理院地図を表示)に目次がある。
 その13までのページでは OpenLayers v3.7.0 で書いてきたが,ここでは OpenLayers 3.16.0 で書かれている。
 前回(その16:OpenLayers v3 で Google Map と地理院地図を重ねる)は,Google Map と地理院地図を重ねてみた。
今回は少し趣きを変えて,経路の距離を測る話について書こう。

 地図上の経路の距離に関しては,その5:テキストデータから折線データ読込みに少し書いた。 そこでは,ヒュベニの公式と呼ばれるものを使っていると書いた。 その部分のみ書き出すと以下のようになる(ここではその13:センターマークと凡例の表示・KML による経路+マーカーの場合でのルーチンを取り出している)。
// ヒュベニの公式で緯度・経度から距離を求めるための定数
var long_r = 6378137.000;     // [m] 長半径
var short_r = 6356752.314245; // [m] 短半径
var rishin = Math.sqrt((long_r * long_r - short_r * short_r)/(long_r * long_r)); // 第一離心率
var a_e_2 = long_r * (1-rishin * rishin);  // a(1-e^2)
var pi = 3.14159265358979;    // Pi
// *******************************************************************
// ヒュベニの公式を使った距離計算
// 2点間の距離
function dist_2pts(lon0, lat0, lon1, lat1) {
    lon0 = lon0 * pi / 180;  lat0 = lat0 * pi / 180; // in radian
    lon1 = lon1 * pi / 180;  lat1 = lat1 * pi / 180; // in radian
    var d_lon = lon1 - lon0;
    var d_lat = lat1 - lat0;
    var ave_lat = (lat1+lat0)/2;
    var Wx = Math.sqrt(1-rishin * rishin * Math.sin(ave_lat) * Math.sin(ave_lat));
    var Mx = a_e_2 /Wx/Wx/Wx;
    var Nx = long_r /Wx;
    var dum = (d_lat * Mx)*(d_lat * Mx) + (d_lon* Nx * Math.cos(ave_lat)) * (d_lon* Nx * Math.cos(ave_lat)); // square of distance
    return Math.sqrt(dum);
}
// ===================================================================
// 経路の長さを求めて表示する
function length_line() {
    var lineLength = 0;
    var features = kml_vector.getSource().getFeatures(); // 通常の JavaScript の array
    for (var j=0; j < features.length; j++) {
        var geometry = features[j].getGeometry(); // featureOverlay の中のgeometry(例:ol.geom.LineString )
        var coordArray = geometry.getCoordinates(); // geometry 変数の中の座標を配列へ(メルカトル座標系)
        for (var i=0; i < coordArray.length; i++) { coordArray[i] = ol.proj.transform(coordArray[i],"EPSG:3857", "EPSG:4326"); }
        for (var i=0; i < (coordArray.length-1); i++) { lineLength = lineLength + dist_2pts(coordArray[i][0],coordArray[i][1],coordArray[i+1][0],coordArray[i+1][1]); }
    }
    document.getElementById("outStr").innerHTML = "  L = "+(Math.floor(lineLength)/1000)+" [km]";
}
 ここで,long とか lon と書いてあるのは,longitude で経度を表している。 また,lat は latitude で緯度を表している。 また,coordArray という配列があるが,それは,[[longitude1, latitude1], [longitude2, latitude2], .....] という2次元配列であり,経路の経度・緯度を表している。

 ヒュベニの公式では,地球を真球ではなく楕円体だと思って距離を求めるものである。 ヒュベニの公式については,GPSのログを間引く(折れ線を間引く)に少し書いたのだが, 私は二地点の緯度・経度からその距離を計算するにある記述通りの計算をさせている。 ヒュベニの公式の精度に関しては,[Mathematica] GeoDistanceとその他の測地線距離算出式の精度によると(そこでの記述が正しいと思って引用している),日本国内に限れば 6桁,ないし7桁の精度はありそうということみたい。 6桁ということは,距離計算の誤差が 100 万分の1程度,ということであり,100 km の距離を求めても 0.1 m 程度の誤差しかない,ということになる。 まぁ,それぐらいなら実用的には全く問題ないと思っている(北極や南極に近づくほど精度は悪くなるらしい)。

 OpenLayers 3 のルーチンの中にもいくつか点の間の距離を求めるものがある。 例えば,ol.geom.LineString クラスに対するメソッド(関数みたいなもの)として,
getLength()
というのがある。 その5:テキストデータから折線データ読込みには書いたのだが,この getLength() を使うと,本来 1592 km の札幌―鹿児島間の距離が 2012 km となってしまい使い物にならなかった。 どうもこれは適切な projection(投影方法)を選ばないといけないみたいなのだが,それ以降まじめに調べていない。

 しかし,この間,OpenLayers の Exampleのリストを見ていると,Measure の Exampleがあるのに気づいた。そこでは,経路を連続する点として定義して,経路の距離を計算する,というものだった。 ソースを見ると,吹き出し(tooltip)を出したり,経路を入力させたりしているので,かなり複雑だったが, ポイントは,半径 6378.137 km の真球(var wgs84Sphere = new ol.Sphere(6378137);として定義)に対して,
wgs84Sphere.haversineDistance(c1, c2);
というルーチンを使って2点間の距離を求める,というものだった。
 ここで,haversine というのは球面三角法の半正矢関数というものらしい。 Wiki 球面三角法によると,
hav x = (1-cos x)/2
というものらしい。 この haversine 関数を使って2点間の距離をだすらしい。

 そこで,この haversineDistance ルーチンを使って距離を求めてみた。 やり方としては,上記のヒュベニの公式を使って経路の距離を求めるルーチンの最後の関数 function length_line() を以下のように変えてみた。
//===================================================================
// 経路の長さを求めて表示する
function length_line() {
    var wgs84Sphere = new ol.Sphere(6378137);
    var measureLength = 0;
    var lineLength = 0;
    var features = kml_vector.getSource().getFeatures(); // 通常の JavaScript の array
    for (var j=0; j < features.length; j++) {
        var geometry = features[j].getGeometry(); // featureOverlay の中のgeometry(例:ol.geom.LineString)
        var coordArray = geometry.getCoordinates(); // geometry 変数の中の座標を配列へ(メルカトル座標系)
        for (var i=0; i < coordArray.length; i++) { coordArray[i] = ol.proj.transform(coordArray[i],"EPSG:3857", "EPSG:4326"); }
        for (var i=0; i < (coordArray.length-1); i++) { lineLength = lineLength + dist_2pts(coordArray[i][0],coordArray[i][1],coordArray[i+1][0],coordArray[i+1][1]); }
        for (var i=0; i < (coordArray.length-1); i++) { measureLength = measureLength + wgs84Sphere.haversineDistance(coordArray[i], coordArray[i+1]); }
    }
    var outstr;
    outstr = (Math.floor(measureLength)/1000)+" [km]";
    document.getElementById("outStr").innerHTML = "  L = "+(Math.floor(lineLength)/1000)+" [km]"+" (ol.measure "+outstr+")";
}
 ここでは,(複数あると仮定した)全経路の経路点のデータを取り出し,経路ごとに coordArray という点の配列に入れ,球面メルカトル図法で与えられている座標を通常の緯度経度による座標系 (WGS 84) に変換したのち,ヒュベニの公式を使った計算 (dist_2pts 関数) と,球面三角法の haversine 関数を用いたルーチンで別々に計算させて,結果を求めている。

 このルーチンを使ったサンプルを「その17」のサンプルとして用意してみた(サンプルページはアクセスログを取るルーチンを組み込んでいます)。 今回の経路は塩屋を始点とした六甲全山縦走の経路である。 このサンプルで,右上にある「経路長の表示」というボタンを押すと,
L = 44.212 [km] (ol.measure 44.257 [km])
という表示がでるはずである。 これは,ヒュベニの公式で求めた経路の距離が 44.212 km だったのに対して,OpenLayers の haversine 関数を用いた計算が 44.257 km になった,というものである。 これを見る限り,3桁目までは同じであり,まずまずの一致を示していると言える。 しかし,楕円体か真球かという違いが 0.1 % の誤差として出てきている。 ヒュベニの公式を用いる場合と haversine 関数を用いる場合で,記述はそんなに違わない(と思っている)ので,個人的にはヒュベニの公式を使おうと思っているが, 面倒だと思われる方は(多少ごちゃごちゃ書かないといけないので),haversine 関数を使うルーチンを使ってみてはいかがだろうか???

その18:Marker アニメーションに続く

0 件のコメント: