😀

地域メッシュコードからポリゴンを生成するPostGIS関数

2024/07/02に公開

postgis

GIS

tech

はじめに

地域メッシュは JIS X 0410:2002 に規定されているようです。緯度経度ベースで隙間なくメッシュ（網目）に分けたものです。メッシュごとに個別のデータが入れます。入るデータは人口などの統計情報であったり、標高データであったりします。

地域メッシュは1次から3次までの基準地域メッシュと、より細かい分割地域メッシュとに分かれます。

メッシュサイズは、緯度経度ベースなので、正確な距離は出ません (緯度に依存します)。1次メッシュは1辺約80kmで、2次は約10km、3次は約1kmです。分割地域メッシュには、3次メッシュを緯度・経度で各2分割した1/2メッシュ、さらに各2分割した1/4メッシュ、さらに分割した1/8メッシュがあります。

詳細は https://www.stat.go.jp/data/mesh/pdf/gaiyo1.pdf などを参考にして下さい。

地域メッシュ単位で提供される統計データには、メッシュコードは必ず入っていますが、それに対応した地物データ (ポリゴン) と関連付けられているとは限りません。メッシュコードは分かってるけど、地図に重ね合わすことができるレベルで、メッシュの位置やサイズが分かるでしょうか。私は分かりません。

また、たとえば https://www.e-stat.go.jp/ では、メッシュコードと属性データを納めたCSVファイルと、メッシュコードと地物データとを組み合わせたGISデータとの個別配布はしてくれていますが、それぞれを手動でダウンロードする必要があります。

もっと手軽にメッシュコードのデータを GIS で扱うには、メッシュコードからポリゴンを生成する関数が必須だと思いますよね? ね?

ということで PostGIS で使える、地域メッシュコードからポリゴンを生成する関数を作成したので、その関数について解説します。

ソースのありか

https://github.com/boiledorange73/pg_jpgrid にあります。

メッシュコードの正規表現

3次メッシュまで

まずは文字列がコードかどうかを確認する必要があるので、計算します。
1次メッシュは2桁、2桁なので ([0-9]{2})([0-9]{2}) となります。
2次メッシュは1次メッシュの後に、1桁、1桁を繋げます。ただし、マイナスをデリミタにする場合もあるのを考慮する必要があります。まとめると ([0-9]{2})([0-9]{2})-?([0-9])([0-9]) となります。
3次メッシュは ([0-9]{2})([0-9]{2})-?([0-9])([0-9])-?([0-9])([0-9]) となりますね。

最後に、先頭から末尾まで余計なものが入らず、かつ、空白は受け付けるようにするための規則を適宜入れると、^\s*([0-9]{2})([0-9]{2})\s*-?\s*([0-9])([0-9])\s*-?\s*([0-9])([0-9])\s*$ となります。

3次より細かくする場合

3次メッシュより細かい場合には、1-4の数字(1桁)だけで南北2分割、東西2分割で、あわせて4分割を行っていきます。マイナスはあってもなくてもいいのと、空白を考慮すると (\s*-?\s*[1-4]\s*) となります。

どこまでも細分できるようにするために ((\s*-?\s*[1-4]\s*)+) としました。この場合は、3次より細かい部分のコードが一つにまとまるので、別途分割する必要があります。

ただし (-?[1-4]) がずっと続き、かつマイナスが意味のない文字ですので、マイナスを取り除いて [1-4]の文字を1文字ずつ見ていく方針にしました。分割自体はせずに、マイナスを消す置換を行い、1文字ずつ見ていくようにしています。

まとめると ^\s*([0-9]{2})([0-9]{2})\s*-?\s*([0-9])([0-9])\s*-?\s*([0-9])([0-9])((\s*-?\s*[1-4]\s*)+)$ となります。

南西隅緯度経度と幅を求める

arr は、3次までは、南北のコード、東西のコードの、二つの数字の組合せなので、配列長は次数の2倍になります。

1次メッシュ

1次メッシュの南西隅緯度経度と幅は次のように求めることにないります。

    reslen := array_length(match_result, 1);

    IF reslen >= 2 THEN
        latsecmin := arr[1] * 2400;
        lngsecmin := (arr[2] + 100) * 3600;
        dlatsec := 2400;
        dlngsec := 3600;
    END IF;

