🗾
地理空間データのポリゴンを簡略化する
地理空間データのポリゴンを簡略化する方法を考えます。
有名なライブラリのturf.jsなどはRamerDouglasPeucker法というアルゴリズムで簡略化を行っています。
この方法は、選択した2点を結んだ線分からの距離が事前に決めた閾値下回る点を間引いていきます。
しかし、この閾値はどうやって決めるとよいでしょうか。
{
"type": "Polygon",
"coordinates": [
[
[0, 30],[1, 30.5],[4, 30],[6, 30.2],[8, 30.1],[9, 30.15],[12, 29.5],[40, 40],[23, 40],[27, 55],[25, 54],[20, 57],[0, 52],[-24, 53],[-24, 50],[-24, 47],[-24, 44],[-24, 41],[-24, 38],[0, 30]
]
]
}
{
"type": "Polygon",
"coordinates": [
[
[0, 30],[0.1, 30.05],[0.4, 30],[0.6, 30.02],[0.8, 30.01],[0.9, 30.015],[1.2, 29.95],[4, 31],[2.3, 31],[2.7, 32.5],[2.5, 32.4],[2, 32.7],[0, 32.2],[-2.4, 32.3],[-2.4, 32],[-2.4, 31.7],[-2.4, 31.4],[-2.4, 31.1],[-2.4, 30.8],[0, 30]
]
]
}
この2つの図形は相似の関係[1]にありますが、同じ閾値を使ってしまうと簡略化後は全く違う形状になってしまうでしょう。
図形の大小に関わらず扱える閾値で簡略化はできないでしょうか。
また簡略化の条件として、ポリゴンの頂点数をX点以下にしたいといったケースもあります。
これらを満たす簡略法が無いものかと調べたのですが、上手に見つけることができませんでした。[2]
仕方ないのでRamerDouglasPeucker法に少し手を加えてみることにします。
RamerDouglasPeucker法は系列データを対象としたアルゴリズムですが、今回はポリゴンに対象を限定するので、面積を求めることができます。
よって打ち切り条件を判定する点が作る三角形の面積/元となる図形の面積 < 閾値とすれば、図形の大小に関わらず同じ閾値を利用することができます。
ポリゴンの頂点数の制限は、閾値以上の点のうち面積でソートして、上限以上は切り捨てるようにします。
function _dist2(p1: number[], p2: number[]): number {
return (p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2
}
function area(linearRing: number[][]): number {
let area = 0.0;
const length = linearRing.length - 1;
for (let i = 0; i < length; i++) {
area =
area +
linearRing[i][0] * linearRing[i + 1][1] -
linearRing[i][1] * linearRing[i + 1][0];
}
return Math.abs(area / 2);
}
function _getCandidates(first: number, last: number, linearRing: number[][], threshold: number): {
index: number;
area: number;
}[] {
let memory = -Infinity;
let index = -1;
for (let i = first + 1; i < last; i++) {
const a = area([linearRing[first], linearRing[i], linearRing[last], linearRing[first]]);
if (a > memory) {
memory = a;
index = i;
}
}
if (memory < threshold) return [];
return [{
index,
area: memory
}, ..._getCandidates(first, index, linearRing, threshold),
..._getCandidates(index, last, linearRing, threshold)
]
}
function validLinearRing(points: number[][]): void {
if (
points.length <= 3 ||
points[0][0] !== points[points.length - 1][0] ||
points[0][1] !== points[points.length - 1][1]
)
throw new Error("Invalid LinearRing");
}
function rdp(
linearRing: number[][],
userOptions?: {
rate?: number;
max?: number;
}
): number[][] {
validLinearRing(linearRing);
const options = Object.assign(
{
rate: 0.001,
max: Infinity
},
userOptions
);
if (options.rate < 0 || options.rate >= 1)
throw new Error("rate must be between 0 and 1.");
if (options.max < 4)
throw new Error("max must be larger than 3.");
const threshold = area(linearRing) * options.rate;
const first = 0;
const last = linearRing.length - 1;
let memory = -Infinity;
let midpoint = -1;
for (let i = first + 1; i < last; i++) {
const d2 = _dist2(linearRing[first], linearRing[i]);
if (d2 > memory) {
memory = d2;
midpoint = i;
}
}
let candidates = [
..._getCandidates(first, midpoint, linearRing, threshold),
..._getCandidates(midpoint, last, linearRing, threshold)
];
if (options.max < linearRing.length) {
candidates.sort(function (a, b) {
return (a.area < b.area ? 1 : -1)
})
candidates = candidates.slice(0, options.max - 3);
}
const candidateIndexes = [first, last, midpoint, ...candidates.map((v) => v.index)];
return linearRing.filter((_, i) => {
return candidateIndexes.includes(i)
})
}
rate: 0.001で間引いた例
rate: 0.01で間引いた例
rate: 0.001で間引いた例
図形の大小でも形状は変わりません。
文京区の例 (「国土数値情報（行政区域データ）」（国土交通省）を加工して作成)
before 893点
rate: 0.0001 -> 156点
お手軽に済ませるには十分な結果だと思います。
面積を閾値に利用するのならVisvalingam algorithmと組み合わせたほうが素直かもしれない。
その場合のコードも載せておきます。
こちらは、隣り合う点同士で三角形を作り、その面積が最小となる点を選んで、面積が閾値以下なら点を削除します。RDPとは上から攻めるか下から攻めるかの違いです。
function vw(
linearRing: Points,
userOptions?: {
rate?: boolean;
threshold?: number;
limit?: number;
}
): Points {
validLinearRing(linearRing);
const options = Object.assign(
{
rate: false,
threshold: 0.001,
limit: Infinity
},
userOptions
);
if (options.rate) {
if (options.threshold < 0 || options.threshold >= 1)
throw new Error("when rate is true, threshold must be between 0 and 1.");
} else {
if (options.threshold < 0)
throw new Error("threshold must be positive number.");
}
if (options.limit < 4)
throw new Error("limit must be larger than 3.");
const threshold = options.rate ? area(linearRing) * options.threshold : options.threshold;
const candidates = linearRing.map((point, i) => {
if (i === 0 || i === linearRing.length - 1)
return {
score: Infinity,
index: i
}
return {
score: area([linearRing[i - 1], point, linearRing[i + 1], linearRing[i - 1]]),
index: i
}
})
while (candidates.length > 4) {
let memory = Infinity;
let eliminate = -1;
for (let i = 1; i < candidates.length - 1; i++) {
if (candidates[i].score < memory) {
memory = candidates[i].score;
eliminate = i;
}
}
if (candidates.length > options.limit || memory < threshold) {
if (eliminate - 2 >= 0) {
const p1 = candidates[eliminate - 2].index;
const p2 = candidates[eliminate - 1].index;
const p3 = candidates[eliminate + 1].index;
candidates[eliminate - 1].score = area([
linearRing[p1], linearRing[p2], linearRing[p3], linearRing[p1]
])
}
if (eliminate + 2 <= candidates.length - 1) {
const p1 = candidates[eliminate - 1].index;
const p2 = candidates[eliminate + 1].index;
const p3 = candidates[eliminate + 2].index;
candidates[eliminate + 1].score = area([
linearRing[p1], linearRing[p2], linearRing[p3], linearRing[p1]
])
}
candidates.splice(eliminate, 1)
} else {
break;
}
}
const candidateIndexes = [...candidates.map((v) => v.index)];
return linearRing.filter((_, i) => {
return candidateIndexes.includes(i)
})
}
面積の計算は平面座標系であることを前提としている。厳密にやりたければこのあたりを参考にareaメソッドを書き換えるとよいが、そこまでの精度が必要になるなら別のアルゴリズムを考えたほうがよいでしょう。
精度をあげたり、隣接するポリゴン同士が簡略後も隙間ができないようにするには、もっと別のアルゴリズムを採用する必要があります。
Discussion