LoginSignup
10
11

More than 5 years have passed since last update.

@kochizufan(Kohei OTSUKA)

GeoHexのHex間距離

Posted at

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

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

geohex_level_0.png

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_level_11.png

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

おまけ：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 ;
10
11
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
10
11