[GMapsJ:2667] ランダムなポイントを並べ変えてきれいなポリゴンを作る方法について

[YON]

ランダムに複数のGLatLngオブジェクトを格納した配列pointsがあります。
その配列ををそのまま、new GPolygon(points)でポリゴン作成しても、
位置がランダムなため順番にならばず、きれいなポリゴンになりません。

そこで、以下のような配列の拡張処理を作り、一応思惑どおりになったのですが、もっと、楽な方法があるような気がします。

思いついた方法は、「まず配列データから中心点の座標を求め、中心点を
基点に各店をNE、SE、SW、NWの4つの象限の配列に分類格納しなお
し、それぞれを東側は緯度の昇順、西側は緯度の降順にソートしたうえで
再び、NE、SE、SW、NWの順でコンカチして、新しい配列に格納する。」
といった手法です。
（南半球、西半球にも対応させるため90と180を先に足して、後で引いてます）
何か、お気づきの点やアドバイスがありましたらよろしくお願いします。

（配列の拡張処理）
Array.prototype.polygonsort=function(){
var wrk=[];var ne=[];var se=[];var sw=[];var nw=[];var result=[];
var sum={"lat":0,"lng":0};var ave={"lat":0,"lng":0};
for(var i=0;i<this.length;i++){wrk.push({"lat":this[i].lat()
+90,"lng":this[i].lng()+180});}
for(var i=0;i<this.length;i++){sum.lat=sum.lat
+wrk[i].lat;sum.lng=sum.lng+wrk[i].lng;}
ave.lat=sum.lat/this.length;ave.lng=sum.lng/this.length;
for(var i=0;i<this.length;i++){
if(wrk[i].lat>ave.lat){
if(wrk[i].lng>ave.lng){ne.push(wrk[i]);}else{nw.push(wrk[i]);}
}else{
if(wrk[i].lng>ave.lng){se.push(wrk[i]);}else{sw.push(wrk[i]);}
}
}
ne.sort(function(a,b){return(b.lat - a.lat)});
se.sort(function(a,b){return(b.lat - a.lat)});
sw.sort(function(a,b){return(a.lat - b.lat)});
nw.sort(function(a,b){return(a.lat - b.lat)});
var all=ne.concat(se,sw,nw);
for(var i=0;i<all.length;i++){result.push(new
GLatLng(all[i].lat-90,all[i].lng-180));}
result.push(result[0]);
result.center=new GLatLng(ave.lat-90,ave.lng-180);
return result;
}

↓のように実装
var map;
function load(){
if (GBrowserIsCompatible()) {
map = new GMap2(document.getElementById("map"),{size:new
GSize(800,600)});
map.setCenter(new GLatLng(35.6851,139.6644), 11);
map.addControl(new GLargeMapControl());
// AddMarker();
}
}
function AddMarker(){
var points=[];
var data=[{"lng":139.6138,"lat":35.6907},{"lng":139.9193,"lat":
35.7251},
{"lng":139.8183,"lat":35.7021},{"lng":139.5104,"lat":
35.7147},
{"lng":139.9223,"lat":35.5805},{"lng":139.4064,"lat":
35.7896}];
for (var i = 0; i < data.length; i++) {
var marker= new GMarker(new GLatLng(data[i].lat,data[i].lng));
points.push(marker.getLatLng());
map.addOverlay(marker);
}
var polygonpoints=points.polygonsort();
var polygon =new GPolygon(polygonpoints);
map.addOverlay(polygon);
map.panTo(polygonpoints.center);
var polygonBounds = polygon.getBounds();
map.setZoom(map.getBoundsZoomLevel(polygonBounds));
}

--
このメールは Google グループのグループ「Google-Maps-API-Japan」の登録者に送られています。
このグループに投稿するには、google-map...@googlegroups.com にメールを送信してください。
このグループから退会するには、google-maps-api-...@googlegroups.com にメールを送信してください。
詳細については、http://groups.google.com/group/google-maps-api-japan?hl=ja からこのグループにアクセスしてください。

