面積計算ツール別の誤差値について

[hachi]

はじめましてhachiと申します。

登録の件ではお騒がせしてしまい申し訳ございませんでした。

以下の事柄について、ご存じでしたら教えていただけますと幸いです。

使用しているQGISは3.10で、座標系はプロジェクト・地物ともにEPSG2447を指定しています。

フィールド演算で面積計算（area($geometry))した場合と

「ジオメトリ属性の追加」で面積計算した場合において、

小数点第6位前後で誤差が生じます。

10個ほど大きさの違う地物を作成し再計算した場合も同様です。

使用したshpを添付させていただきます。

なお、念のためCADで確認したところ、同様に小数点第6位前後で誤差が確認されました。

この際の誤差は上記の計算で出した値とはどれとも違う値となります。

フィールド演算機と「ジオメトリ属性の追加」で使用している求積法が違う、もしくは浮動小数点の誤差のどちらかではと思っていたのですが、QGISの設定や面積計算式に過不足があり、それらが影響した結果、誤差が生じているという可能性はあるでしょうか。

そうではなかった場合、どちらのツールを使った方がより正確な値なのでしょうか。

㎡単位の小数点第6位で生じる誤差が、地上で行うほとんどの計測業務に影響を及ぼすものではないと承知の上でお聞きしております。

何卒宜しくお願いいたします。

[ありた]

こんにちは

まず、 area($geometry) と「ジオメトリ属性の追加」の差ですが、2進法でデータを

管理しているコンピュータ上で実数を扱う場合、その有効数字は16桁程度です。

整数部で10桁近い数値の場合、小数部の有効桁は6桁程度しかありません。

加えて、2進数で数値を扱うため、そもそも数値を正確に表現できないこともあります。

ですので、今回の場合は実数の比較上ほぼ同値とみなしてよいと思います。

（基本的な算出方法は同一でも、加減算や乗算の順番によって誤差が生じて、同一の

結果が得られないことがあります。）

次に ArcGIS の結果ですが、なんとなく座標値もしくは計算結果に対して丸め処理を

行っているように思います。算出の面積値がキレイな数値すぎます。

[hachi]

ありた様

夜分に申し訳ございません。ご確認いただきありがとうございます。

ArcGISでの計算結果については仰る通りで、桁が多いものは自動で丸め処理してしまうようです。

フィールド演算を使えば指定が効いたのかもしれません。

この、添付画像のid　1･2･3の行のように、比較的整数値の桁数が少ない場合に、arc_mと q_area（どちらも㎡単位）の両方が小数部9桁程度まで近似値を取ったことに対し、q_pro_mは明らかに小数部6桁目以降の数値が変わっています。

有効桁数の範囲内かと思うのですが、これは単純に、「ジオメトリ属性値の追加」での算出方法（加減算や乗算の順番など）がArcGISやフィールド演算での方法とは異なっている可能性がある、と思えば良いでしょうか。

また、算出方法が違っていると思われた場合、設定等で統一させることは可能でしょうか。

2021年9月16日木曜日 0:49:12 UTC+9 ありた:

[yoichi...@gmail.com]

かやまと申します

QGISの面積計算の関数は2種類あります

$area   と area  です

https://docs.qgis.org/3.16/ja/docs/user_manual/working_with_vector/functions_list.html#area

$area  は

The area calculated by this function respects both the current project's ellipsoid setting and area unit settings

 area  のほうは

Returns the area of a geometry polygon object. Calculations are always planimetric in the Spatial Reference System (SRS) of this geometry, and the units of the returned area will match the units for the SRS.

とあって

This differs from the calculations performed by the $area function, which will perform ellipsoidal calculations based on the project's ellipsoid and area unit settings.

とあります

楕円体を考慮した面積計算なのかどうかの違いがあります。

ご指定の操作でどちらの関数が使われているかはわかりませんが、この2種類の面積計算の違いの可能性もあります。

2021年9月16日木曜日 3:28:43 UTC+9 hachi:

[ありた]

たしかに ID: 1, 2, 3 も小数点以下6～7桁目程度の差が生じていますね。

たとえば ID: 2 のジオメトリ情報（座標値）を確認してみると以下のとおりでした。

MultiPolygon (((

   53727.09608578203915386 -103390.91826425890030805,

   53822.01992106921534287 -103359.27698582984157838,

   53822.01992106921534287 -103472.28155164791678544,

   53736.13645104748138692 -103476.80173428064153995,

   53727.09608578203915386 -103390.91826425890030805

)))

ファイルは Shapefile 形式であり、仕様上 double （倍精度浮動小数点数）として

保存されているため、有効な桁数は約16桁です（2進数で53ビット分）。

（表記として16桁以上表示されていますが、有効桁としては約16桁）

https://www.esrij.com/cgi-bin/wp/wp-content/uploads/documents/shapefile_j.pdf

平面上の多角形の面積の求め方でよく使われる式は  となります。

