GeoHexのHex間距離

GithubでGeoHex3のpythonライブラリのissuesに、英語でGeoHexの1辺と距離に関する質問がつきました。
これに英語でスラッと答えられるほど英語力ないのですが、上記の関係については以前、現社内でまとめたことがあるので、そっちから引用してみようと思います。
当該の質問者には、自動翻訳か何かで見てもらう感じで…。

赤道上でLevel0の時の、GeoHexのHex間距離

Level0の時、GeoHexは地球の緯線上に（互い違いに）18個並びます。
18個並ぶという事は、互い違いに並ぶので、 Hexの1辺計算では27辺並ぶ 事になります。

赤道上での地球の1周距離 は、正確ではありませんがGeoHexがGoogleメルカトルと呼ばれる地球の球体近似（赤道半径、極半径とも6378137mの真球）を行っているので、それに従って計算すると、 40075016.68557849m

赤道上に27辺並ぶわけなので、 1辺の長さは 1484259.8772436476m

辺の中心からHexの中心までの長さ は、正三角形の高さ計算なので、1484259.8772436476 * sqrt(3) / 2 = 1285406.7595109711m

Hex中心間距離 は、この2倍なので 2570813.5190219423m

異なるレベルで赤道上の時の、GeoHexのHex間距離

GeoHexはレベルが1つ上がるたびに、 Hexの面積は9分割 されます。

面積が9分割されるという事は、 距離的なスケールは3分の1 になります。

よって、 レベルNの時の赤道上でのHex間距離 は、 2570813.5190219423 / 3^N m になります。

例： レベル7での赤道上でのHex間距離 は、 2570813.5190219423 / pow(3,7) = 1175.497722460879 m

異なる緯度での、GeoHexのHex間距離

（真球近似の元で） 地球上の緯線の長さ は、赤道上を1とすると、 緯度θのところでは、cosθ

それがメルカトル図法のためどこの緯度でも同じ緯線長さ（=1）になっているので 実際の長さに直すには、cosθを掛ければよい

よって、 レベルNの時の、緯度θ付近でのHex間距離 は、 2570813.5190219423 / 3^N * cosθ m になります。

日本周辺での値

日本周辺での、レベル11の時のHex間距離を求めると、

北方 緯度45度付近： 2570813.5190219423 / pow(3,11) * cos(45*π/180)　＝ 10.261758158289256 m

南方 緯度25度付近： 2570813.5190219423 / pow(3,11) * cos(25*π/180)　＝ 13.152626413705448 m

平均: 11.707192285997351 m （約11.7m）

一般的な雑居ビルのスケール感と比較してみると、

まあそんなもんなんすかね。

おまけ：GeoHexの「へクス間個数」

GeoHexはその性質上、 へクス間のへクスに沿って辿っていった個数 （大戦略のユニットのように）は、比較的簡単に求められます。
（性質は、 http://geogames.net/labs/geohex こちらのページで説明されています。このページのはGeoHex V.1のため、コード体系は異なりますが、ヘクス座標系に関する考え方はV3と一緒です。）
なんらかのGeoHex実装（後述のストアドファンクションgetXYByCode等）を用いれば、GeoHexコードをヘックス座標系のX,Y値に直せますが、

2つのGeoHex Hex1, Hex2があるとして、そのへクス座標が(X1, Y1), (X2, Y2)とした場合、その間のへクス個数を出す計算は、

• まず、X, Yそれぞれの差、ΔX = X1 - X2, ΔY = Y1 - Y2 を出します。
• ΔX と ΔYの正負が同じ（或いは片方が0）の場合
  • 絶対値の大きいほうが間のへクス数になります。
    （たとえば、ΔX > 0, ΔY > 0で、abs(ΔX) > abs(ΔY) ならば、abs(ΔX)が間のヘクス個数）
• ΔX と ΔYの正負が異なる場合
  • 絶対値の和が間のヘクス数になります。
    （たとえば、ΔX > 0, ΔY < 0 ならば、abs(ΔX) + abs(ΔY)が間のヘクス個数）

確認方法

http://geohex.net/v3.html のページで確認できます。
ヘクス指定も、
http://geohex.net/v3.html?code=XM02432111132
のような形ででき、ヘクスX座標, Y座標も表示できます。

