凹型船体
K-最近傍アプローチを使用したクラスター境界の作成
のダニエル・イアンによる「灰色の空の下の緑の草地の白いボート」数ヶ月前、私はここに英国の交通事故のホットスポットのマッピングに関する記事をMediumに書きました。 私は、地理的データに対するDBSCANクラスタリングアルゴリズムの使用を説明することについて主に心配していました。 記事では、報告された交通事故に関する英国政府が発行した地理情報を使用しました。 私の目的は、交通事故が最も頻繁に報告されている地域を見つけるために、密度ベースのクラスタリングプロセスを実行することでした。 最終的な結果は、これらの事故のホットスポットを表すジオフェンスのセットの作成でした。
特定のクラスター内のすべてのポイントを収集することで、クラスターがマップ上でどのように見えるかを知ることができますが、クラスターの外形という重要な情報が不足しています。 この場合、地図上でジオフェンスとして表現できる閉じた多角形について話しています。 ジオフェンス内の任意のポイントは、この形状を興味深い情報にするクラスタに属していると仮定することができます。 ポリゴンの内側にあるすべての新しくサンプリングされたポイントは、対応するクラスターに属していると仮定できます。 私が記事で示唆したように、あなたはあなた自身のサンプリングされたGPS位置を分類するためにそれらを使用することによって、あなたの運転リ
ここでの問題は、特定のクラスターを構成する点の雲から意味のある多角形を作成する方法です。 最初の記事での私のアプローチはやや素朴で、本番コードですでに使用していた解決策を反映していました。 この解は、クラスタの各点を中心に円を配置し、すべての円をマージして雲の形をした多角形を形成することを必要としました。 結果は非常に素晴らしい、また現実的ではありません。 また、最終的な多角形を構築するための基本形状として円を使用することにより、これらはより合理化された形状よりも多くのポイントを持ち、そ
一方、このアプローチは、Shapelyのcascaded_union関数を使用してすべての円をマージするため、（少なくとも開発者の観点から）計算が単純であるとい もう1つの利点は、クラスター内のすべてのポイントを使用してポリゴンの形状が暗黙的に定義されることです。
より洗練されたアプローチのためには、クラスターの境界点、つまり点群の形状を定義するように見える点を何とか特定する必要があります。 興味深いことに、いくつかのDBSCAN実装では、クラスタリングプロセスの副産物として実際に境界点を回復することができます。 残念ながら、この情報は（明らかに）SciKit Learnの実装では利用できないため、実行する必要があります。
最初に頭に浮かんだアプローチは、点の集合の凸包を計算することでした。 これはよく理解されたアルゴリズムですが、このような凹面形状を処理しないという問題に苦しんでいます:
この形状は、基礎となる点の本質を正しくキャプチャしません。 それが識別器として使用された場合、いくつかのポイントは、そうでないときにクラスター内にあると誤って分類されます。 別のアプローチが必要だ
凹型船体の代替
幸いなことに、この状態に代わるものがあります：凹型船体を計算することができます。 前の画像と同じ点のセットに適用したときの凹型の船体の外観は次のとおりです:
または多分これ:
凸面の船体とは対照的に、一連の点の凹面の船体が何であるかの単一の定義はありません。 私がここで提示しているアルゴリズムでは、船体をどのように凹ませたいかの選択は、単一のパラメータ：k—船体計算中に考慮される最近傍の数。 これがどのように機能するか見てみましょう。
アルゴリズム
私がここで提示しているアルゴリズムは、ポルトガルのミンホ大学のAdriano MoreiraとMaribel Yasmina Santosによって十年以上前に記述されました。 要約から:
与えられた点によって占有される面積を表す非凸殻上に凸を生成する平面内の点の集合の包絡線を計算するアルゴリズムについて述べた。 提案したアルゴリズムは，唯一のアルゴリズムパラメータであるｋの値を用いて最終解の”滑らかさ”を制御するｋ-最近傍法に基づいている。
このアルゴリズムを地理情報に適用するため、角度や距離を計算するときにいくつかの変更を加えなければなりませんでした。 しかし、これらはアルゴリズムの要点を変更するものではなく、次のステップで広く説明することができます:
- 最も低いy（緯度）座標を持つ点を見つけて、それを現在のものにします。
- 現在の点に最も近いk点を求めます。
- k個の最も近い点から、前の角度から最大の右回転に対応する点を選択します。 ここでは、ベアリングの概念を使用し、270度（真西）の角度から始めます。
- 成長している線の文字列に新しい点を追加することによって、それ自体が交差しないかどうかを確認します。 そうである場合は、k-最も近い点から別の点を選択するか、より大きな値kで再起動します。
- 新しい点を現在の点にして、リストから削除します。
- k回の反復の後、最初の点をリストに戻します。
- は2番にループします。
アルゴリズムは非常に単純なようですが、特に地理座標を扱っているため、多くの詳細に注意する必要があります。 距離と角度は別の方法で測定されます。
私が公開しているコード
は、前の記事のコードの適応バージョンです。 同じクラスタリングコードと同じ雲の形のクラスタージェネレータが見つかります。 更新されたバージョンには、geomath.hullsという名前のパッケージが含まれており、ConcaveHullクラスを見つけることができます。 あなたの凹面の外皮を作成するには、次のようにしてください:
上記のコードでは、pointsは次元(N,2)の配列であり、行には観測点が含まれ、列には地理座標(経度、緯度)が含まれています。 結果の配列はまったく同じ構造を持ちますが、クラスターのポリゴンシェイプに属するポイントのみが含まれます。 種類のフィルター。
私たちは配列を処理するので、numpyを争いに持ち込むのは自然なことです。 すべての計算は可能な限り正当にベクトル化され、配列から項目を追加および削除するときのパフォーマンスを向上させる努力がなされました（スポイラー： 不足している改善点の1つは、コードの並列化です。 しかし、それは待つことができます。
翻訳中にいくつかの最適化が行われましたが、論文で公開されているようにアルゴリズムの周りにコードを整理しました。 このアルゴリズムは、論文によって明確に識別されるいくつかのサブルーチンを中心に構築されているので、今すぐそれらを邪魔してみましょう。 あなたの読書の快適さのために、私は紙で使用されているのと同じ名前を使用します。
CleanList-ポイントのリストのクリーニングは、クラスコンストラクタで実行されます:
ご覧のとおり、ポイントのリストは、パフォーマンス上の理由からNumPy配列として実装されています。 リストのクリーニングは、10行目で実行され、一意のポイントのみが保持されます。 データセット配列は、行の観測値と2つの列の地理座標で編成されます。 私はまた、メインデータセット配列にインデックスを付けるために使用される13行目にブール配列を作成していることに注意してください。 私はNumPyのドキュメントで”マスク”と呼ばれるこの技術を見てきましたが、それは非常に強力です。 素数については、後で説明します。
FindMinYPoint—これには小さな関数が必要です:
この関数は、dataset配列を引数として呼び出され、緯度が最も低い点のインデックスを返します。 行は、最初の列に経度、2番目の列に緯度でエンコードされていることに注意してください。
RemovePoint
AddPoint—これらはindices配列を使用するため、非常に簡単です。 この配列は、アクティブなインデックスをメインデータセット配列に格納するために使用されるため、データセットから項目を削除するのは簡単です。
論文で説明されているアルゴリズムは、船体を構成する配列に点を追加する必要がありますが、これは実際には次のように実装されています:
後で、行の文字列が交差していないと見なされると、test_hull変数はhullに割り当てられます。 しかし、私はここで試合に先んじています。 データセット配列からポイントを削除するのは、次のように簡単です:
self.indices = False
それを元に戻すことは、同じインデックスの配列値をtrueに戻すことだけです。 しかし、このすべての利便性は、インデックスに私たちのタブを維持することの価格が付属しています。 これについては後で詳しく説明します。
NearestPoints—平面座標を扱っていないので、ここで物事が面白くなり始めるので、Pythagorasで、Haversineで出てください:
2番目と3番目のパラメーターはデータセット形式の配列で、1番目の列の経度と2番目の列の緯度です。 ご覧のとおり、この関数は、2番目の引数の点と3番目の引数の点との間の距離の配列をメートル単位で返します。 これらが得られれば、k-最近傍を簡単に得ることができます。 しかし、そのための特殊な機能があり、それはいくつかの説明に値する:
この関数は、基本インデックスを持つ配列を作成することから始めます。 これらは、データセット配列から削除されていない点のインデックスです。 たとえば、10ポイントのクラスターで最初のポイントを削除することで開始した場合、基本インデックス配列は次のようになります。 次に、距離を計算し、結果の配列インデックスを並べ替えます。 最初のkが抽出され、ベースインデックスを取得するためのマスクとして使用されます。 それは一種のねじれだが、動作します。 ご覧のとおり、この関数は座標の配列を返しませんが、データセット配列へのインデックスの配列を返します。
SortByAngle—単純な角度を計算しているのではなく、ベアリングを計算しているので、ここではさらに問題があります。 これらは、角度が時計回りに増加して、真北ゼロ度として測定されます。 ベアリングを計算するコードのコアは次のとおりです:
この関数は、最初の引数にインデックスがある点から、3番目の引数にインデックスがある点まで測定されたベアリングの配列を返します。 ソートは簡単です:
この時点で、candidates配列には、方位の降順でソートされたk-最も近い点のインデックスが含まれています。
IntersectQ—私自身の線交点関数を転がすのではなく、私は助けを求めてShapelyに目を向けました。 実際、ポリゴンを構築する際には、基本的にライン文字列を処理し、前のものと交差しないセグメントを追加します。 このためのテストは簡単です：私たちは、建設中の船体配列をピックアップし、それを形の良い線の文字列オブジェクトに変換し、それが単純（非自己交差）
一言で言えば、整形された行の文字列は自己交差すると複雑になるため、is_simple述語はfalseになります。 落ち着け
PointInPolygon—これは実装するのが最も難しいことが判明しました。 最終的なハルポリゴン検証を実行するコードを見て説明してください（クラスターのすべてのポイントがポリゴンに含まれているかどうかを確認し):
交差と包含をテストするShapelyの関数は、最終的なハルポリゴンがクラスターのすべての点と重なるかどうかを確認するのに十分であったはずですが、そうではありませんでした。 どうして？ Shapelyは座標に依存しないため、緯度と経度で表現された地理座標をデカルト平面上の座標とまったく同じ方法で処理します。 しかし、あなたが球上に住んでいて、角度（またはベアリング）が測地線に沿って一定ではないとき、世界は異なった振る舞いをします。 バグダッドと大阪を結ぶ測地線の参照の例は、これを完全に示しています。 いくつかの状況下では、アルゴリズムは方位基準に基づいて点を含めることができますが、後でShapelyの平面アルゴリズムを使用して、ポリゴンのわずか それが小さな距離補正がそこでやっていることです。
これを理解するのにしばらく時間がかかりました。 私のデバッグヘルプは、QGIS、フリーソフトウェアの素晴らしい作品でした。 疑わしい計算のすべてのステップで、データをwkt形式でcsvファイルに出力して、レイヤーとして読み込むことにしました。 本当の命の恩人！
最後に、ポリゴンがクラスタのすべてのポイントをカバーできない場合、唯一のオプションはkを増やして再試行することです。 ここで私は自分の直感を少し追加しました。
Prime k
この記事では、kの値を1ずつ増やし、アルゴリズムをゼロから再度実行することを示唆しています。 このオプションを使用した初期のテストはあまり満足のいくものではありませんでした。 これはkの遅い増加によるものだったので、私は別の増加スケジュールを使用することにしました：素数のテーブル。 アルゴリズムはすでにk=3から始まっているので、素数のリストで進化させるのは簡単な拡張でした。 これはあなたが再帰呼び出しで起こっているのを見るものです:
私は素数のためのものを持っています、あなたが知っている…
Blow Up
このアルゴリズムによって生成された凹型の船体ポリゴンは、船体の内側の点のみを判別するため、さらに処理する必要がありますが、それには近くありません。 解決策は、これらのスキニークラスターにいくつかのパディングを追加することです。 ここで私は前に使用されたのとまったく同じ手法を使用しています、そしてここでそれがどのように見えるかです:
ここでは、Shapelyのbuffer関数を使用してトリックを行いました。
この関数は形状のある多角形を受け入れ、それ自身の膨張したバージョンを返します。 第二のパラメータは、追加されたパディングのメートル単位の半径です。
コードの実行
GitHubリポジトリからコードをローカルマシンにプルすることから始めます。 実行するファイルは、メインディレクトリのShowHotSpots.pyです。 最初の実行時に、コードは2013年から2016年までの英国の交通事故データを読み取り、それをクラスター化します。 結果は、その後の実行のためにCSVファイルとしてキャッシュされます。
次に、2つのマップが表示されます:1つ目は雲の形をしたクラスターを使用して生成され、2つ目はここで説明する凹状クラスタリングアルゴリズムを使用します。 ポリゴン生成コードの実行中に、いくつかの障害が報告されることがあります。 アルゴリズムが凹型のハルを作成できない理由を理解するために、コードはクラスターをCSVファイルにdata/out/failed/ディレクトリに書き込みます。 いつものように、QGISを使用してこれらのファイルをレイヤーとしてインポートできます。
本質的に、このアルゴリズムは、自己交差せずに形状を”一周”するのに十分な点が見つからない場合に失敗します。 これは、これらのクラスターを破棄するか、またはそれらに別の処理（凸包または合体気泡）を適用する準備ができている必要があることを意味します。
この記事では，ＤＢＳＣＡＮが生成した地理的クラスタを凹形状に後処理する方法を示した。 この方法は、他の選択肢と比較して、クラスタのためのより良いフィッティング外部ポリゴンを提供することができます。
読み、コードをいじり楽しむためにありがとう！
Kryszkiewicz M.,Lasek P.(2010)TI-DBSCAN:三角不等式によるDBSCANによるクラスタリング。 In:Szczuka M.,Kryszkiewicz M.,Ramanna S.,Jensen R.,Hu Q.(eds)Rough Sets and Current Trends in Computing. 2010年、RSCTCに移籍した。 コンピュータサイエンスの講義ノート、vol6086。 Springer,Berlin,Heidelberg
Scikit-learn:Pythonでの機械学習,Pedregosa et al.,JMLR12,pp. 2825-2830,2011
Moreira,A.and Santos,M.Y.,2007,Concave Hull:a k-nearest neighbors approach for the computation by a region of a set of points
緯度/経度ポイント間の距離、方位などを計算する
GitHub repository
