GIS
foss4g
Turf.js
Maplat
mapshaper
FOSS4G Day 15

mapshaperでトポロジーエラー検出

Maplatのkochizufanです。
去年まではその2もあったのに今年はその1もなかなか埋まらない感じなので、空気読まず2トピックめ書いてみます。
もし「俺が書くつもりだったのに!」な人がいたらごめんなさい。

内容的に前回の続きです。
前回のturf.js同様、Maplatで小便利に使ってるけど多分あまり知られていないプロダクトを紹介します。

turf.jsのおかげで古地図/正確地図対応点群の三角形分割が実現でき、古地図/正確地図間の座標変換ができるようになったMaplatですが、もし三角形分割のポリゴンにトポロジーエラーが存在すると、せっかくの売りである座標系全体の1対1変換が実現できません。
トポロジーエラーというとわからない人もいるかもしれませんが、こういう感じのものです。
トポロジーエラー
本来三角網は重なりなく隙間なく領域を満たさないといけないのですが、上の画像だと三角形の辺が重なってしまっているのが見て取れます。
こういうエラーがあるとうまく座標変換できないので、検知して修正を促してやらないといけないのですが、これがturf.jsだけだとうまくできない。

turf.jsにも重なりを検知するturf.kinksという関数はあるのですが、残念ながらこいつは単独のポリゴンの中での捻れしか検知できません。
マルチポリゴンの中で複数ポリゴン間の辺の重なり合いを検知することはしてくれないのです。

で、それをやってくれるライブラリを探していて見つけたのが、mapshaperというライブラリでした。
このライブラリ、メインの機能はShapefile、GeoJson、TopoJson等のGISデータフォーマットを相互変換してくれるライブラリです。
が、なんでかよくわかりませんが、その変換ライブラリにおまけ?としてジオメトリ内の辺の交わりを検出してくれるサブ機能がついてくるのです。
mapshaper Web GUI
上の画像はmapshaperのWeb GUIでのデータ読み込みダイアログですが、オプションに「detect line intersections」というのが見えます。
これをオンにしてデータを読み込ませると、辺の交わりを検出して画面上にエラーとして出力してくれます。

が、おまけ機能だけあって?なかなか十分な情報がない。
Web GUIに機能があるということは、ライブラリにも機能が実装されているはずですが、どう関数を呼べばいいのかよくわからない。
Web GUIのソースコードとにらめっこしたり、どうしてもわからなくて作者にGithubのissueで質問したりしてようやく、ライブラリからの呼び出し方がわかりました。

トポロジーエラー検出手順
// 辺の交わり検出関数はmapshaper.internalという属性の下にあります
var internal = mapshaper.internal; 
// 点群から三角網マルチポリゴンのGeoJsonを生成したとします。
var tins = turf.tin(points); 

// mapshaperは直接GeoJsonを扱えないので、マルチポリゴンの頂点だけを配列にします。
var coords = tins.features.map(function(poly) { return poly.geometry.coordinates[0]; });

// 全然意味はわからないんですが、下記の順番で呼ぶと交点の配列が生成されます。
// ただし、個々の辺からみた交点の一覧なので、同じ点が重複して何度か含まれます。
var arcs = new internal.ArcCollection(coords);
var xy = internal.findSegmentIntersections(arcs);
var retXy = internal.dedupIntersections(xy);

// 重複の解消と、ついでにマルチポイントのGeoJsonに変換して結果を得ます。
var retGeoJson = retXy.reduce(function(prev, point, index, array) {
    if (!prev) prev = {};
    prev[point.x + ':' + point.y] = point;
    if (index != array.length - 1) return prev;
    return Object.keys(prev).map(function(key) {
        return turf.point([prev[key].x, prev[key].y]);
    });
}, []);

この手順で呼び出すと、mapshaperを使ってトポロジーエラーの検出と、存在する場合その辺の交点座標を得ることができます。
トポロジーエラーを検出しないといけないユースケースをお持ちの方がどれほどおられるかわかりませんが、もし機会があれば使ってみてください。

あと最後のおまけとして...
mapshaperはメイン機能がGISデータ形式変換なので、それなりにライブラリサイズが大きくトポロジーエラー検出のためだけに使うにはちょっとライブラリサイズがヘビーなのですが、私の方でライブラリサイズを削減するために、トポロジーエラー検出に必要なコードだけをピックアップした縮小版を作っています。
あらゆる例外を網羅できてるかはなはだ怪しいところがある代物ですが、自己責任で使いたい人がおられたら使ってみてください。