前回までで、円の大きさが地球を１周するような大きなものにならない限りは等距離の円を同心円状に描くことが出来るようになりました。
今回は、
・円の半径を地球を１周するような規模に設定した場合、本来とは逆の領域が選択されてしまう場合がある
この不具合に対応したいと思います。
ポリゴンを描く際に選択する領域をA、選択していない方の領域をBとした時、Google Maps APIではAとBの面積の差があまりないような場合には、A・Bのうち北側に位置する方の領域を選択する傾向にあるようです。そのため、南半球にある地点を円の中心点にする場合はAの方が南側になりますからBが選択されてしまいます。また、例えば東京から15000kmの等距離円などBの面積が明らかにAよりも小さくなる場合にはBが選択されます。
google.maps.Polygonメソッドに「指定した緯度・経度がポリゴン内部に含まれるように領域選択する」オプションがあればいいのになあ・・・と思うのですが、残念ながら備わっていません。
でも、そのオプションの有無を探している際に、google.maps.geometryライブラリのポリゴン関数にcontainsLocationというものがあるのを見つけました。
・Maps Javescript API ジオメトリ ライブラリより抜粋
containsLocation(point:LatLng, polygon:Polygon)
この関数は、指定の地点がポリゴンの内部または境界に存在すれば、true を返します。
これを利用すれば、次のようなロジックでポリゴン描画すれば、上記ケースにも対応出来るのではないかと考えました。
1. 半径の大きい円のポリゴンをダミーで描く（visible=falseとして表示させない）
2. 円の中心点がダミーポリゴンの内部に含まれているかどうかをcontainsLocation関数で調べる
3. 内部に含まれている場合はそのままでOKなのでダミーと同じものをもう一度描き表示させる。
4. 内部に含まれていない場合は逆の領域Bがポリゴンとして選択されてしまったことになる。
領域Aを選択させるには、一旦地球（地図）の全領域を覆うポリゴンを描いた後、領域Bを中抜け
領域としてくり抜いてしまえば良い。
コードにすると次のようになります
//地図全体を覆うポリゴン用座標データ（反時計回り）
var polyWholeWorld = [
new google.maps.LatLng(0, 180),
new google.maps.LatLng(-90, 180),
new google.maps.LatLng(-90, 90),
new google.maps.LatLng(-90, 0),
new google.maps.LatLng(-90, -90),
new google.maps.LatLng(-90, -180),
new google.maps.LatLng(0, -180),
new google.maps.LatLng(90, -180),
new google.maps.LatLng(90, -90),
new google.maps.LatLng(90, 0),
new google.maps.LatLng(90, 90),
new google.maps.LatLng(90, 180),
new google.maps.LatLng(0, 180)
];
//ダミーのポリゴンを作成
var dummyCircle = new google.maps.Polygon({
paths: polyCoords, //ポリゴン領域の位置情報（緯度・経度）の配列（時計回り）
visible: false
});
//中心点（myPoint）がダミーポリゴン内部に含まれているかどうかをチェック
var innerFlag = google.maps.geometry.poly.containsLocation(myPoint, dummyCircle);
if (innerFlag) {
//内部にある場合は、そのままのデータでOK
polyPath = [polyCoords];
} else {
//内部にない場合は、地球全領域からダミーポリゴンの領域をくり抜く
polyPath = [polyWholeWorld, polyCoords]; //polyWholeWorldとpolyCoordsは回りが逆
}
//等距離円の描画
myCircle = new google.maps.Polygon({
paths: polyPath,
strokeWeight: 0,
fillColor: red,
fillOpacity: 0.4
});
myCircle.setMap(map);
後は前回同様に内側の円領域を切り抜く処理を加える必要があります。通常のケースでは単に１つ内側の円領域データ（配列）を逆順にすれば良いのですが、「内側の円の段階で本来とは逆の領域が選択されてしまっている場合にはどうすればいいのだろうか？」という疑問が浮かびましたが、このケースでも正しい領域を描くために使用したデータ（地球全体を選択して領域Bをくり抜くパターン）をそのまま全て逆順にする次のような手法で上手くいきました。
地球全体領域－外円の逆領域（領域B1）－｛地球全体領域－内円の逆領域（領域B2）｝
この式を展開すれば
内円の逆領域（領域B1）－外円の逆領域（領域B1）
これで良かったのか・・・でも考え方としては上式の方が分かりやすい気がするので、結果は同じですが今回は上式を採用しました。
計算・描画結果 ボツワナの首都（ハボロネ）から3000km,5000km,8000km,10000kmの同心円
上手く出来たようです。しかし、ボツワナから10,000kmでは日本に全然届かない・・・やっぱり遠いなあ。
前回、一応それっぽい同心円を描いてみたものの、等距離円の距離が世界地図規模（千km単位）となると計算誤差が大きくなるため、今回はもう少し正確性の高い計算方法を探索します。
以前に位置情報（緯度・経度）の分かっている２点間の距離を知りたいと思った時には、その求め方が記載されたサイトはWeb上に多数存在し、結構簡単に見つけられたのですが、今回必要としている「ある１点からの距離がｒである点の緯度・経度を求める方法」はなかなか見つけられず、一番頼りになる国土地理院のサイトにも記載なし。「もうここはおとなしくgoogle.maps.Circle利用で妥協するか（同心円の色が重なってしまいますが）」と諦めかけたところで漸く見つけました・・・Vincenty法（順解法）。一旦見つけてしまえば「何で今まで見つけられなかったんだろう？」と思うのですが、探し物にはよくあることですね。
始点の緯度・経度、始点からの距離・方位角を、このVincenty法（順解法）の式に代入することで到達点の緯度・経度が求まります。
同心円を描くには始点の周り360度の点を全て求める必要があるので、以下のように関数化しました。
//Vincenty法（順解法）
function VincentyDirect(p1,deg1,s) {
/*
引数
p1 : 中心となる点（点１）の緯度・経度（google.maps.LatLngクラス）
deg1 : 求める点の方角（真北からの時計回りの角度。北:0、東:90、南：180、西：270）
s : 点１からの距離（ｍ）
*/
var p2 //求める点（点２）の緯度・経度（google.maps.LatLngクラス）
var a = 6378137.0; //長半径（赤道半径）
var f = 1.0 / 298.257223563; //扁平率（WGS84）
var b = a * (1 - f); //短半径（極半径）
var phi1,phi2; //点１、点２の緯度
var L1,L2; //点１、点２の経度
var L; //点１、点２の経度差
var alpha1,alpha2;
var alpha;
var U1,U2;
var lambda
var sigma;
var sigma1,formersigma;
var sigmam,deltasigma;
var epsilon = 1.0e-12; //計算精度（sigmaの許容誤差）
var sinalpha,cosalpha2;
var u2,A,B,C;
var cnt=0;
phi1 = calRadian(p1.lat());
L1 = calRadian(p1.lng());
alpha1 = calRadian(deg1);
U1 = Math.atan((1 - f)*Math.tan(phi1));
sigma1 = Math.atan2(Math.tan(U1),Math.cos(alpha1));
sinalpha = Math.cos(U1)*Math.sin(alpha1);
cosalpha2 = 1-Math.pow(sinalpha,2);
u2 = cosalpha2*((Math.pow(a,2)-Math.pow(b,2))/Math.pow(b,2));
A = 1+u2/16384*(4096+u2*(-768+u2*(320-175*u2)));
B = u2/1024*(256+u2*(-128+u2*(74-47*u2)));
sigma = s/(b*A);
do {
sigmam = (2*sigma1+sigma)/2;
deltasigma = B*Math.sin(sigma)
*(Math.cos(2*sigmam)
+1/4*B*(Math.cos(sigma)*(-1+2*Math.pow(Math.cos(2*sigmam),2))
-1/6*B*Math.cos(2*sigmam)*(-3+4*Math.pow(Math.sin(sigma),2))*(-3+4*Math.pow(Math.cos(2*sigmam),2))));
formersigma = sigma;
sigma = s/(b*A)+deltasigma;
cnt++;
if (cnt>100) {
alert("error");
return null;
}
} while (Math.abs(sigma-formersigma)>epsilon);
phi2 = Math.atan2((Math.sin(U1)*Math.cos(sigma)+Math.cos(U1)*Math.sin(sigma)*Math.cos(alpha1))
,((1-f)*Math.sqrt(Math.pow(sinalpha,2)+Math.pow((Math.sin(U1)*Math.sin(sigma)-Math.cos(U1)*Math.cos(sigma)*Math.cos(alpha1)),2))));
lambda = Math.atan2(Math.sin(sigma)*Math.sin(alpha1),(Math.cos(U1)*Math.cos(sigma)-Math.sin(U1)*Math.sin(sigma)*Math.cos(alpha1)));
C = f/16*cosalpha2*(4+f*(4-3*cosalpha2));
L = lambda-(1-C)*f*sinalpha*(sigma+C*Math.sin(sigma)*(Math.cos(2*sigmam)+C*Math.cos(sigma)*(-1+2*Math.pow(Math.cos(2*sigmam),2))));
L2 = L + L1;
alpha2 = Math.atan(sinalpha/(-1*Math.sin(U1)*Math.sin(sigma)+Math.cos(U1)*Math.cos(sigma)*Math.cos(alpha1)));
p2 = new google.maps.LatLng(calDegree(phi2),calDegree(L2));
return p2;
}
アークタンジェントの関数はMath.atan()とMath.atan2()の２種類あり、これらを使い分ける必要がありますが、Wikipediaではその辺も親切に解説されています。
Vincenty法での緯度・経度を計算するサイト（数箇所）での計算結果と比較することで検証しましたが、誤差はほとんどなかったので上記ソースで多分大丈夫かと思います。なお、扁平率はWGS84のものを採用しましたが、GRS80の扁平率を採用しても計算結果はほとんど変わりません。
と、ここまで作成した段階で残念な発見・・・同様の計算をする関数がgoogle.maps.geometryライブラリに用意されていました。
・Maps Javescript API ジオメトリ ライブラリより抜粋
computeOffset() を使用すると、特定の角度、出発地点、移動距離（メートル単位）から、目的地の座標を計算できます。
Google Maps APIを読み込む際に、
<script type="text/javascript" src="http://maps.google.com/maps/api/js?language=jp&key=取得したAPIキー&libraries=geometry"></script>
といった具合にlibraries=geometryのパラメータを追記するだけで利用可能です。google.maps.Circleも当関数を利用しているようです。
computeOffset()と自作の上記関数とで誤差があるか検証したところ、低緯度の鉛直方向（北半球なら南側）の誤差が大きくなるようで、ハボロネ（ボツワナの首都）から5000km地点計算時には約31kmのズレが出ました。5000kmで31kmならば誤差は1%未満、地図上ではほとんど分からないレベルですが、原因が気になるのでいろいろ調べてみるとgoogle.maps.geometryライブラリでは地球を完全な球体とみなして計算していることが分かりました（自作関数の扁平率をゼロにして再度検証した結果は両者一致しました）。
大差ないとはいえ、自作関数の方が扁平率を考慮している分だけより正確ということになるのか・・・わざわざ作った甲斐がありました。computeOffset()ではなく、こちらを採用します。
ここまで出来てしまえば、地図に落とす処理は簡単です。
今回のテスト用html、Javascriptのソースを以下に記載します。
・concentriccircles_test2.html
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
<html>
<head>
<title>concentric circles on the map</title>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<meta name="viewport" content="initial-scale=1.0, user-scalable=no" />
<meta http-equiv="Content-Style-Type" content="text/css">
<script type="text/javascript" src="http://maps.google.com/maps/api/js?language=jp&libraries=geometry"></script>
<style type="text/css">
<!--
.circletitle {
font-size:12px;
color:#FFFF00;
padding-left:20px;
white-space: nowrap;
}
.circletitle input {
width:40px;
}
.circletitle input[type="button"] {
width:60px;
font-size:14px;
}
-->
</style>
<script type="text/javascript" src="js/concentriccircles_test2.js"></script>
</head>
<body style="margin:0px;" onload="initialize()">
<div id="top_bar" style="background-color:#333399;padding:5px;">
<table>
<tr>
<th style="color:#FFFFFF;text-align:left;padding-left:20px;" colspan="7" style="font-size: 16px;">中心となる地点から一定距離にある点を求めて地図上に同心円を描く（テスト２）～Vincentyの順解法利用～</th>
</tr>
<tr>
<td class="circletitle">circle1</td>
<td class="circletitle">circle2</td>
<td class="circletitle">circle3</td>
<td class="circletitle">circle4</td>
<td class="circletitle">circle5</td>
<td colspan="2"> </td>
</tr>
<tr>
<td class="circletitle">
<input type='text' id="f_distance0">km
<input type='text' id="f_color0">
</td>
<td class="circletitle">
<input type='text' id="f_distance1">km
<input type='text' id="f_color1">
</td>
<td class="circletitle">
<input type='text' id="f_distance2">km
<input type='text' id="f_color2">
</td>
<td class="circletitle">
<input type='text' id="f_distance3">km
<input type='text' id="f_color3">
</td>
<td class="circletitle">
<input type='text' id="f_distance4">km
<input type='text' id="f_color4">
</td>
<td class="circletitle">
<input type="button" id="btnDrawCircle" name="btnDrawCircle" onclick="drawCircle()" value="draw">
</td>
<td class="circletitle">
<input type="button" id="btnResetCircle" name="btnResetCircle" style="color:#ff0000;" onclick="resetCircle()" value="reset">
</td>
</tr>
</table>
</div>
<div id="map_canvas" style="margin:5px;height:600px"></div>
</body>
</html>
・concentriccircles_test2.js
var map;
var svg;
var Gaborone = new google.maps.LatLng(-24.6091799,25.7900739);
var iniZoom = 4;
var pointMarker = new google.maps.Marker();
var circleInfo = new Array();
var myCircle = new Array();
var polygons = new Array();
var drawFlag;
function initialize() {
drawFlag = false;
var dColor = ["blue","green","red","yellow","magenta"];
var dDistance = [500,1000,1500];
for (var i = 0; i <= 4; i++) {
if (dDistance[i]) {
document.getElementById("f_distance"+i).value = dDistance[i];
}
if (dColor[i]) {
document.getElementById("f_color"+i).value = dColor[i];
}
}
var myOptions = {
zoom: iniZoom,
center: Gaborone,
mapTypeControl: false,
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 setCircleInfo() {
var curDistance,curColor,curPplText;
circleInfo = new Array();
for (var i = 0; i <= 4; i++) {
curDistance = parseFloat(document.getElementById("f_distance"+i).value)*1000;
curColor = document.getElementById("f_color"+i).value;
if (!curColor) {
curColor = "#ffffff";
}
if (isFinite(curDistance)) {
circleInfo.push({"distance":curDistance,"color":curColor});
}
}
circleInfo.sort(function(a,b){
if(a.distance<b.distance) return -1;
if(a.distance > b.distance) return 1;
return 0;
});
}
function resetCircle() {
var i;
for (var i = 0; i < myCircle.length; i++) {
if (myCircle[i]) {
myCircle[i].setMap(null);
}
}
myCircle.length = 0;
if (pointMarker.getMap()) {
pointMarker.setMap(null);
}
drawFlag = false;
}
function drawCircle() {
var polyCoords = new Array();
var polyPath = new Array();
var degpitch = 0.5;//0.5度ピッチで計算
if (!pointMarker.getMap()) {
alert("中心となる点にマーカーを設置して下さい。");
return 0;
}
if (drawFlag) {
return 0;
}
var myPoint = pointMarker.getPosition();
setCircleInfo();
for (var i = 0; i < circleInfo.length; i++) {
var polyCoord = new Array();
//中心点から指定の距離だけ離れた地点の座標を計算し、配列に収納（ポリゴン用）
for (var theta = 0; theta < 360; theta += degpitch) { //thetaは真北からの角度（時計回り）
polyCoord.push(VincentyDirect(myPoint,theta,circleInfo[i].distance));
}
polyCoord.push(polyCoord[0]);
polyCoords[i] = polyCoord;
polyPath = [polyCoords[i]];
//最内以外の円描画時は、１つ内側の円の座標を反対周りにした配列を付加することでドーナツ型ポリゴンにする
if (i > 0) {
polyPath.push(polyCoords[i-1].slice().reverse());
}
myCircle[i] = new google.maps.Polygon({
paths: polyPath,
strokeWeight: 0,
fillColor: circleInfo[i].color,
fillOpacity: 0.4
});
myCircle[i].setMap(map);
}
drawFlag = true;
}
function calRadian(deg){
return deg * Math.PI / 180.0;
}
function calDegree(rad){
return rad * 180 / Math.PI;
}
//Vincenty法（順解法）
function VincentyDirect(p1,deg1,s) {
/*
引数
p1 : 中心となる点（点１）の緯度・経度（google.maps.LatLngクラス）
deg1 : 求める点の方角（真北からの時計回りの角度。北:0、東:90、南：180、西：270）
s : 点１からの距離（ｍ）
*/
var p2 //求める点（点２）の緯度・経度（google.maps.LatLngクラス）
var a = 6378137.0; //長半径（赤道半径）
var f = 1.0 / 298.257223563; //扁平率（WGS84）
var b = a * (1 - f); //短半径（極半径）
var phi1,phi2; //点１、点２の緯度
var L1,L2; //点１、点２の経度
var L; //点１、点２の経度差
var alpha1,alpha2;
var alpha;
var U1,U2;
var lambda
var sigma;
var sigma1,formersigma;
var sigmam,deltasigma;
var epsilon = 1.0e-12; //計算精度（sigmaの許容誤差）
var sinalpha,cosalpha2;
var u2,A,B,C;
var cnt=0;
phi1 = calRadian(p1.lat());
L1 = calRadian(p1.lng());
alpha1 = calRadian(deg1);
U1 = Math.atan((1 - f)*Math.tan(phi1));
sigma1 = Math.atan2(Math.tan(U1),Math.cos(alpha1));
sinalpha = Math.cos(U1)*Math.sin(alpha1);
cosalpha2 = 1-Math.pow(sinalpha,2);
u2 = cosalpha2*((Math.pow(a,2)-Math.pow(b,2))/Math.pow(b,2));
A = 1+u2/16384*(4096+u2*(-768+u2*(320-175*u2)));
B = u2/1024*(256+u2*(-128+u2*(74-47*u2)));
sigma = s/(b*A);
do {
sigmam = (2*sigma1+sigma)/2;
deltasigma = B*Math.sin(sigma)
*(Math.cos(2*sigmam)
+1/4*B*(Math.cos(sigma)*(-1+2*Math.pow(Math.cos(2*sigmam),2))
-1/6*B*Math.cos(2*sigmam)*(-3+4*Math.pow(Math.sin(sigma),2))*(-3+4*Math.pow(Math.cos(2*sigmam),2))));
formersigma = sigma;
sigma = s/(b*A)+deltasigma;
cnt++;
if (cnt>100) {
alert("error");
return null;
}
} while (Math.abs(sigma-formersigma)>epsilon);
phi2 = Math.atan2((Math.sin(U1)*Math.cos(sigma)+Math.cos(U1)*Math.sin(sigma)*Math.cos(alpha1))
,((1-f)*Math.sqrt(Math.pow(sinalpha,2)+Math.pow((Math.sin(U1)*Math.sin(sigma)-Math.cos(U1)*Math.cos(sigma)*Math.cos(alpha1)),2))));
lambda = Math.atan2(Math.sin(sigma)*Math.sin(alpha1),(Math.cos(U1)*Math.cos(sigma)-Math.sin(U1)*Math.sin(sigma)*Math.cos(alpha1)));
C = f/16*cosalpha2*(4+f*(4-3*cosalpha2));
L = lambda-(1-C)*f*sinalpha*(sigma+C*Math.sin(sigma)*(Math.cos(2*sigmam)+C*Math.cos(sigma)*(-1+2*Math.pow(Math.cos(2*sigmam),2))));
L2 = L + L1;
alpha2 = Math.atan(sinalpha/(-1*Math.sin(U1)*Math.sin(sigma)+Math.cos(U1)*Math.cos(sigma)*Math.cos(alpha1)));
p2 = new google.maps.LatLng(calDegree(phi2),calDegree(L2));
return p2;
}
各距離の円の色が重ならないよう、外側の円のポリゴンはドーナツ状にしています。Google Maps APIのポリゴンでは、内側領域のデータを外側とは逆回りの配列で追記することで中抜けの形を描画することが出来るので、外側の円の描画時にはその１つ内側のポリゴン配列を逆回りにしてPolygonクラスにデータを渡しています（上記ソース125～129行目）。
そして、上のHTMLを開いて描画させた結果が以下の図です。
計算・描画結果 ボツワナの首都（ハボロネ）から1000km,2000km,3000kmの同心円
同心円の色も重ならず、また半透明なので下の地図が見えなくなってしまうこともなくバッチリです。
「これで完成！」と思ったのですが・・・・・・。
距離を10000kmなど大きくして試してみると・・・
計算・描画結果 ボツワナの首都（ハボロネ）から3000km,5000km,10000kmの同心円
このように一番外側の距離円（もはや円ではありませんが）のポリゴンが本来とは逆の領域を選択してしまっています。その関係で内側の円には外側の円の色（今回の場合は赤）が重なってしまいました。
中心点や距離をいろいろ変えて描画させてみたところ、どうやらGoogle Maps APIのポリゴンは地球を１周して南北に２分割するような大きな領域を描く際には北半球側を選択するようです。
ところが悔しいことに、google.maps.Circleで同条件の同心円を描いた場合は以下の図のように南半球でもきちんと選択されています
同条件でにてgoogle.maps.Circleを利用した結果
「色が重なる程度の小さなことは気にせず、おとなしくgoogle.maps.Circleを使いなさい！」ということなんでしょうか？・・・いえいえ、まだ諦めません。次回へつづく
