OSM(OpenStrreMap)の海岸線データを繋ぎ合わせて、陸地ポリゴンを形成できた。 しかし、特に大陸は巨大なポリゴンとなるため、 これを zoom 9, 10程度のタイルの大きさを目安として、多数の小さなポリゴンに分割したい。
先日の最初のテストでは、簡単な方法では分割できないという結果であったが、 これは交差判定の引数の数値が非常に大きいことが影響したようである。
自前のプログラムで、改めて水平分割を試すと、そこそこに分割候補線があることが分かった。
水平分割の場合、東西を分割線で結び途中海岸線と交差がない場所を選ぶ。分割により、上部に一つの ポリゴンがスライスされる。
水平分割を終えた後、大きいポリゴンは垂直分割を行う。これを２度繰り返すと、 大きいポリゴンは少なくなるはずである。
矩形のポリゴンはそれ以上細分しなくてもいいであろう。
ここの分割操作では、ポリゴンを２分割するだけであり、３つ以上のポリゴンが生まれることはないため、 プログラムはシンプルである。
分割線と海岸線の２箇所の交点で、ポリゴンを二つに分割する。分割により、新しいポリゴンができるだけでなく、 元のポリゴンは小さくなるため、ノードの削除も必要となる。 また、交点が新しいノードとなるが、ノードの挿入ではなく、元のノードの値の変更として扱う。
以上のような特徴のあるポリゴンのため、独自 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 となる。
しかし、下のプログラムでは１件だけ逆があった。たった１件のため、この原因の究明は先に延ばす。
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;
}
ポリゴン分割の見込みはついた。スマートにポリゴンレコードを出力するために、 ここで、一からプログラムを書き直す。
分割により、時計回りになったのではないか？カスピ海が陸地として描画されている。
元のポリゴンの先頭がどこにあったかが問題である。始点と終点が同じというのがガンになる。 毎回、元のファイルを作り変えることは避けたい。
始点が分割される側か残る側かにより、処理を分けねばならないだろう。 最初に始点が一番南に位置するようにデータをずらす方が分かりやすい気がする。
rotateメソッドを導入、iXmaxや iYmax が 0 になるようにする。当然、他方の修正も必要である。
修正により、ファイルサイズは妥当な値になった。
プログラムを大幅に書き換えたため、やはり、バグが入り込んだ。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 が抜けたかも知れない。