OSMデータ(xml形式)は node, way, relation という階層構造である。座標値は node にある。 way は line や polygon を表す。ここには座標値はなく node id 列がある。 relation のメンバーは node, way, relation である。
構造的には極めてシンプルである。全データをメモリに取り込める場合、 node、way、relation の構造体を定義して、そこにデータを取り込めばよい。 relation は メンバーに relation を含むため、再帰構造となる。
しかし、圧縮ファイルでも 50GB の大きさである。 上記の構造体に全データを取り込むには数百GB～数TBのメモリが必要となる。
つまり、構造はシンプルであるが、メモリに取り込めないため、その操作が非常に厄介である。
間引きでは Ramer-Douglas-Peuckerアルゴリズムが著名で、 PostGISのST_SimplifyPreserveTopologyでもこれが使われている。 PostGISには複数のpolygonが交差しているとき、交差のない部分だけを取り出す機能があり、 これを使って実質上交差を回避している。この実装は間引きよりも難しい。 特に multipolygon の交差回避が難しい。
Ramer-Douglas-Peuckerアルゴリズム（またはDouglas-Peuckerアルゴリズム）は 下に示すように、極めてシンプルなものである。
simplySection関数は線分 points[i]-points[j] の間引きを行う関数である。 許容誤差を epsilon とする。
線分 points[i]-points[j]と途中の点 points[k] との垂直距離を調べて、 その最大値を maxDistance とする。
この最大値が epsilon を超えなければ、全ての中間点は間引いて、処理は終わりである。
最大値が許容誤差を超える場合には、 この最大値を持った点 points[maxIndex] は間引かないことに決定する。 この点を境として、ふたつの線分 points[i]-points[maxIndex]および points[maxIndex]-points[j] に分けて、 それぞれについて simplySection関数を再帰コールする。
[C/C++言語]
// Returns the distance from point p to the line between p1 and p2
double perpendicular_distance(point_t p, point_t p1, point_t p2) {
double dx = p2.x - p1.x;
double dy = p2.y - p1.y;
double d = sqrt(dx * dx + dy * dy);
return fabs(p.x * dy - p.y * dx + p2.x * p1.y - p2.y * p1.x)/d;
}
void simplifySection(point_t *points, size_t i, size_t j, double epsilon) {
if ((i + 1) == j) return;
double maxDistance = -1.0;
size_t maxIndex = i;
for (size_t k = i + 1; k < j; k++) {
double distance = perpendicular_distance(points[k], points[i], points[j]);
if (distance > maxDistance) {
maxDistance = distance;
maxIndex = k;
}
}
if (maxDistance <= epsilon) {
for (size_t k = i + 1; k < j; k++) {
points[k].unused = true; // 間引く
}
} else {
simplifySection(points, i, maxIndex, epsilon);
simplifySection(points, maxIndex, j, epsilon);
}
}
[Java]
// １度＝約100km ⇒ 1m＝約0.00001度
public static double PerpendicularDistance(long p, long p1, long p2) {
double E7 = 10000000.0;
double px = getX(p) / E7;
double py = getY(p) / E7;
double p1x = getX(p1) / E7;
double p1y = getY(p1) / E7;
double p2x = getX(p2) / E7;
double p2y = getY(p2) / E7;
double dx = p2x - p1x;
double dy = p2y - p1y;
double d = Math.sqrt(dx*dx + dy*dy);
return Math.abs(px*dy - py*dx + p2x*p1y - p2y*p1x)/d;
}
public static void SimplifySection(long[] pnts, int i, int j, double epsilon) {
if ((i + 1) == j) return;
double maxDistance = -1.0;
int maxIndex = i;
for (int k = i + 1; k < j; k++) {
double distance = PerpendicularDistance(pnts[k], pnts[i], pnts[j]);
if (distance > maxDistance) {
maxDistance = distance;
maxIndex = k;
}
}
if (maxDistance <= epsilon) {
for (int k = i + 1; k < j; k++) {
pnts[k] = Long.MIN_VALUE; // 間引く
}
} else {
SimplifySection(pnts, i, maxIndex, epsilon);
SimplifySection(pnts, maxIndex, j, epsilon);
}
}
このように、再帰を使うことにより、効率よく間引きを行う。