pyplotのcontorfの結果を取得する

# はじめに
GISで利用するラスターデータからコンター図を得たいと思い「python使えば簡単じゃん」と思ったが以外に情報がなく手間取ったのでここに記録する。

# ラスターデータからコンター図を作る
pyplotを使えば何の不自由もなくできる。

``` python:sample1.py
import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage
from scipy.interpolate import splprep, splev

#サンプルデータ
sample_data = [
 [ 0.,  0.,  0.,  0., 0. ],
 [ 0.,  4.,  7.,  4., 0. ],
 [ 0.,  7.,  7.,  7., 0. ],
 [ 0.,  4.,  7.,  4., 0. ],
 [ 0.,  0.,  0.,  0., 0. ],
 ]

#コンターレベル
levels = [0, 5, 10]
r = plt.contourf(sample_data, levels)
plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/22962/aa541cae-76ed-3409-1b3a-e5bb8be6f101.png)

この真ん中に浮かぶ島の境界線(コンター)をエクスポートしたい。

# contourfからコンターのポリゴン情報を取得する

contourfは戻り値で[CntourSet](https://matplotlib.org/api/contour_api.html#matplotlib.contour.ContourSet)を返す。この中にコンターのポリゴンデータが入っているようだ。ダンプしてみよう。

``` python:sample2.py
import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage
from scipy.interpolate import splprep, splev

#サンプルデータ
sample_data = [
 [ 0.,  0.,  0.,  0., 0. ],
 [ 0.,  4.,  7.,  4., 0. ],
 [ 0.,  7.,  7.,  7., 0. ],
 [ 0.,  4.,  7.,  4., 0. ],
 [ 0.,  0.,  0.,  0., 0. ],
 ]

#コンターレベル
levels = [0, 5, 10]
r = plt.contourf(sample_data, levels)
for seg in r.allsegs:
    print(seg)
```

実行結果

```
[array([[1.        , 0.        ],
       [2.        , 0.        ],
       [3.        , 0.        ],
       [4.        , 0.        ],
       [4.        , 1.        ],
       [4.        , 2.        ],
       [4.        , 3.        ],
       [4.        , 4.        ],
       [3.        , 4.        ],
       [2.        , 4.        ],
       [1.        , 4.        ],
       [0.        , 4.        ],
       [0.        , 3.        ],
       [0.        , 2.        ],
       [0.        , 1.        ],
       [0.        , 0.        ],
       [1.        , 0.        ],
       [1.33333333, 1.        ],
       [1.        , 1.33333333],
       [0.71428571, 2.        ],
       [1.        , 2.66666667],
       [1.33333333, 3.        ],
       [2.        , 3.28571429],
       [2.66666667, 3.        ],
       [3.        , 2.66666667],
       [3.28571429, 2.        ],
       [3.        , 1.33333333],
       [2.66666667, 1.        ],
       [2.        , 0.71428571],
       [1.33333333, 1.        ]])]
[array([[2.        , 0.71428571],
       [2.66666667, 1.        ],
       [3.        , 1.33333333],
       [3.28571429, 2.        ],
       [3.        , 2.66666667],
       [2.66666667, 3.        ],
       [2.        , 3.28571429],
       [1.33333333, 3.        ],
       [1.        , 2.66666667],
       [0.71428571, 2.        ],
       [1.        , 1.33333333],
       [1.33333333, 1.        ],
       [2.        , 0.71428571]])]
```

list の中に2つのarrayデータが入っている。contourfでlevelsを[0,5,10]と指定しているので0:5の範囲が0番目、5:10の範囲が1番目の値のようだ。

となれば、r.allsegs[1] の値を出力すればこの件はおしまいだ。

# 2つ以上のコンターを持っている場合

r.allsegs[1] の内容は座標データarrayのarrayとなっている。てことは、5:10範囲の複数のコンターを格納しているということかな。

``` python:sample3.py

import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage
from scipy.interpolate import splprep, splev

