float32で緯度経度を表したときの誤差はどれくらいか(float64だとメモリがもったいないの巻)

# はじめに

数ギガバイト単位の float64 の 緯度経度データをオンメモリで持つ羽目になり、メモリ削減のため float32 が使えないかと考えた際のメモです。

## IEEE754

コンピュータ上での数値を扱う際、その値は最終的には何らかの形で2進数で表す必要があります。

小数を表す方法には
 - 固定小数点(整数を1000など決まった値で割るなどして小数にする)
 - 固定小数10進表記(16進数1桁に10進1桁を割り当てる)
などいくつかありますが、演算効率や各種の言語でのサポートの厚さからよく使われているのは浮動小数点（小数点の位置が決まっていない）とよばれる形式です。

この浮動小数点の2進数エンコード方法について世界標準として使われているのが  [IEEE754](https://ja.wikipedia.org/wiki/IEEE_754) という規格です。

IEEE754 では 16bit、32bit、64bit などでの浮動小数点数の表現方法が規定されています。下記にその一部を示します。

| ビット数 | 名称 | 日本語名称 |
|---------|------|-----------|
| 16      |binary16|半精度浮動小数点数|
| 32      |binary32|単精度浮動小数点数|
| 64      |binary64|倍精度浮動小数点数|

## float32 の精度

さて、本件は float32 でどれくらいの緯度経度の精度が出るのかを調べたいということでした。

float32 は 32bit、IEEE754 での単精度浮動小数点数です。

32bit の内訳は、符号 s が 1ビット、指数部 e が 8ビット、仮数部 f が 23ビット、全部で 32ビット あります。

|31|30|29|28|27|26|25|24|23|22|21|20|19|18|17|16|15|14|13|12|11|10|09|08|07|06|05|04|03|02|01|00|
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|s|e7|e6|e5|e4|e3|e2|e1|e0|f22|f21|f20|f19|f18|f17|f16|f15|f14|f13|f12|f11|f10|f9|f8|f7|f6|f5|f4|f3|f2|f1|f0|

それぞれの意味を見ていきます

### 符号 s
0:正、1:負 の意味です。

### 仮数部 f22-f0
ざっくりいうと有効ケタ数に相当する部分の値を2進数で表したものになります。1で始まるようにしてやり、その1自体は自明なのでエンコードしません。ということで 仮数部の最初の桁が1だと 0.5、次の桁が 0.25、その次が 0.125 ... というように 23ビットまで表現できます。

f22-f0 の各桁に対応する10進小数の値を表にしました。あわせて、その桁までで何bitあるのかも示しています。[^1]

[^1]: ~~JavaScript~~ BigDecimal は 10進小数をそのまま扱ってくれるのでこういう表を作るときに誤差を気にしなくてよいのが便利です。

```
~% kotlinc
Welcome to Kotlin version 1.5.31 (JRE 14.0.1+7)
Type :help for help, :quit for quit
>>> var v = java.math.BigDecimal("0.5"); var d2 = java.math.BigDecimal(2);  for (i in 0..22) { println("f${22-i}: ${i+1}bit: ${v.toPlainString()}"); v = v.divide(d2) }
f22: 1bit: 0.5
f21: 2bit: 0.25
f20: 3bit: 0.125
f19: 4bit: 0.0625
f18: 5bit: 0.03125
f17: 6bit: 0.015625
f16: 7bit: 0.0078125
f15: 8bit: 0.00390625
f14: 9bit: 0.001953125
f13: 10bit: 0.0009765625
f12: 11bit: 0.00048828125
f11: 12bit: 0.000244140625
f10: 13bit: 0.0001220703125
f9: 14bit: 0.00006103515625
f8: 15bit: 0.000030517578125
f7: 16bit: 0.0000152587890625
f6: 17bit: 0.00000762939453125
f5: 18bit: 0.000003814697265625
f4: 19bit: 0.0000019073486328125
f3: 20bit: 0.00000095367431640625
f2: 21bit: 0.000000476837158203125
f1: 22bit: 0.0000002384185791015625
f0: 23bit: 0.00000011920928955078125
```

ということで 23bit の場合、10進表記だとおおよそ 0.0000001刻み、$10^{-6}$ ～ $10^{-7}$ くらいの精度があることがわかりました。

### 指数部 e7-e0 

指数部とは、仮数部を何ビットシフトさせるか（ざっくりいうと2倍ないし2で割ることを何回繰り返すか) を示したもので、8ビットあります。ということは 0 ～ 255 までシフト幅を表現できるのですが、127 のオフセットがついているので -128 ～ 127 までのシフト幅ということになります。
（浮動小数点数の値としては $2^{-128}$ ～ $2^{127}$ の範囲を表せる）

