【BigQuery 】緯度経度から50mメッシュコードを算出する方法

ぐぐっても出なかったので自分でSQL書きました。

ただ、未検証なのでご了承のほどよろしくお願いいたします。

BigQuery

WITH
  t0 AS (
  SELECT
    139.6917337 AS longitude,
    35.6895014 AS latitude ),
  t1 AS (
  SELECT
    latitude,
    longitude,
    /* FLOOR X 以下で最大の整数値を返します。*/ 
    /*round X のみが存在する場合、X を最も近い整数に丸めます。N が存在する場合、X を小数点以下の N 桁に丸めます。*/ 
    /*MOD(X, Y) Y による X の除算の剰余を返します。*/ 
    /* 1次メッシュ 辺の長さ：約80km 経度差：1度 緯度差：40分*/ 
    CONCAT( SAFE_CAST(SAFE_CAST(FLOOR(latitude*60/40) AS INT64) AS STRING), SAFE_CAST(SAFE_CAST(FLOOR(longitude)-100 AS INT64) AS STRING) ) AS mesh1d,
    /* 2次メッシュ 辺の長さ：約10km 経度差：7分30秒 緯度差：5分*/ 
    CONCAT( SAFE_CAST(SAFE_CAST(FLOOR(MOD(SAFE_CAST(ROUND(latitude*60, 7) AS NUMERIC), 40.0)/5) AS INT64) AS STRING), SAFE_CAST(SAFE_CAST(FLOOR(MOD(SAFE_CAST(longitude-100 AS NUMERIC), 1.0)*60/7.5) AS INT64) AS STRING) ) AS mesh2d,
    /*3次メッシュ 辺の長さ：約1km 経度差：45秒 緯度差：30秒 */ 
    CONCAT( SAFE_CAST(SAFE_CAST(FLOOR(MOD(MOD(SAFE_CAST(ROUND(latitude*60, 7) AS NUMERIC),40.0),5.0)* 60/30) AS INT64) AS STRING), SAFE_CAST(SAFE_CAST(FLOOR(MOD(MOD(SAFE_CAST(longitude AS NUMERIC),1.0)*60,7.5) * 60/45) AS INT64) AS STRING) ) AS mesh3d,
    /* 2分の1地域メッシュ 辺の長さ：約500m 経度差：22.5秒 緯度差：15秒 三次メッシュを緯線方向、経線方向に2等分してできる区域*/ 
    CONCAT( SAFE_CAST(SAFE_CAST(FLOOR(MOD(MOD(MOD(SAFE_CAST(ROUND(latitude*60, 7) AS NUMERIC),40.0),5.0),0.5)*4)*2+ FLOOR(MOD(MOD(MOD(SAFE_CAST((longitude-100) AS NUMERIC),1.0)*60,7.5),0.75)* 60/22.5)+1 AS INT64) AS STRING) ) AS mesh4d,
    /* 4分の1地域メッシュ 辺の長さ：約250m 経度差：11.25秒 緯度差：7.5秒 */ 
    CONCAT( SAFE_CAST(SAFE_CAST(FLOOR(MOD(MOD(MOD(MOD(SAFE_CAST(ROUND(latitude*60, 7) AS NUMERIC),40.0),5),0.5),0.25)*8)*2+ FLOOR(MOD(MOD(MOD(MOD(SAFE_CAST((longitude-100) AS NUMERIC),1)*60,7.5),0.75),0.375)*60/11.25)+1 AS INT64) AS STRING) ) AS mesh5d,
    /* 8分の1地域メッシュ 辺の長さ：約125m 経度差：5.625秒 緯度差：3.75秒 */
     CONCAT( SAFE_CAST(SAFE_CAST(FLOOR(MOD(MOD(MOD(MOD(MOD(SAFE_CAST(ROUND(latitude*60, 7) AS NUMERIC),40.0),5.0),0.5),0.25),0.125)*16)*2+ FLOOR(MOD(MOD(MOD(MOD(MOD(SAFE_CAST((longitude-100) AS NUMERIC),1)*60,7.5),0.75),0.375),0.1875)* 60/5.625)+1 AS INT64) AS STRING) ) AS mesh6d,
    /* 10分の1地域メッシュ 辺の長さ：約100m 経度差：4.5秒 緯度差：3秒 三次メッシュを緯線方向、経線方向に10等分してできる区域*/ 
    CONCAT( SAFE_CAST(SAFE_CAST(FLOOR(MOD(MOD(MOD(SAFE_CAST(ROUND(latitude*60, 7) AS NUMERIC),40.0),5.0),0.1)*20)*2+ FLOOR(MOD(MOD(MOD(SAFE_CAST((longitude-100) AS NUMERIC),1.0)*60,7.5),0.75)* 60/4.5)+1 AS INT64) AS STRING) ) AS mesh7d,
    /* 20分の1地域メッシュ 辺の長さ：約50m 経度差：2.25秒 緯度差：1.5秒 1/10メッシュを緯線方向、経線方向に2等分してできる区域*/ 
    CONCAT( SAFE_CAST(SAFE_CAST(FLOOR(MOD(MOD(MOD(MOD(SAFE_CAST(ROUND(latitude*60, 7) AS NUMERIC),40.0),5.0),0.1),0.05)*40)*2+ FLOOR(MOD(MOD(MOD(MOD(SAFE_CAST((longitude-100) AS NUMERIC),1.0)*60,7.5),0.75),0.075)* 60/2.25)+1 AS INT64) AS STRING) ) AS mesh8d
  FROM
    t0 )
SELECT
  latitude,
  longitude,
  CONCAT(mesh1d,mesh2d,mesh3d,mesh4d,mesh5d,mesh6d,mesh7d,mesh8d) AS mesh50m
FROM
  t1

参考

• 緯度経度から1/8（6次）メッシュを求めるSQL