#サンプルデータ
sample_data = [
[ 0.,  0., 0.,  0.,  0.,  0., 0. ],
[ 0.,  7., 7.,  0.,  7.,  7., 0. ],
[ 0.,  7., 7.,  0.,  7.,  7., 0. ],
[ 0.,  0., 0.,  0.,  0.,  0., 0. ],
[ 0.,  7., 7.,  0.,  7.,  7., 0. ],
[ 0.,  7., 7.,  0.,  7.,  7., 0. ],
[ 0.,  0., 0.,  0.,  0.,  0., 0. ],
]

#コンターレベル
levels = [0, 5, 10]
r = plt.contourf(sample_data, levels)
print("コンターの数:{}".format(len(r.allsegs[1])))
plt.show()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/22962/e691a9b8-3cba-1929-c577-0bb7b1a8421c.png)


```
コンターの数:4
```

allsegs[1] に4つポリゴンが存在している。複数個まとめてallsegs[1]に格納されているようだ。

# 穴あきコンターの場合

ここまでは順調。しかし、ちょっと厄介なのが **ドーナツのように穴が開いたコンター** の場合。

``` python:sample4.py
import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage
from scipy.interpolate import splprep, splev

#サンプルデータ
sample_data = [
[ 0.,  0., 0.,  0.,  0.,  0., 0. ],
[ 0.,  7., 7.,  7.,  7.,  7., 0. ],
[ 0.,  7., 0.,  0.,  0.,  7., 0. ],
[ 0.,  7., 0.,  0.,  0.,  7., 0. ],
[ 0.,  7., 0.,  0.,  0.,  7., 0. ],
[ 0.,  7., 7.,  7.,  7.,  7., 0. ],
[ 0.,  0., 0.,  0.,  0.,  0., 0. ],
]

#コンターレベル
levels = [0, 5, 10]
r = plt.contourf(sample_data, levels)
print(r.allsegs[1])
plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/22962/83e8346b-3853-4225-b559-659b7184835c.png)

```
[array([[1.        , 0.71428571],
       [2.        , 0.71428571],
       [3.        , 0.71428571],
       [4.        , 0.71428571],
       [5.        , 0.71428571],
       [5.28571429, 1.        ],
       [5.28571429, 2.        ],
       [5.28571429, 3.        ],
       [5.28571429, 4.        ],
       [5.28571429, 5.        ],
       [5.        , 5.28571429],
       [4.        , 5.28571429],
       [3.        , 5.28571429],
       [2.        , 5.28571429],
       [1.        , 5.28571429],
       [0.71428571, 5.        ],
       [0.71428571, 4.        ],
       [0.71428571, 3.        ],
       [0.71428571, 2.        ],
       [0.71428571, 1.        ],
       [1.        , 0.71428571],
       [1.28571429, 2.        ],
       [1.28571429, 3.        ],
       [1.28571429, 4.        ],
       [2.        , 4.71428571],
       [3.        , 4.71428571],
       [4.        , 4.71428571],
       [4.71428571, 4.        ],
       [4.71428571, 3.        ],
       [4.71428571, 2.        ],
       [4.        , 1.28571429],
       [3.        , 1.28571429],
       [2.        , 1.28571429],
       [1.28571429, 2.        ]])]
```

穴あきポリゴンなので外側の境界線と穴の境界線が存在するはずだが、穴がない場合と構造上の違いがないように見える。
実はこれ、外側の境界線が閉じた後に内側の境界線の情報が続いているのだ。
0番目の[1., 0.71428571]から21番目の[1., 0.71428571]が外側の境界線で、内側はそれ以降のデータとなる。
つまり、外側の境界線が閉じる(始点と同じ値)となったら、後続は内側の境界線となる。
内側の穴が複数ある場合は、1つの内側の境界線が閉じたら後続が内側の境界線となる。
また、外側の境界線は反時計回り、内側の境界線は時計回り順で格納されている。

# 複数コンターと穴あきを考慮してコンターの境界を出力

そんなわけで「コンターの外側と内側同士はポリゴンが閉じたところで分割する」というルールにのっとって修正してみる。

```python:sample5.py
import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage
from scipy.interpolate import splprep, splev

