NumPy 2D ndarray の、NaN だらけ外周を取り除く問題

かっこつけて英語のタイトルにしようとして果てる、の巻。

英語のタイトル付けることが多いのは、元ネタが StackOverflow のことが多いからなの。そいで、StackOverflow からのネタが多いのは、これとかこれが理由。

つーか、日本語でも伝わりにくいなぁ…。立体地形図を作るために色んな領域を扱ってるとね、「海なんかいらないのに!」となることが多いのね。陸域だけに興味があるわけだからさ。例えば伊豆半島南端ぴったりを南端に、銚子をぴったり東端に、とするのに、地理院地図WEBや Google Map で調べるのは、なかなかに苦痛だ。

海データは NaN にしている、ワタシは。つまり、入力のこれ:

 1 >>> import numpy as np
 2 >>> a = np.full((10, 11), np.nan)
 3 >>> a = np.triu(a) + np.arange(10*11).reshape((10, 11))
 4 >>> a[-2:,:2] = np.nan
 5 >>> a[-4:-2,2:4] = np.nan
 6 >>> a
 7 array([[  nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan],
 8        [  11.,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan],
 9        [  22.,   23.,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan],
10        [  33.,   34.,   35.,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan],
11        [  44.,   45.,   46.,   47.,   nan,   nan,   nan,   nan,   nan,   nan,   nan],
12        [  55.,   56.,   57.,   58.,   59.,   nan,   nan,   nan,   nan,   nan,   nan],
13        [  66.,   67.,   nan,   nan,   70.,   71.,   nan,   nan,   nan,   nan,   nan],
14        [  77.,   78.,   nan,   nan,   81.,   82.,   83.,   nan,   nan,   nan,   nan],
15        [  nan,   nan,   90.,   91.,   92.,   93.,   94.,   95.,   nan,   nan,   nan],
16        [  nan,   nan,  101.,  102.,  103.,  104.,  105.,  106.,  107.,   nan,   nan]])

から、これを抽出したいわけである:

 1 >>> # 何か操作…
 2 >>> a
 3 array([[  11.,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan],
 4        [  22.,   23.,   nan,   nan,   nan,   nan,   nan,   nan,   nan],
 5        [  33.,   34.,   35.,   nan,   nan,   nan,   nan,   nan,   nan],
 6        [  44.,   45.,   46.,   47.,   nan,   nan,   nan,   nan,   nan],
 7        [  55.,   56.,   57.,   58.,   59.,   nan,   nan,   nan,   nan],
 8        [  66.,   67.,   nan,   nan,   70.,   71.,   nan,   nan,   nan],
 9        [  77.,   78.,   nan,   nan,   81.,   82.,   83.,   nan,   nan],
10        [  nan,   nan,   90.,   91.,   92.,   93.,   94.,   95.,   nan],
11        [  nan,   nan,  101.,  102.,  103.,  104.,  105.,  106.,  107.]])

ググってもピンポイントではヒットせんなぁ。

仕方ないので自力でみたが。もうね、一撃でとかそんなん無理みたいで。まっとうにやるしかなかったわ:

 1 >>> def get_bound(arr2d, f=np.isnan):
 2 ...      #
 3 ...      def _sl(arr2d, i, axis):
 4 ...           if axis == 0:
 5 ...                return arr2d[i,:]
 6 ...           else:
 7 ...                return arr2d[:,i]
 8 ...      #
 9 ...      def _gb(arr2d, axis, f):
10 ...           idx_l = 0
11 ...           for i in range(a.shape[axis]):
12 ...                idx_l = i
13 ...                if not np.all(f(_sl(a, i, axis))):
14 ...                     break
15 ...           idx_h = a.shape[axis]
16 ...           for i in range(a.shape[axis], 1, -1):
17 ...                idx_h = i
18 ...                if not np.all(f(_sl(a, i - 1, axis))):
19 ...                     break
20 ...           return slice(idx_l, idx_h)
21 ...      #
22 ...      return _gb(arr2d, 0, f), _gb(arr2d, 1, f)
23 ...
24 >>> get_bound(a)
25 (slice(1, 10, None), slice(0, 9, None))
26 >>> a[get_bound(a)]
27 array([[  11.,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan],
28        [  22.,   23.,   nan,   nan,   nan,   nan,   nan,   nan,   nan],
29        [  33.,   34.,   35.,   nan,   nan,   nan,   nan,   nan,   nan],
30        [  44.,   45.,   46.,   47.,   nan,   nan,   nan,   nan,   nan],
31        [  55.,   56.,   57.,   58.,   59.,   nan,   nan,   nan,   nan],
32        [  66.,   67.,   nan,   nan,   70.,   71.,   nan,   nan,   nan],
33        [  77.,   78.,   nan,   nan,   81.,   82.,   83.,   nan,   nan],
34        [  nan,   nan,   90.,   91.,   92.,   93.,   94.,   95.,   nan],
35        [  nan,   nan,  101.,  102.,  103.,  104.,  105.,  106.,  107.]])
36 >>> 
37 >>> # ややエレガントにみえるバージョン
38 >>> # (ndimsが1でもいける(確認済み)。3でもいけるはず? (3は未確認。))
39 >>> def get_bound(arr, f=np.isnan):
40 ...     #
41 ...     def _gb(arr, axis, f):
42 ...         _sl = [slice(None), ] * arr.ndim
43 ...         idx_l = 0
44 ...         for i in range(arr.shape[axis]):
45 ...             idx_l = i
46 ...             _sl[axis] = slice(i, i + 1)
47 ...             if not np.all(f(arr[_sl])):
48 ...                 break
49 ...         idx_h = arr.shape[axis]
50 ...         for i in range(arr.shape[axis], 1, -1):
51 ...             idx_h = i
52 ...             _sl[axis] = slice(i - 1, i)
53 ...             if not np.all(f(arr[_sl])):
54 ...                 break
55 ...         return slice(idx_l, idx_h)
56 ...     #
57 ...     return [_gb(arr, axis, f) for axis in range(arr.ndim)]
58 ...
59 >>> get_bound(a)
60 [slice(1, 10, None), slice(0, 9, None)]
61 >>> a[get_bound(a)]
62 array([[  11.,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan],
63        [  22.,   23.,   nan,   nan,   nan,   nan,   nan,   nan,   nan],
64        [  33.,   34.,   35.,   nan,   nan,   nan,   nan,   nan,   nan],
65        [  44.,   45.,   46.,   47.,   nan,   nan,   nan,   nan,   nan],
66        [  55.,   56.,   57.,   58.,   59.,   nan,   nan,   nan,   nan],
67        [  66.,   67.,   nan,   nan,   70.,   71.,   nan,   nan,   nan],
68        [  77.,   78.,   nan,   nan,   81.,   82.,   83.,   nan,   nan],
69        [  nan,   nan,   90.,   91.,   92.,   93.,   94.,   95.,   nan],
70        [  nan,   nan,  101.,  102.,  103.,  104.,  105.,  106.,  107.]])

こんなん、かなり一般的なタスクなんじゃないかって気がするんだけどなぁ…。例えば画像処理でまわりの黒い部分を取り除きたい、とかさ。

なお、「自力た」言うても、これはヒントにした。

ほいだば、実践だす。

ラフに選択した範囲ではこうなるとして:

1 # elevs はどこかからもってきた標高データ (上の絵)
2 elevs = elevs[get_bound(elevs)]  # 海ばっかのとこはいらない

うむ。よか。


ちょっとだけバリエーション。elevs.T (転置) と np.fliplr, np.flipud だけで欠けの方向のバリエーションをすぐに作れるのでね。










島だろうとヘコたれない。