このため乗算の結果で整数部が11桁生じ、小数部の精度は5桁となります。

これは (53800, -103400) 付近に存在することで発生するものですので、あらかじめ

差し引いて (0, 0) 近くに並行移動させ、乗算を行っても整数部の桁数が増えないようにし

できるだけ有効桁数を保つコーディングテクニックなどもあります。

こういった実装の差や、場合によっては処理系の最適化で結果が変化する余地があります。

ちなみに、検算のために浮動小数点数ではない手法（有理数型）で ID: 2 の面積を

誤差なく計算したところ 9071.830658727436 になりました。

他のポリゴンについては確認していないのと、 QGIS の実装そのものを見れてないので

必ずしもそうだとは言い切れませんが、この例だと area($geometry) がもっとも

正確な結果を示しました。

なお、かやまさんがおっしゃられている地図平面上の面積 area($geometry) と

地球楕円体上の $area は異なる結果を示し、 GIS ではこちらの差異の方が、よく

トラブルや勘違いのもとになることがあります。

（今回は処理結果を確認すると地図平面上の面積における誤差のようでした）

[hachi]

かやま様　ありた様

ありがとうございます。

御礼が大変遅くなり申し訳ございません。

かやま様へ

回答ありがとうございます。

既にありた様の方で補足いただいているのですが、

今回は地図平面上の面積計算ツール同士で生じた微細な誤差の理由について質問させていただいておりました。

説明が不足しており誠に申し訳ございません。

ありた様へ

解説、検算までしていただきありがとうございます。

＞乗算の結果で整数部が11桁生じ、小数部の精度は5桁となる

この解釈として、計算過程において整数部が増加すると、整数部に引っ張られて有効桁数が繰り上がるので、最終の計算結果の有効桁数も繰り上がった桁数となる、という理解で合っていますでしょうか。

考え方）

①53727.〇〇〇…　小数部11桁まで有効

②103390.〇〇〇…　小数部10桁まで有効

①×②＝5554834530.〇〇〇…　小数部5桁まで有効

…

最終＝9071.830658727436　　小数部5桁まで有効

（※途中計算で有効桁数が小数部5桁を維持したと想定）

ArcGISジオメトリ演算とQGISフィールド演算結果、有理数型で検算いただいた計算結果が有効桁数以上でも近似したことを考えると、同じ座標値を用いたであろうプロセシングツールでの計算結果が有効桁数以上で異なる原因については、一つ前の返信で仰っていただいたように算出方法や計算順序が微妙に違っていたのかなと思います。

今回、誤差が検査基準桁数まで拡大することはなさそうだと考えられたため、プロセシングツールで使用されている算出方法までは追わなくて良いと思っておりますが、ツールを使用される中で違和感など感じられたことはあったでしょうか。

2021年9月17日金曜日 0:05:04 UTC+9 ありた:

[ありた]

考え方としてはあってます。

浮動小数点数として、コンピュータ内部では（正確には2進数としてですが）

5.372709608578204 × 10^4　　×　　-1.033592769858298 × 10^5　＝

　-5.553193805974640 × 10^9　という計算が行われています。

確かに、今回のご質問で複数の処理間で計算結果が異なることは少し盲点ではありました。

ただ計算機を使う以上、誤差や桁落ちなどは「あるある」なことです。

また地理学としても、地図の精度として常に誤差は身近にあるもので、たとえば

縮尺1/1000相当の地図情報レベル1000では精度（標準偏差）0.7m以内と定められています。

そのなかで、計算機上で有効桁内にあったとしても各点の座標値としてマイクロメートル

レベルの精度があるのか、その座標値をもとにした面積が真値なのか。ということは

目的次第ですが考えたほうがいいのかもしれません。

[hachi]

ありた様

度々返信遅れまして申し訳ございません。

ご回答ありがとうございます。

仰られるとおり、計算機の限界による誤差については発生せざるを得ないものかと思います。

当方で取り扱う数値上でも、今回の誤差は精度的に全く問題ありません。

今回質問させていただいた発端は、ツールによって生まれる誤差が、＄areaのように

そもそも基準とする座標系が違っているような、大きな誤差を生むものであったら困る、といった理由からでした。

フィールド演算とプロセシングツールの演算で起こる誤差が、計算機上必然的に起こる誤差と言えるのであれば気にしなくても良さそうです。

CADでは複数の求積計算が可能ですが、同じ座標値を使用しても、求積法ごとにマイクロ単位で違う値を算出します。

今回の件はそれと同様の事象であると認識いたしました。

QGISは沢山の方のアイデアが組み込まれて作られているので、ツールによってアルゴリズムが違うことも当然あり得ることでしたよね。

これまではGISで取り扱うような計算は様式が統一されたものと思い込んでおりましたので、今回の件で新たに認識できてよかったです。

誤差発生に対して丁寧に回答下さり本当にありがとうございました。

2021年9月19日日曜日 9:12:53 UTC+9 ありた: