トップ地図システム基礎技術 > 大きなポリゴンを多数の小さなポリゴンに分割する

大きなポリゴンを多数の小さなポリゴンに分割する

前ページ

はじめに

OSM(OpenStrreMap)の海岸線データを繋ぎ合わせて、陸地ポリゴンを形成できた。 しかし、特に大陸は巨大なポリゴンとなるため、 これを zoom 9, 10程度のタイルの大きさを目安として、多数の小さなポリゴンに分割したい。

先日の最初のテストでは、簡単な方法では分割できないという結果であったが、 これは交差判定の引数の数値が非常に大きいことが影響したようである。

自前のプログラムで、改めて水平分割を試すと、そこそこに分割候補線があることが分かった。

水平分割の場合、東西を分割線で結び途中海岸線と交差がない場所を選ぶ。分割により、上部に一つの ポリゴンがスライスされる。

水平分割を終えた後、大きいポリゴンは垂直分割を行う。これを2度繰り返すと、 大きいポリゴンは少なくなるはずである。

矩形のポリゴンはそれ以上細分しなくてもいいであろう。

ここの分割操作では、ポリゴンを2分割するだけであり、3つ以上のポリゴンが生まれることはないため、 プログラムはシンプルである。

ポリゴンの管理

分割線と海岸線の2箇所の交点で、ポリゴンを二つに分割する。分割により、新しいポリゴンができるだけでなく、 元のポリゴンは小さくなるため、ノードの削除も必要となる。 また、交点が新しいノードとなるが、ノードの挿入ではなく、元のノードの値の変更として扱う。

以上のような特徴のあるポリゴンのため、独自 class Polygon を導入する。

最初のバージョンを以下に示す。

ノードは例えば class Node とする案もあるが、プログラムを複雑にしないため、X座標列とY座標列として管理する。 Listにすることも考えたが、総合的な分かりやすさから、配列とした。

OSMのポリゴンでは、先頭ノードと末尾ノードが同じである。OSMと違えるのはトラブルの原因となる恐れがあるため、 class Polygon もこれに合わせた。

配列 xs、ys は循環配列として操作する。先頭と末尾が同じため、通常の循環配列と少し異なる。 配列の長さは length であるが、循環配列として操作するときは、最後の要素は無視して、 配列の長さは length - 1 とする。

ポリゴンは時計回りであることは保証されているが、ノード (xs[0]、ys[0]) がどこに位置するかは決まっていない。 このため、X座標あるいはY座標が最大となるノード番号 iXmax、iYmax を求めている。 分割処理では、水平分割の時は iXmax、垂直分割の時は iYmax から循環配列を一巡するのが分かりやすい。 

    class Polygon {
        int length;         // 実際のノード数は length - 1
        int size;           // 配列のサイズ
        int[] xs;           // xs[0] == xs[length-1]
        int[] ys;           // ys[0] == ys[length-1]
        int xmin, ymin, xmax, ymax;
        int iXmax, iYmax;   // xs[iXmax] == xmax, ys[iYmax] == ymax

        public Polygon(int size) {
            this.size = size;
            xs = new int[size];
            ys = new int[size];
            length = 0;
        }

        public Polygon(int[] xs, int[] ys) {
            length = size = xs.length;
            this.xs = new int[size];
            this.ys = new int[size];
            System.arraycopy(xs, 0, this.xs, 0, length);
            System.arraycopy(ys, 0, this.ys, 0, length);
            setMinMax();
        }

        int next(int i) { return i < length-2 ? i + 1 : i == length-2 ? 0 : 1; }
        
        void setMinMax() {
            xmin = xmax = xs[0];
            ymin = ymax = ys[0];
            iXmin = iYmin = 0;
            for (int i = 1; i < length - 1; i++) {
                if (xs[i] < xmin) { xmin = xs[i]; }
                if (xs[i] > xmax) { xmax = xs[i]; iXmax = i; }
                if (ys[i] < ymin) { ymin = ys[i];  }
                if (ys[i] > ymax) { ymax = ys[i]; iYmax = i;}
            }
        }

        // 値を変更する
        void set(int ix, int x, int y) {
            if (ix >= 0 && ix < length) {
                xs[ix] = x;
                ys[ix] = y;
            }
        }

        void add(int x, int y) {    // 末尾に追加
            if (length < size) {
                xs[length] = x;
                ys[length] = y;
                length++;
            }
        }

        // from は含む、to は含まない
        void add(Polygon poly, int from, int to) {
            if (length + (to - from) <= size) {
                for (int i = from; i < to; i++) {
                    xs[length] = poly.xs[i];
                    ys[length] = poly.ys[i];
                    length++;
                }
            }
        }

        // from は含む、to は含まない
        void remove(int from, int to) {
            System.arraycopy(xs, to, xs, from, length-to);
            System.arraycopy(ys, to, ys, from, length-to);
            length -= to - from;
        }
    }

ポリゴンの分割処理

水平分割のときは、Y座標が小さい上からスライスするため、iYmax から始めれば、 スタートノードは常に下にある。

OSMのポリゴンは反時計回りが基準のため、分割線(水平線)との交点のX座標を順に x1、x2 とすると、 x1 > x2 となる。

