メッシュコードについて

メッシュとは緯度経度に基づいて、地域を格子状に分割したもので、
メッシュコードとは、この地域メッシュに割り当てられたユニークな識別番号になります。
測地系は平成14年4月1日以前は日本測地系でしたが、それ以降は世界測地系に準じているようです。
この記事では、私がよく使うPostgreSQLとAWS Athena(Presto SQL)で緯度経度からメッシュコードを計算するクエリを紹介いたします。

計算方法の参考文献：http://www.stat.go.jp/data/mesh/pdf/gaiyo1.pdf#page=7

クエリ

紹介するクエリは以下のSQLで動作することを確認しております。

PostgreSQL
PrestoSQL(AWS Athena)

今回はlon, latがともに世界測地系であることを前提としています。

1次メッシュ

辺の長さ：約80km
経度差：1度
緯度差：40分

WITH t AS (
  SELECT
  139.71475 AS lon,
  35.70078 AS lat
)
SELECT
    CONCAT(
      CAST(CAST(floor(lat*60/40) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(lon-100) AS INTEGER) AS VARCHAR)
    ) AS mcode
FROM t

> 5339

2次メッシュ

辺の長さ：約10km
経度差：7分30秒
緯度差：5分

WITH t AS (
  SELECT
  139.71475 AS lon,
  35.70078 AS lat
)
SELECT
    CONCAT(
      CAST(CAST(floor(lat*60/40) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(lon-100) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor((lat*60)%40/5) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor((lon-floor(lon))*60/7.5) AS INTEGER) AS VARCHAR)
    ) AS mcode
FROM t

> 533945

3次メッシュ

辺の長さ：約1km
経度差：45秒
緯度差：30秒

WITH t AS (
  SELECT
  139.71475 AS lon,
  35.70078 AS lat
)
SELECT
    CONCAT(
      CAST(CAST(floor(lat*60/40) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(lon-100) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor((lat*60)%40/5) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor((lon-floor(lon))*60/7.5) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(((lat*60)%40)%5*60/30) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(((lon-floor(lon))*60)%7.5*60/45) AS INTEGER) AS VARCHAR)
    ) AS mcode
FROM t

> 53394547

2分の1地域メッシュ

辺の長さ：約500m
経度差：22.5秒
緯度差：15秒

WITH t AS (
  SELECT
  139.71475 AS lon,
  35.70078 AS lat
)
SELECT
    CONCAT(
      CAST(CAST(floor(lat*60/40) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(lon-100) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor((lat*60)%40/5) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor((lon-floor(lon))*60/7.5) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(((lat*60)%40)%5*60/30) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(((lon-floor(lon))*60)%7.5*60/45) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(((lat*60)%40)%5*60%30/15)*2 + floor(((lon-floor(lon))*60)%7.5*60%45/22.5)+1 AS INTEGER) AS VARCHAR)
    ) AS mcode
FROM t

> 533945471

4分の1地域メッシュ

辺の長さ：約250m
経度差：11.25秒
緯度差：7.5秒

WITH t AS (
  SELECT
  139.71475 AS lon,
  35.70078 AS lat
)
SELECT
    CONCAT(
      CAST(CAST(floor(lat*60/40) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(lon-100) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor((lat*60)%40/5) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor((lon-floor(lon))*60/7.5) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(((lat*60)%40)%5*60/30) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(((lon-floor(lon))*60)%7.5*60/45) AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(((lat*60)%40)%5*60%30/15)*2 + floor(((lon-floor(lon))*60)%7.5*60%45/22.5)+1 AS INTEGER) AS VARCHAR),
      CAST(CAST(floor(((lat*60)%40)%5*60%30%15/7.5)*2 + floor(((lon-floor(lon))*60)%7.5*60%45%22.5/11.25)+1 AS INTEGER) AS VARCHAR)
    ) AS mcode
FROM t

> 5339454711

プログラミングのように変数にできれば、もっと綺麗に計算できるのですがクエリだとどうしても力技になってしまいます。
しかし、テーブルにメッシュコードのカラムを付与することができれば簡単にメッシュごとの統計処理を行えるようになります。
是非ご検討ください。

SQLで緯度経度からメッシュコードを計算する