久々に本職のシステム関連の投稿を・・・。

 ある日、配属先の同僚から次のような相談を受けました。
「今、ここに首都(ハボロネ)からの距離が500km、1000km、1500kmの同心円になっている地図があるんだけど、この中心をハボロネに限らず任意の点から作図出来るようなアプリ作れない?」
 地図上の距離に関しては以前に「緯度・経度情報のある2点間の距離」を計算する関数を作成しているので「ああ、出来るよ!」と気軽に答えたものの・・・細部にこだわったこともあり(相手はそこまで求めていないのですが)、結構苦労しました。

 まず、地図に関しては中心とする地点の選びやすさを考慮してGoogle Maps APIを利用することにしました。Google Mapにはgoogle.maps.Circleという便利なクラスがあり、これを使えば簡単に同心円が作れます。なので、円を描くだけならそれでいいのです。
 しかし、見本には各円に色が塗ってある。もちろんCircleクラスでも円を塗りつぶすことが出来ますが、複数の円を描く同心円では色が重なってしまいます。透過度をゼロにして内側から描けば色は重なりませんが、今度は地図が見えなくなってしまい場所が何処だか分からなくなってしまうので・・・。結局自分で計算することにしました。

 「点Aからの距離がrである点の座標(緯度・経度)」を求める計算(点Bは無数に存在)。最初考えたのはAから真東にだけ離れた点(同緯度)と真北にrだけ離れた点(同経度)を最初に求めて、その2点を通る楕円上の点の座標を求める方法です(下のイメージ図)。

点Aからの距離がrである点の座標(緯度・経度)を求める計算手法(1)~イメージ~

点Aからの距離がrである点の座標(緯度・経度)を求める計算手法(1)~イメージ~

 以前作成した2点間の距離を計算する関数を使えば、点Aから緯度1度分離れた点及び経度1度分離れた点の距離が計算出来るので、それを基準に上図の点B・点Cの位置情報を求め、後は三角関数を使用して楕円上の点を計算していくというものです。実際、ネットで調べてみると同じような理論で計算されている方もちらほらといらっしゃるので、これでOKかな?と目論んでいたものの・・・。
 以下のコードが当手法で計算しポリゴンとして地図上に落としたものです。

  var map;
  var mapCenter = new google.maps.LatLng(-23, 23);
  var pointMarker = new google.maps.Marker();
  var polygons = new Array();

  function initialize() {

     var myOptions = {
        zoom: 3,
        center: mapCenter,
        mapTypeControl: true,
        mapTypeId: google.maps.MapTypeId.ROADMAP
     };
     map = new google.maps.Map(document.getElementById("map_canvas"), myOptions);

     google.maps.event.addListener(map, "click", function(clickevent) {
        var myPoint = new google.maps.LatLng(clickevent.latLng.lat(), clickevent.latLng.lng());
        setPointMarker(myPoint);
     });


  }

  function setPointMarker(location, pile) {

    if(!pointMarker.getMap()){
       pointMarker = new google.maps.Marker({
          position: location,
          map: map,
          draggable: true,
       });
    }else{
       pointMarker.setPosition(location);
    }
  }

  function drawCircle() {

     if (!pointMarker.getMap()) {
        alert("中心となる点にマーカーを設置して下さい。");
        return 0;
     }

     var myPoint = pointMarker.getPosition();
     var pN = new google.maps.LatLng(myPoint.lat()-0.5,myPoint.lng()); //中心点から北へ0.5度
     var pS = new google.maps.LatLng(myPoint.lat()+0.5,myPoint.lng()); //中心点から南へ0.5度
     var pW = new google.maps.LatLng(myPoint.lat(),myPoint.lng()-0.5); //中心点から西へ0.5度
     var pE = new google.maps.LatLng(myPoint.lat(),myPoint.lng()+0.5); //中心点から東へ0.5度
     var dLat = getDistance(pN,pS); //緯度1度分の距離
     var dLng = getDistance(pW,pE); //距離1度分の距離
     var curLat,curLng;
     var cLat,cLng;
     var polyCircle = new Array;
     var r = 3000000;//今回は3000kmで計算
     var degpitch = 0.5; //今回は同心円上の点は0.5度ピッチで計算

     for (var theta = 0; theta < 360; theta += degpitch) {
        curLng = myPoint.lng() + (r / dLng) *  Math.cos(calRadian(theta));  //経度(X座標)
        curLat = myPoint.lat() + (r / dLat) *  Math.sin(calRadian(theta));  //緯度(Y座標)
        polyCircle.push(new google.maps.LatLng(curLat,curLng));
     }
     polyCircle.push(polyCircle[0]);

     var myCircle = new google.maps.Polygon({
        paths: polyCircle,
        fillColor: '#00FFFF',
        fillOpacity: 0.5,
        strokeColor: '#00FFFF',
        strokeOpacity: 1.0,
        strokeWeight: 1
     });     
     myCircle.setMap(map);

     //比較用 Google Maps API Circleクラスでの描画
     var myCircleB = new google.maps.Circle({
        map: map,
        center: myPoint,
        radius: r,
        fillOpacity: 0,
        strokeColor: "#FF0000",
        strokeOpacity: 1,
        strokeWeight: 2
     });

  }

  //度数→ラジアン変換
  function calRadian(deg){
     return deg * Math.PI / 180.0;
  }

  //2点間の距離を計算(GRS80)
  function getDistance(p1, p2) {
     var a = 6378137.0;             //長半径
     var f = 1.0 / 298.257222101;   //扁平率
     var b = a * (1 - f);           //短半径
     var e2 = (Math.pow(a, 2) - Math.pow(b, 2)) / Math.pow(a, 2);
     var dy = calRadian(p1.lat() - p2.lat());
     var dx = calRadian(p1.lng() - p2.lng());
     var my = calRadian((p1.lat() + p2.lat()) / 2.0);
     var w = Math.sqrt(1.0 - e2 * Math.pow(Math.sin(my), 2));
     var m = (a * (1 - e2)) / Math.pow(w, 3);
     var n = a / w;

     return Math.sqrt(Math.pow(dy * m, 2) + Math.pow(dx * n * Math.cos(my), 2));
  }

 一見良さげなのですが、Circleクラスで描いた円(下図の赤い線)と比較すると結構ズレがあります。
 真東・真西の位置はぴったり合っているようなのですが、特に南東・南西方向のズレは大きく100km以上です(下の拡大図参照)。これだと都市地図等の近距離なら気にならない程度の誤差で収まるかもしれませんが、世界地図レベルとなると実用には耐えられそうもありません。