しかし、下のプログラムでは1件だけ逆があった。たった1件のため、この原因の究明は先に延ばす。

    static List<Polygon> splitPolygonHorizontally(List<Polygon> src, int span) {
        List<Polygon> dst = new ArrayList<>();
        for (Polygon poly : src) {
            int bgn = ((poly.ymin + span) / span) * span;
            for (int y = bgn; y < poly.ymax; y += span) {  // y : 水平分割線のY座標
                int cntCross = 0;
                int i1 = -1, i2 = -1;
                int k = poly.iYmax;
                do {
                    if ((poly.ys[k] - (long)y) * (poly.ys[k+1] - y) < 0) {  // 交差判定
                        cntCross++;
                        if (i1 < 0) i1 = k;  // 最初の交点 i1とi1+1の間
                        else i2 = k;            // 最後の交点 i2とi2+1の間
                    } 
                    k = poly.next(k);           // 末尾の次は先頭
                } while (k != poly.iYmax) ;     // 一巡するまで

                if (cntCross == 2 && i2 - i1 >= 2) { // 水平分割可能
                    int x1 = crossHorizonX(poly.xs[i1], poly.ys[i1], poly.xs[i1+1], poly.ys[i1+1], y);
                    int x2 = crossHorizonX(poly.xs[i2], poly.ys[i2], poly.xs[i2+1], poly.ys[i2+1], y);
                    // (x1,y): 線分i1ーi1+1と水平線yの交点, (x2,y): 線分i2ーi2+1と水平線yの交点
                    if (x2 > x1) {
                        System.out.printf("Error x1=%d x2=%d\n", x1, x2);
                        continue;
                    }

                    //分割で生まれるポリゴン x1, i1+1, ... , i2, x2, x1
                    Polygon pN = new Polygon(i2 - i1 + 3);
                    pN.add(x1, y);
                    for (int i = i1+1; i <= i2; i++) {
                        pN.add(poly.xs[i], poly.ys[i]);
                    }
                    pN.add(x2, y);
                    pN.add(x1, y);
                    pN.setMinMax();
                    dst.add(pN);

                    // 元のポリゴン(小さくなる)
                    poly.set(i1+1, x1, y);
                    poly.set(i2, x2, y);
                    poly.remove(i1+2, i2);  // i1+2, i1+3, ... , i2-1 を削除
                    poly.setMinMax();
                }
            }
            dst.add(poly);	// 最後の分割で残った南側のポリゴンをリストに追加する
        }
        return dst;
    }

おわりに

ポリゴン分割の見込みはついた。スマートにポリゴンレコードを出力するために、 ここで、一からプログラムを書き直す。

リファレンス

[1] Shapely を使用して、グリッドを形成する同じサイズの長方形にポリゴンを分割できます。
[2] osmcode / osmcoastline
[3] ST_Subdivide
[4] ポリゴン パーティション
[5] Splitting an arbitrary polygon by a line
[6] java.awt.geom.Line2D

来歴およびノート

2023.2.3 日本地図は描画されるようになったが、大陸が描画されない。

分割により、時計回りになったのではないか?カスピ海が陸地として描画されている。

元のポリゴンの先頭がどこにあったかが問題である。始点と終点が同じというのがガンになる。 毎回、元のファイルを作り変えることは避けたい。

始点が分割される側か残る側かにより、処理を分けねばならないだろう。 最初に始点が一番南に位置するようにデータをずらす方が分かりやすい気がする。

rotateメソッドを導入、iXmaxや iYmax が 0 になるようにする。当然、他方の修正も必要である。

2023.2.3 分割によりファイルサイズが少し小さくなった。
少し大きくなるべきのため、どこかにバグがある。

修正により、ファイルサイズは妥当な値になった。

2023.2.3 より分かりやすく大幅修正
改めてまず分割なしでテスト。メモリ使用量の制約から、面積の小さいものだけをレンダリングした。 日本の場合、本州を除いて正しくレンダリングできた。よって、Polygon#writePolygon メソッドは問題なしであろう。
2023.2.2 プログラム大幅書き換え

プログラムを大幅に書き換えたため、やはり、バグが入り込んだ。Androidスマホでレコードをチェックする。

まず、空間検索の入り口で、ヘッダとレコード長をチェックした。 ヘッダーは 16進数で 0x22 のため10進では 34 でよい。この二つは OK である。

I/System.out: #2 /storage/5BF6-842C/Map/lands1/1/1.dat 3ms 5KB
D/Debug: ix=2 head=34 length=1975
 ix=1979 head=34 length=3925
 ix=5906 head=34 length=79

次に、空間検索で抽出されたレコードの境界ボックスを調べた。 zoom 2 ではつぎのようになった。 XY座標値は 0~230 の範囲にあるので、日本ではこのような値になるのであろう。

 (167758265 368809538) (167763213 368813514)
 (167903256 369087940) (167904402 369089086)

zoom 9では伊豆大島など島だけは正常に描画された。よって、バグは多くあるにせよ、 何も描画されないわけではない。

分割のいらない小さな島については問題ないようだ。つまり、Polygonの出力は問題なく、 分割処理に問題があるようである。

ファイルサイズ上は妥当な値である。データに誤りがあるのであろう。

念のため、分割を止めてみる。

分割するかデータを間引かないとメモリ不足がおきる。一瞬表示されたので問題ないだろう。 一部のポリゴンが逆回りかもしれない。reverse したはずだが、修正でこの reverse が抜けたかも知れない。