単純な閉曲線上の点を表す緯度と経度のペアの任意のセットがあるとします。デカルト空間では、グリーンの定理を使用して、このような曲線で囲まれた面積を簡単に計算できます。球の表面の面積を計算するための類似のアプローチは何ですか？私が求めているのは、Matlabのareaint関数の背後にあるアルゴリズム（の近似値でさえ）だと思います。
28091 次
6 に答える
13
これを行うにはいくつかの方法があります。
1）緯度帯からの寄与を統合します。ここで、各ストリップの面積は（Rcos（A）（B1-B0））（RdA）になります。ここで、Aは緯度、B1とB0は開始経度と終了経度、すべての角度はラジアンです。
2）表面を球面三角形に分割し、ジラールの定理を使用して面積を計算し、これらを合計します。
3）ここでJames Schekが提案したように、GIS作業では、平坦な空間への投影を維持する面積を使用して、そこでの面積を計算します。
データの説明から、最初の方法が最も簡単なように聞こえます。（もちろん、私が知らない他のもっと簡単な方法があるかもしれません。）
編集–これら2つの方法を比較します。
最初の検査では、球面三角形アプローチが最も簡単に見えるかもしれませんが、一般的にはそうではありません。問題は、領域を三角形に分割するだけでなく、球面三角形、つまり、辺が大円弧である三角形に分割する必要があることです。たとえば、緯度の境界は修飾されないため、これらの境界は、大円の弧により近いエッジに分割する必要があります。そして、これは、大円が球面角度の特定の組み合わせを必要とする任意のエッジに対して行うのがより困難になります。たとえば、球の周りの中央のバンドをどのように分割するかを考えてみましょう。たとえば、緯度0から45度までのすべての領域を球面三角形に分割します。
結局、各方法で同様のエラーを使用してこれを適切に行う場合、方法2では三角形の数は少なくなりますが、決定するのは困難になります。方法1はより多くのストリップを提供しますが、それらを決定するのは簡単です。したがって、より良いアプローチとして方法1をお勧めします。
于 2009-08-27T17:51:49.293 に答える
10
MATLABの「areaint」関数をJavaで書き直しましたが、これはまったく同じ結果になります。「areaint」は「単位あたりの表面」を計算するので、答えに地球の表面積（5.10072e14平方メートル）を掛けました。
private double area(ArrayList<Double> lats,ArrayList<Double> lons)
{
double sum=0;
double prevcolat=0;
double prevaz=0;
double colat0=0;
double az0=0;
for (int i=0;i<lats.size();i++)
{
double colat=2*Math.atan2(Math.sqrt(Math.pow(Math.sin(lats.get(i)*Math.PI/180/2), 2)+ Math.cos(lats.get(i)*Math.PI/180)*Math.pow(Math.sin(lons.get(i)*Math.PI/180/2), 2)),Math.sqrt(1- Math.pow(Math.sin(lats.get(i)*Math.PI/180/2), 2)- Math.cos(lats.get(i)*Math.PI/180)*Math.pow(Math.sin(lons.get(i)*Math.PI/180/2), 2)));
double az=0;
if (lats.get(i)>=90)
{
az=0;
}
else if (lats.get(i)<=-90)
{
az=Math.PI;
}
else
{
az=Math.atan2(Math.cos(lats.get(i)*Math.PI/180) * Math.sin(lons.get(i)*Math.PI/180),Math.sin(lats.get(i)*Math.PI/180))% (2*Math.PI);
}
if(i==0)
{
colat0=colat;
az0=az;
}
if(i>0 && i<lats.size())
{
sum=sum+(1-Math.cos(prevcolat + (colat-prevcolat)/2))*Math.PI*((Math.abs(az-prevaz)/Math.PI)-2*Math.ceil(((Math.abs(az-prevaz)/Math.PI)-1)/2))* Math.signum(az-prevaz);
}
prevcolat=colat;
prevaz=az;
}
sum=sum+(1-Math.cos(prevcolat + (colat0-prevcolat)/2))*(az0-prevaz);
return 5.10072E14* Math.min(Math.abs(sum)/4/Math.PI,1-Math.abs(sum)/4/Math.PI);
}
于 2015-02-10T05:21:17.777 に答える
9
タグの1つで「地理」について言及しているので、ジオイドの表面上のポリゴンの領域の後にいるとしか考えられません。通常、これは地理座標系（つまり、lon / lat）ではなく投影座標系を使用して行われます。lon / latで行う場合、返される測定単位は球の表面のパーセントになると思います。
より「GIS」フレーバーでこれを実行したい場合は、領域の測定単位を選択し、領域を保持する適切な投影法を見つける必要があります（すべてがそうであるとは限りません）。任意のポリゴンの計算について話しているので、ランベルト正積図法のようなものを使用します。投影の原点/中心をポリゴンの中心に設定し、ポリゴンを新しい座標系に投影してから、標準の平面手法を使用して面積を計算します。
地理的領域で多くのポリゴンを作成する必要がある場合は、他の投影法が機能する可能性があります（または十分に近くなります）。たとえば、UTMは、すべてのポリゴンが単一の子午線の周りにクラスター化されている場合の優れた近似です。
これがMatlabのareaint関数の動作と関係があるかどうかはわかりません。
于 2009-08-28T15:42:35.570 に答える
5
Matlabの機能については何も知りませんが、ここで説明します。球形のポリゴンを球形の三角形に分割することを検討してください。たとえば、頂点から対角線を描画します。球面三角形の表面積は次の式で与えられます。
R^2 * ( A + B + C - \pi)
ここRで、は球の半径、、、AおよびBはC三角形の内角（ラジアン）です。括弧内の量は「球形過剰」として知られています。
側面のnポリゴンはn-2三角形に分割されます。すべての三角形を合計し、の共通因子を抽出し、R^2すべてを\piまとめると、ポリゴンの面積は次のようになります。
R^2 * ( S - (n-2)\pi )
ここSで、はポリゴンの角度の合計です。括弧内の量は、ポリゴンの球形の超過分です。
[編集]これは、ポリゴンが凸面であるかどうかに関係なく当てはまります。重要なのは、三角形に分割できることです。
あなたは少しのベクトル数学から角度を決定することができます。3つの頂点、、がAありB、Cの角度に関心があるとしBます。したがって、B大円セグメント（ポリゴンのエッジ）に沿った点から球への2つの接線ベクトル（それらの大きさは無関係）を見つける必要があります。のためにそれを解決しましょうBA。大円は、OAとOBで定義される平面にあります。ここOで、は球の中心であるため、法線ベクトルに垂直である必要がありますOA x OB。OBまた、そこに接しているので、垂直である必要があります。したがって、そのようなベクトルはによって与えられOB x (OA x OB)ます。右手の法則を使用して、これが適切な方向にあることを確認できます。これはに簡略化されることにも注意してくださいOA * (OB.OB) - OB * (OB.OA) = OA * |OB| - OB * (OB.OA)。
次に、古き良き内積を使用して、辺間の角度を見つけることができます。BA'.BC' = |BA'|*|BC'|*cos(B)ここで、BA'およびは、辺に沿ったからおよびへBC'の接線ベクトルです。BAC
[これらは接線ベクトルであり、点間の文字通りではないことを明確にするために編集されました]
于 2009-08-27T17:51:25.050 に答える
0
上記の回答に大まかに触発されたPython3の実装を次に示します。
def polygon_area(lats, lons, algorithm = 0, radius = 6378137):
"""
Computes area of spherical polygon, assuming spherical Earth.
Returns result in ratio of the sphere's area if the radius is specified.
Otherwise, in the units of provided radius.
lats and lons are in degrees.
"""
from numpy import arctan2, cos, sin, sqrt, pi, power, append, diff, deg2rad
lats = np.deg2rad(lats)
lons = np.deg2rad(lons)
# Line integral based on Green's Theorem, assumes spherical Earth
#close polygon
if lats[0]!=lats[-1]:
lats = append(lats, lats[0])
lons = append(lons, lons[0])
#colatitudes relative to (0,0)
a = sin(lats/2)**2 + cos(lats)* sin(lons/2)**2
colat = 2*arctan2( sqrt(a), sqrt(1-a) )
#azimuths relative to (0,0)
az = arctan2(cos(lats) * sin(lons), sin(lats)) % (2*pi)
# Calculate diffs
# daz = diff(az) % (2*pi)
daz = diff(az)
daz = (daz + pi) % (2 * pi) - pi
deltas=diff(colat)/2
colat=colat[0:-1]+deltas
# Perform integral
integrands = (1-cos(colat)) * daz
# Integrate
area = abs(sum(integrands))/(4*pi)
area = min(area,1-area)
if radius is not None: #return in units of radius
return area * 4*pi*radius**2
else: #return in ratio of sphere total area
return area
ここで、もう少し明示的なバージョン（および、より多くの参照とTODOを含む...）を見つけてください。
于 2020-04-12T19:11:41.043 に答える