計算・描画結果

計算・描画結果 ボツワナの首都(ハボロネ)から3000kmの円
Google Maps APIのCircleクラス利用で描いたもの(赤い線)とは結構ズレがあります。
この地図を別窓で表示


誤差の大きい部分を拡大(中心から南東方向)

ズレの大きい部分を拡大(中心から南東方向)。100km以上の誤差があります。これではちょっと・・・


 考えてみたら、経度1度当たりの距離は高緯度になるほど短くなるのに、それを考慮せずに1度当たりの距離を過大評価してしまっている状態ですから、等距離の円も高緯度では実際よりも萎んだ形となってしまうのも当然です。
 また、地球は真球ではなく楕円体なので同じ「緯度1度分」であっても北へ1度と南へ1度とでは距離が違うのですね。つまり、緯度差と距離は比例しないので比率換算では緯度は求められない(上記の方法はそこが間違い)。
 参考までに、上記の2点間の距離を計算する関数を利用して、東京から北へ10度離れた点と南へ10度離れた点までの距離を調べてみたところ、

 ・東京(35.69,139.69)から北へ1度離れた点までの距離:1110479.2m
 ・東京(35.69,139.69)から南へ1度離れた点までの距離:1108641.5m

 2km近い距離差となりました。意外だったのは北側の方が長かったことです。地球は赤道半径の方が極半径より長いので、「低緯度側(北半球の場合は南側)の方が当然長くなるだろうと計算前には思っていたのですが結果は逆でした。プログラムが間違っているのかとも思いましたが、国土地理院のサイトでの計算でもほぼ同じ結果を得られたので間違いないようです。

 緯度は地球の中心からの角度で決めたもの(地心緯度)ではなく、楕円体面の法線と赤道の角度で決めたもの(地理緯度)だそうで、そのため高緯度ほど緯度1度当たりの距離が長くなるようで・・・、私にはなんだか難しいですが、いずれにしても上記の方法では駄目なので、考え直しが必要です。次回へつづく

コメント   

 プラ(Pula)とはボツワナの現地語(ツワナ語)で雨のことです。
 ケッペンの気候区分によればボツワナの国土の大半はステップ気候に属しており(一部は砂漠気候)、雨が少なく作物が育ちにくい環境にあります。

ケッペンの気候区分(ボツワナ)

ケッペンの気候区分(ボツワナ)


 それだけに、たまに降る雨(プラ)はとても貴重で天からの贈り物。通貨がプラであることや、国旗のデザインに水色が採用されていることからも、ボツワナの人々にとって雨がいかに大切であるかが分ります。
100プラ紙幣

通貨もプラ。上写真は100プラ紙幣(約1,000円)


ボツワナ国旗

ボツワナ国旗。水色は雨(水)を表します。

 10月から4月くらいまでが雨季に当たりますが、首都ハボロネの年間降水量は400㎜程度と東京の4分の1ほどで、雨季といってもそんなに降るわけではありません(だからステップ気候なわけですが)。
 しかし、今シーズンは雨が多いです。ラニーニャの影響だそうですが、特にこの2ヶ月くらいは雨が降らない日の方が多いです。昨年は干ばつで農業が大打撃を受けたようなので、こちらの人々は「今年は雨が多くていいね!」と言っていますが、もともと雨が少ない国なので排水の設備が整っておらず、各所で水害が発生しています。
 私は自転車通勤なのでレインコート必携ですが、雨のおかげでさほど暑くならないのは助かります。昨年の夏は記録的な暑さで気温46度の日も何日かあったそうなので・・・。

ここ最近は雨ばかり

ここ最近は雨ばかり

コメント