これによると、XM02432111132は、 http://geohex.net/v3.html?code=XM02432111132 でヘクス座標は805913/-395120、
XP21344423331は http://geohex.net/v3.html?code=XP21344423331 でヘクス座標は984068/-301764

ΔX = 805913 - 984068 = -178155
ΔY = -395120 - (-301764) = -93356

正負が一緒なので、絶対値が大きいほうを採用して、ヘクス間ヘクス数178155個
レベル11でのヘクス間距離が11.707192285997351 mなので、ヘクスに沿った距離：11.707192285997351 * 178155 = 2085694.841711858 = 約2085kmくらいでしょうか。

mySQL用ストアドファンクション getXYByCode, getDistance

mySQLのストアドファンクションにgetXYByCodeとHEX単位の距離を求める関数を実装しました。

sql

delimiter //

drop procedure if exists proc_getXYByCode;
create procedure proc_getXYByCode(in code varchar(15), out h_x int, out h_y int)
begin
  declare H_KEY text default 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz';
  declare H_BASE real default 20037508.34;
  declare lv int;
  declare h_dec9, h_dec3 varchar(100) default '';
  declare i, l, xx, yy, dif int;
  declare h_dec0 varchar(2);
  declare h_pow real;
  declare max_hsteps, hsteps int;

  set h_x = 0;
  set h_y = 0;
  set lv = length(code) - 2;
  set h_dec9 = concat((instr(H_KEY, substr(code, 1, 1)) - 1)* 30 + (instr(H_KEY, substr(code, 2, 1)) - 1), substr(code, 3));
  if (h_dec9 regexp '^5[^125][^125]') then
    set h_dec9 = concat('7', substr(h_dec9, 2));
  elseif (h_dec9 regexp '^1[^125][^125]') then
    set h_dec9 = concat('3', substr(h_dec9, 2));
  end if;

  set h_dec9 = right(concat('0000', h_dec9), lv+3);
  set i = 0, l = length(h_dec9);
  while ( i < l ) do
    set h_dec3 = concat(h_dec3, right(concat('00', conv(substr(h_dec9, i+1, 1), 10, 3)), 2));
    set i = i + 1;
  end while;
  set i = 0, l = lv + 2;
  while ( i <= l ) do
    set h_pow = pow(3, lv + 2 - i);
    set xx = substr(h_dec3, i * 2 + 1, 1);
    set yy = substr(h_dec3, i * 2 + 2, 1);
    case xx 
      when '0' then set h_x = h_x - h_pow;
      when '2' then set h_x = h_x + h_pow;
      else begin end;
    end case;
    case yy
      when '0' then set h_y = h_y - h_pow;
      when '2' then set h_y = h_y + h_pow;
      else begin end;
    end case;
    set i = i + 1;
  end while;

  -- ここから _adjustXY()
  set max_hsteps = pow( 3, lv + 2);
  set hsteps = abs( h_x - h_y );

  if (hsteps = max_hsteps and h_x > h_y) then
    set xx = h_x;
    set h_x = h_y;
    set h_y = xx;
  elseif (hsteps > max_hsteps) then
    set dif = hsteps - max_hsteps;
    set xx = floor(dif / 2);
    set yy = dif - xx;
set @xx = xx;
set @yy = yy;
    if (h_x > h_y) then
      set dif = h_x;
      set h_x = h_y + yy + xx;
      set h_y = dif - xx - yy;
    elseif ( h_y > h_x ) then
      set dif = h_x;
      set h_x = h_y - yy - xx;
      set h_y = dif + xx + yy;
    end if;
  end if;
  -- 戻値 h_x, h_y;
end//

drop function if exists func_getDistance//
create function func_getDistance(code1 varchar(15), code2 varchar(15)) returns integer
begin
  declare x1, y1, x2, y2 integer;
  call proc_getXYByCode(code1, x1, y1);
  call proc_getXYByCode(code2, x2, y2);

  set x1 = x1 - x2;
  set y1 = y1 - y2;

  if (sign(x1) = sign(y1)) then
    case when abs(x1) > abs(y1) then return x1;
         else return y1;
    end case;
  else
    return abs(x1) + abs(y1);
  end if;
end//

delimiter ;