[W少年]

難しすぎてよく解りません、と言いながらコメントしていますが私も似たようなことを考えたことがあります。
生涯学習で位置天文学のサイトをつくっているので、コーディングはお手上げですがこの辺の理論は少し強いと思います。
富士山頂のDEMデータから尾根線上のポイントを抽出して最後に円形につないで山頂の形状をグラフィックしようと考えてDEMデータから山頂付近をタイ
ル状に並べて検討しました。
結果は小さな丘陵が沢山あって環状に近似することができなくてDEMデータを全て計算することになったのですが、抽出の仕上げとして以下の方法を考えて
いました。
まず、全データの中心をどこにするかが物議をかもす点ではありますが厳密には各点の座標を緯度、経度から三次元座標に変換して、その平均の座標から緯
度、経度を計算すると厳密にポリゴンの中心を求めたことになります。
e^2=1-(dp/de)^2
X=ρ*cos(φ')=(de/sqrt(1-e^2*sin(φ)^2)+h)*cos(φ)
Z=ρ*sin(φ')=(de*(1-e^2)/sqrt(1-e^2*sin(φ)^2)+h)*sin(φ)
ポイントの座標は
(X*cos(λ),X*sin(λ),Z)
座標から緯度、経度を求めるには上記座標から抽出したλ、φ'から
φ=atan((de/dp)^2*tan(φ'))
により中心点の緯度、経度を求めます。
その後、全ての点の座標から中心点の座標の値を引いて中心点の経度λで-λ回転したのち中心点の緯度φで90°-φ回転します。
どういう回転演算かは誤植があるとバグの原因になるので考えてみてください。
すると、X座標の符号を逆にしたのちXの正の方向、つまりX座標が負の値の方向からY座標が正の方向へまわるように角度をとるとそれが中心点から各点へ
の方角になります。
これをソートの条件式にして方角でソートするときれいな星型ができるはずです。

[W少年]

ですが、この計算をまともにやるとかなり大変なので、以下の近似法が有効でしょう。
平面型のグーグルマップはメルカトル図法なので緯度と経度の線は長方形をなしています。
よって単純に緯度、経度のそれぞれの平均値を中心としてそこからの方角を緯度、経度の中心との差から求めると回転演算せずに二次元の計算だけで方角が
求
められます。
それにより求めた方角でソートすると大体の場合上記の方法とあまり変わらない結果となるはずです。

[YON]

W少年様、ありがとうございます。
中心からの方位を求めて、それでソートするというやり方がありましたね。ぱっと思いつきません
でした、感謝します。とりあえず、高度とかゆがみは無視して地球を滑らかな真球とみなせば、
（緯度、軽度）を極座標になおして球面三角法で計算するのですよね。
球面三角法で、2点間の距離方位を取得するコードは、以前↓を参考に作った事があるので、
理解できます。

http://www-seis1.kugi.kyoto-u.ac.jp/visual/grams/index4.html
http://gpscycling.net/fland/make/ckyoho.htm

W少年様の中心点を算出する数式も参考にさせてもらってさっそく試してみます。

　ただ、精度はあがりそうですが、javascript的には三角関数、逆三角関数、平方根、除算を
多用する事になるので、計算誤差と処理速度を考えると、単純簡易方法でもよいような
気もします。ケースバイケースですね。

[kino]

シンプルで読みやすいコードですね。

凸包問題を解きたいのかと思いましたが、コードを見ると凹多角形ですか・・・
凹多角形は一意に定まらないので、「意見はありますか」という質問には「どうしたいんですか？」
としか返せないと思います。

例えば、提案方式では「☆」や「凹」の形は再現できません。
そこを「どのように再現したいか」を定めずにプログラムテクニックの詳細に踏み込むと
単に話が発散してしまうと思います。

凸包や凹多角形の話は単純そうに見えて奥が深いので、Google Mapsで議論するよりは
数学関連のところの方がいいのでは？　という気がします。

凸包だと一意に定まるので話は単純です。
http://www.geocodezip.com/map-markers_ConvexHull_Polygon.asp

[W少年]

ひとまず、先の投稿では地球楕円体面上で中心を求める方法を真面目に解説しましたが、
極地を含む場合とか緯度経度の分布範囲が広い場合はまともに三次元演算が必要かもしれませんが、
市街地図の中などの局所的な範囲でやる場合はそこまで大がかりな計算をする必要はないと思います。
上のリンクは真面目に計算しているように見えますが、下のリンクは局所的な近似が入っています。
私もダイヤモンド富士の計算を始めたころは球面三角法で計算していましたが、
太陽や月の出没計算以外では無駄に複雑化するだけであまりメリットがないので、
位置天文学的な三次元回転を多用する方法に変えたら保守効率が格段に上がりました。
簡易計算法をもう少し解説するとポイント(φ1、λ1)、(φ2、λ2)、(φ3、λ3)・・・
に対して中心点を
(Σφ/n、Σλ/n)
方角を
azimuth=Math.acos((φk-Σφ/n)/Math.sqrt((φk-Σφ/n)*(φk-Σφ/n)+(λk-Σλ/
n)*(λk-Σλ/n)))
if (λk-Σλ/n<0) azimuth=2*Math.PI-azimuth
ただし、(Σφ/n、Σλ/n)は緯度、経度の各々の相和をデータ数で割った単純平均の値で、経度が180度をまたがる場合は+か-のどちらかにそろえ
てから行わないと結果が出鱈目になります。
ここで求められる方角はソートの引数としての方角で厳密な方角の値ではありません。
なお、現在使われている赤道半径deと極半径dpの値は以下の値です。
de=6378137m
dp=6356752m

[YON]

ありがとうございます。
凸包や凹多角形ですか、確かに凹多角形は一意に定まらないですね。見落としました。
確かにポリゴン問題は奥が深いです。
中心点座標と、作ったポリゴンに内接する円の中心座標の関係とか、
中心点から各頂点に線を引いた時の交点の数とかで、何とか、一意に選ぶ条件を
考えられそうな気がしますが、難しそうですね。数値幾何学の研究が目的でないので、
断念し、適当な方法で妥協します。

[YON]

中心点から等距離にある地点の緯度・経度を球面三角法（dp=deとして）で360個（0度～359度）計算し
その点からW少年様ご提示の簡易計算法
「azimuth=Math.acos((φk-Σφ/n)/Math.sqrt((φk-Σφ/n)*(φk-Σφ/n)+(λk-Σλ/
n)*(λk-Σλ/n)))
if (λk-Σλ/n<0) azimuth=2*Math.PI-azimuth 」
を使って、中心点からの方位を再計算し、さらに各点（緯度・経度）から中心点までの方位と距離を
逆計算してみました。

そして結果の方位差・距離差ををエクセルで、0～360度にわたってプロットしてみると
方位差は見事な4/ΠずれたΠ周期の正弦曲線、
距離差はΠ周期の余弦曲線になりました。

なんでそうなるかは、さっぱり解りませんが、数学は不思議で美しいです。
理由については、どの方位を0度に選んだかによる影響のような気がしますが、．．．
Google Mapsで議論すべき事じゃないですよね。

[W少年]

> 距離差はΠ周期の余弦曲線になりました。
それは多分計算の過程で緯度と経度が入れ替わっているので直したやつを送りました。
なお、簡略式でも緯度方向と経度方向に各々の曲率半径をかけると球面上の実際の方角に近くなります。
球であれば経度の差に緯度の余弦をかければいいのですが、実際の地球の形に即した回転楕円体で計算しようとするともう少し複雑になります。