2次メッシュと3次メッシュ

2次メッシュは、南西隅緯度経度は1次メッシュの結果に加算して導き、幅は上書きするだけです。

    IF reslen >= 4 THEN
        latsecmin := latsecmin + arr[3] * 300;
        lngsecmin := lngsecmin + arr[4] * 450;
        dlatsec := 300;
        dlngsec := 450;
    END IF;

3次メッシュも同様です。

    IF reslen >= 6 THEN
        latsecmin := latsecmin + arr[5] * 30;
        lngsecmin := lngsecmin + arr[6] * 45;
        dlatsec := 30;
        dlngsec := 45;
    END IF;

3次より細かい場合

3次より細かい場合には、先に match_result から除いて frac に、数字だけでなる文字列入れています。

    IF reslen >= 7 THEN
        frac := regexp_replace(match_result[7], '-', '', 'g');
        match_result := match_result[1:6];
        reslen := 6;
    END IF;

数字だけでなっているので、1文字ずつ substr で抽出して処理しています。
幅は半分に割っていき、南西隅緯度経度を、1なら変えず、2なら経度だけ変え、3`緯度だけ変え、4``なら緯度経度ともに変える、というのを繰り返しています。

    -- 4桁目以降の処理
    frac_len := length(frac);
    -- RAISE NOTICE 'frac % frac_len %', frac, frac_len;
    FOR frac_n IN 1..frac_len LOOP
        div := (1 << frac_n)::DOUBLE PRECISION;
        -- RAISE NOTICE 'div %', div;
        dlatsec := 30.0 / div;
        dlngsec := 45.0 / div;
        frac_one := substr(frac, frac_n, 1)::INTEGER;
        -- 2,4 -> east (1,3 -> west)
        IF mod(frac_one, 2) = 0 THEN
            lngsecmin := lngsecmin + dlngsec;
        END IF;
        -- 3,4 -> north (1,2 -> south)
        IF frac_one >= 3 THEN
            latsecmin := latsecmin + dlatsec;
        END IF;
    END LOOP;

秒単位から度単位に変換する

この関数では、緯度、経度の単位は秒で計算して、最後に 3600 で割って度単位にしています。

ポリゴンを生成する

関数の宣言は次のようになっています。

CREATE OR REPLACE FUNCTION jpgridcode2polygon(code TEXT, srid INTEGER DEFAULT 4326) RETURNS GEOMETRY(POLYGON) AS $$
...

srid のデフォルト値を 4326 (WGS84) としています。ST_SetSRID()で付けているので、変換はしません。
また、関数の返り値の型に GEOMETRY(POLYGON) と、ポリゴンであることを示す型修飾をつけています。あと、SRIDの型修飾も可能です。

この関数のやってることは、南西隅緯度経度とメッシュ幅からポリゴンを作っているだけです。

    RETURN ST_SetSRID(
        ST_MakePolygon(
            ST_MakeLine(ARRAY[
                ST_MakePoint(arr[2],arr[1]),
                ST_MakePoint(arr[4],arr[1]),
                ST_MakePoint(arr[4],arr[3]),
                ST_MakePoint(arr[2],arr[3]),
                ST_MakePoint(arr[2],arr[1])
            ])
        ),
        srid
    );

こういうのができました

1/8メッシュごとの人口総数をGIS上に表示して見ました。

人口総数の図

おわりに

いかがだったでしょうか。

ここでは、地域メッシュコードの正規表現と、南西隅緯度経度と南北・東西方向のメッシュ幅を求める方法を示しました。また、そこからポリゴンを作る関数の作成までを紹介しました。

これで、メッシュコードを含む非GISデータとともにこの関数を呼ぶビューを作ると、非GISデータを GIS でサッと見ることができるようになります。

もちろんこれだけだとあまり意味がないので、統計GISから落としたデータに組み合わせて使ってみたいところですね。そのうち何かします。

出展

政府統計の総合窓口(e-Stat) を使用しました。
地理院地図を使用しました。

本記事のライセンス

クリエイティブ・コモンズ・ライセンス
この記事はクリエイティブ・コモンズ表示 4.0 国際ライセンスの下に提供されています。