#サンプルデータ
sample_data = [
[ 0.,  0., 0., 0., 0., 0., 0., 0. ],
[ 0.,  7., 7., 0., 7., 7., 7., 0. ],
[ 0.,  7., 7., 0., 7., 0., 7., 0. ],
[ 0.,  7., 7., 0., 7., 7., 7., 0. ],
[ 0.,  7., 7., 0., 7., 0., 7., 0. ],
[ 0.,  7., 7., 0., 7., 7., 7., 0. ],
[ 0.,  0., 0., 0., 0., 0., 0., 0. ],
]
EPS = 1e-20
#穴を含むコンターを複数のリストに分割する
def normalize(pts):
    results = []
    points = []
    sp = None
    for p in pts:
        points.append(p)
        if sp is not None:
            if abs(sp[0] - p[0]) < EPS and abs(sp[1] - p[1]) < EPS:
                results.append(np.array(points))
                points = []
                sp = None
        else:
            sp = p

    return results

#コンターレベル
levels = [0, 5, 10]
r = plt.contourf(sample_data, levels)
for i, c in enumerate(r.allsegs[1]):
    d = normalize(c)
    print("外側境界({})".format(i))
    print(d[0])
    if 1 < len(d):
        print("内側境界({})".format(i))
        print(d[1:])
plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/22962/868b536f-d9b7-3e38-906c-f34806dfff14.png)

```
外側境界(0)
[[1.         0.71428571]
 [2.         0.71428571]
 [2.28571429 1.        ]
 [2.28571429 2.        ]
 [2.28571429 3.        ]
 [2.28571429 4.        ]
 [2.28571429 5.        ]
 [2.         5.28571429]
 [1.         5.28571429]
 [0.71428571 5.        ]
 [0.71428571 4.        ]
 [0.71428571 3.        ]
 [0.71428571 2.        ]
 [0.71428571 1.        ]
 [1.         0.71428571]]
外側境界(1)
[[4.         0.71428571]
 [5.         0.71428571]
 [6.         0.71428571]
 [6.28571429 1.        ]
 [6.28571429 2.        ]
 [6.28571429 3.        ]
 [6.28571429 4.        ]
 [6.28571429 5.        ]
 [6.         5.28571429]
 [5.         5.28571429]
 [4.         5.28571429]
 [3.71428571 5.        ]
 [3.71428571 4.        ]
 [3.71428571 3.        ]
 [3.71428571 2.        ]
 [3.71428571 1.        ]
 [4.         0.71428571]]
内側境界(1)
[array([[4.28571429, 2.        ],
       [5.        , 2.71428571],
       [5.71428571, 2.        ],
       [5.        , 1.28571429],
       [4.28571429, 2.        ]]), array([[4.28571429, 4.        ],
       [5.        , 4.71428571],
       [5.71428571, 4.        ],
       [5.        , 3.28571429],
       [4.28571429, 4.        ]])]
```

これで、穴が開いたコンターでも対応可能となった。

# B-spline補間する

この手のデータは「いい感じに出してよ」といわれるので定番の B-spline で補完しよう。