### 特別な値

不明な数 `NaN` や 正負の無限大 `Inf` `-Inf` を表すために特別な値があります

| s | e7-e0 | f22-f0 | 特別な値 |
|---|-------|--------|---------|
| 1 | 全部1 | 全部0 | +Inf |
| 0 | 全部1 | 全部0 | -Inf |
| 1 | どれかが1 | 全部0 | +0 |
| 0 | どれかが1 | 全部0 | -0 |
| - | 全部1 | どれかが1 | NaN |
| - | 全部0 | どれかが1 | 非正規化数 - 仮数部の最上位に1がないものとして そのまま扱う |

## 緯度経度

ようやく本題となる、上記を緯度経度に適用した場合の精度について考えます。

経度は、日本では通常 120から150あたりまでの範囲をとり、具体的には 136.398375383923 のような10進小数で表されます。緯度は同様に 20から40あたりの範囲をとります。

float32 で 整数部分の 130 前後の数値を表そうとすると `1000 0010` と 2進数で 8ケタほど必要です。仮数部は 23ビットだったので、先頭の1を除いて7ビットを緯度経度の整数部分に取られてしまい、小数部分の表現に使える残りは16ビットになります。

16ビットで経度の小数部分を表した場合の精度を考えてみます。

上の表から 16bit の場合の 10進での精度は 0.000015.. でした。

この精度が実際の距離としてどれくらいかを概算します。

赤道の周囲が約4万km なので、赤道上ではざっくり1度 ≒ 111km と考えられます。日本に相当する緯度35度上では、111km x cos35° ≒ 約91km。

ということは 

|度|経度での距離|緯度での距離|
|---|---|---|
|1度|91km|111km|
|0.1度|9.1km|11.1km|
|0.01度|910m|1.11km|
|0.001度|91m|111m|
|0.0001度|9.1m|11.1m|
|0.00001度|91cm|1.11m|
|0.000001度|9.1cm|11.1cm|

 「程度である」ととらえてよいと思います。

つまり float32 で経度を表したときは 1〜2m 程度の精度しかないということです。

一方、緯度は 整数部分が 32、2進数で `0010 0000` なので、整数部は先頭の1を除いて5bit。小数部分には 18bit があてられるので、上の表から 0.0000038.. あたりの精度、30〜40cm くらいの精度はとれそうです。

これくらいならば 「道路をほかの道路と間違えないくらいの精度」としては、ありといえばありではないかというのが結論でした。

# まとめ

ざっくり、float32 ならば経度は 1-2m、緯度は30-40cm程度の誤差。float64 でメモリ使用量が多すぎるときは float32 にすることでメモリ使用量を半分にできる

なお日本に限れば、経度からオフセット 135 を引いてやれば整数部分が 緯度同様の2桁に収まるので、有効桁数が10進で1桁稼げ、距離に換算して30-40cm前後のもう少し高い精度がもてそうです。

# 補足

[S2](https://s2geometry.io/devguide/s2cell_hierarchy) という緯度経度の表現方法があります。

これは地球の球面を立方体に投影し、各面を4分割することを最大30分割くりかえす、その分割されたセルに番号をつけるというものです。

しかも、[ヒルベルト曲線](https://ja.wikipedia.org/wiki/%E3%83%92%E3%83%AB%E3%83%99%E3%83%AB%E3%83%88%E6%9B%B2%E7%B7%9A)というフラクタル的なカクカクの曲線上の位置(1次元)を使って採番しているため、「近い番号は近い場所を表す」という性質があり、検索や範囲指定に使いやすいのです。

その番号(S2 Cell ID)は64bit。

- 6面のうちどの面を表すのかに3bit(0～5)
- 各面ごとに4分割(2bit)を最大30回で60bit
- 分割をどこまでやるかを示す1の表現に1bit
を使いまです。


これを使えば 64bit で緯度経度それぞれ 9mm くらいの精度が出せるとのこと。こっちのほうがいいんじゃないか。。。