トップ地図システム基礎技術 > 極座標で表された多角形の面積計算

極座標で表された多角形の面積計算

openlayersにおける面積計算

openlayers[1]では 文献[2]に基づいて、次のプログラムで面積を求めている。

/**
 * Returns the spherical area for a list of coordinates.
 *
 * [Reference](https://trs-new.jpl.nasa.gov/handle/2014/40409)
 * Robert. G. Chamberlain and William H. Duquette, "Some Algorithms for
 * Polygons on a Sphere", JPL Publication 07-03, Jet Propulsion
 * Laboratory, Pasadena, CA, June 2007
 *
 * @param {Array} coordinates List of coordinates of a linear
 * ring. If the ring is oriented clockwise, the area will be positive,
 * otherwise it will be negative.
 * @param {number} radius The sphere radius.
 * @return {number} Area (in square meters).
 */
function getAreaInternal(coordinates, radius) {
  let area = 0;
  const len = coordinates.length;
  let x1 = coordinates[len - 1][0];
  let y1 = coordinates[len - 1][1];
  for (let i = 0; i < len; i++) {
    const x2 = coordinates[i][0];
    const y2 = coordinates[i][1];
    area +=
      toRadians(x2 - x1) *
      (2 + Math.sin(toRadians(y1)) + Math.sin(toRadians(y2)));
    x1 = x2;
    y1 = y2;
  }
  return (area * radius * radius) / 2.0;
}

多角形の面積計算

Google Map の緯度からY座標への変換は複雑なため、XY平面座標の矩形であっても、 面積は 縦 x 横 で求まるわけではない。

地球は球体(より正確には楕円体)であるから、表面上の多角形の面積は極座標の方が求めやすい。 厳密な計算は複雑であるが、通常の地図ではそれほど精密な値はいらない。 その場合には比較的簡単なプログラムで地球上の多角形の面積を求めることができる。

ここで、多角形は、森林、公園、建物などを表す。地球の全表面に比べればごく狭いものとする。

下のプログラムで、楕円高は丁度 100 というわけではない。 どのみち、長半径に比べればごく小さい値なので、この数値は無視してもよい。

座標値は 整数化した経度、緯度の順で1次元配列で与えている。 引数を calcWayArea(double[] lons, double[] lats) とすればプログラムはもっと簡単になる。

メソッドが戻す面積の単位は ㎡ である。

  double calcWayArea(int[] lonlats) {
      double Radius = 6378137 + 100;      // 長半径 + 楕円高(仮)
      double[] xs = new double[lonlats.length/2];
      double[] ys = new double[lonlats.length/2];
      for (int n = 0; n < lonlats.length/2; n++) {
          xs[n] = lonlats[n*2]/10000000.0;
          ys[n] = lonlats[n*2+1]/10000000.0;
      }
      double wayarea = 0;
      double x1 = xs[xs.length-2];
      double y1 = ys[ys.length-2];
      for (int n = 0; n < xs.length-1; n++) {
          double x2 = xs[n];
          double y2 = ys[n];
          wayarea += toRadian(x2-x1) * 
                     (2+Math.sin(toRadian(y1))+Math.sin(toRadian(y2)));
          x1 = x2;
          y1 = y2;
      }
      return wayarea * Radius * Radius / 2.0;
  }   // 反時計回りが正、時計回りは負となる

この計算結果は反時計回りが正、時計回りは負となる。地図システムではこの情報も重要であり、 これを使って、多角形のデータの並びを一方向に統一する。

文献[2]では時計回りが正となっている。因みに OSMの海岸線は殆どが反時計回りであり、面積が負になる。

参考記事

[1] openlayers/src/ol/sphere.js
[2] Some Algorithms for Polygons on a Sphere