```python:sample6.py
import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage
from scipy.interpolate import splprep, splev

#サンプルデータ
sample_data = [
[ 0.,  0., 0., 0., 0., 0., 0., 0. ],
[ 0.,  7., 7., 0., 7., 7., 7., 0. ],
[ 0.,  7., 7., 0., 7., 0., 7., 0. ],
[ 0.,  7., 7., 0., 7., 7., 7., 0. ],
[ 0.,  7., 7., 0., 7., 0., 7., 0. ],
[ 0.,  7., 7., 0., 7., 7., 7., 0. ],
[ 0.,  0., 0., 0., 0., 0., 0., 0. ],
]
EPS = 1e-20
#穴を含むコンターを複数のリストに分割する
def normalize(pts):
    results = []
    points = []
    sp = None
    for p in pts:
        points.append(p)
        if sp is not None:
            if abs(sp[0] - p[0]) < EPS and abs(sp[1] - p[1]) < EPS:
                results.append(np.array(points))
                points = []
                sp = None
        else:
            sp = p

    return results

#コンターレベル
levels = [0, 5, 10]
r = plt.contourf(sample_data, levels)
for i, c in enumerate(r.allsegs[1]):
    d = normalize(c)
    is_innser = False
    for pp in d:
        tck, u = splprep(pp.T, u=None, s=0, per=1) 
        u_new = np.linspace(u.min(), u.max(), 100)
        x_new, y_new = splev(u_new, tck, der=0)
        plt.plot(x_new, y_new, 'r-')
        if not is_innser:
            print("外側境界({})".format(i))
        else:
            print("内部境界({})".format(i))
        print(x_new,y_new)
        is_innser = True
plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/22962/461cbf9c-87bd-8704-84fb-8bf9527cf457.png)

…なんか長い直線が変な感じに補完されてる。
これは、contourfが直線を整数位置の座標で分割をしてしまうためだ。なので、「傾きが同じだった点は捨てる」処理を入れると期待通りになる。

``` python:sample7.py
import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage
from scipy.interpolate import splprep, splev

#サンプルデータ
sample_data = [
[ 0.,  0., 0., 0., 0., 0., 0., 0. ],
[ 0.,  7., 7., 0., 7., 7., 7., 0. ],
[ 0.,  7., 7., 0., 7., 0., 7., 0. ],
[ 0.,  7., 7., 0., 7., 7., 7., 0. ],
[ 0.,  7., 7., 0., 7., 0., 7., 0. ],
[ 0.,  7., 7., 0., 7., 7., 7., 0. ],
[ 0.,  0., 0., 0., 0., 0., 0., 0. ],
]
EPS = 1e-20
#傾き算出
def get_a(p1, p2):
    d0 = p2[0] - p1[0]
    d1 = p2[1] - p1[1]
    if d0 == 0:
        return 0
    return d1 / d0
#穴を含むコンターを複数のリストに分割する
def normalize(pts):
    results = []
    points = []
    sp = None
    pp = None
    prev_a = None
    for p in pts:
        if pp is not None:
            a = get_a(pp, p)
            if prev_a is not None and abs(prev_a - a) < EPS:
                points.pop()
            prev_a = a
        pp = p
        points.append(p)
        if sp is not None:
            if abs(sp[0] - p[0]) < EPS and abs(sp[1] - p[1]) < EPS:
                results.append(np.array(points))
                points = []
                sp = None
        else:
            sp = p

    return results

#コンターレベル
levels = [0, 5, 10]
r = plt.contourf(sample_data, levels)
for i, c in enumerate(r.allsegs[1]):
    d = normalize(c)
    is_innser = False
    for pp in d:
        tck, u = splprep(pp.T, u=None, s=0, per=1) 
        u_new = np.linspace(u.min(), u.max(), 100)
        x_new, y_new = splev(u_new, tck, der=0)
        plt.plot(x_new, y_new, 'r-')
        if not is_innser:
            print("外側境界({})".format(i))
        else:
            print("内部境界({})".format(i))
        print(x_new,y_new)
        is_innser = True
plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/22962/effed9af-608b-c5c3-bab8-d0c50bc33ff8.png)

意図したとおりに補間された。

# まとめ
pyplot.contourf の結果をB-spline補間して座標値リストとして取得することができた。外部システムとのI/Fをきちんと作ればラスターデータから任意の閾値のコンターを作成し、B-spline補間した結果を受け取ることが可能となる。

強力で生産性の高いpyplotを外部システムから利用する価値は大きい。しかし、contourfの結果構造がなかなか理解が難しかったので記事とした。