Python
Amazon
画像処理
DeepLearning
Kaggle

Kaggle まとめ: Planet, Understanding the Amazon from Space

More than 1 year has passed since last update.

1. はじめに

過去に参加したKaggleの情報をアップしていきます.
ここでは,Planet: Understanding the Amazon from Spaceのデータ紹介とフォーラムでの目立った議論をピックアップします.
優勝者のアプローチなどの紹介は,別の記事で紹介します.

2. 背景

Screen Shot 2017-07-17 at 21.44.54.png

アマゾン盆地では大量の森林が日々失われています.これにより世界的な生態系の破壊,居住地の喪失,気象変動,その他多くの壊滅的影響が発生しています.森林での自然破壊や人間の侵入をデータ化することで,政府や地域のステークホルダーが迅速かつ高価的にこのような問題に反応することができるようになります.
今回のコンペでは,3-5米での解像度を持つ衛星写真を用いることで,小規模な森林破壊や人的・偶発的といった判断が可能となります.
そこで,今回のコンペでは,衛星写真から地表面のラベリングを行います.

今回の特徴的な点としては,以下が挙げられます.

  • 衛星写真
  • 多クラスラベリング

3. 評価指標

今回の評価指標はF2スコアになります.
スコアリングには一般的にF scoreを使用しますが,F2スコアはPrecisionよりRecallの重みを大きくした指標です.
Screen Shot 2017-07-17 at 21.54.54.png

4. データの紹介

今回のデータは,アマゾン(ブラジル,ペルー,ウルグアイ,コロンビア,ベネズエラ,ガイアナ,ボリビア,エクアドル)の衛星写真です.
各写真には3グループのラベル(1. 気象状態,2. 一般的な地表,3. 珍しい地表)が一つ以上ついています.

chips.jpg

以下では各ラベルを紹介していきます.

4.1. 曇りラベル

衛星写真では雲が主要な障壁となります.ここでは3通りの曇りラベルが存在します.

4.1.1. Cloudy Scenes (曇り)

cloudy_1.jpg

4.1.2. Partly Cloudy Scenes (部分的な曇り)

pc1.jpg

4.1.3. Hazy Scenes (うす曇り)

haze1.jpg

4.2. 一般的な地理ラベル

7つの一般的なラベルを紹介します.

4.2.1. 熱帯雨林

データの中で最も多いラベルです.
primary.jpg

4.2.2. 水 (川,湖)

アマゾン盆地で最も重要な特徴となる,川,貯水池,三日月湖を示すラベル.

river.jpg

4.2.3. 居住地

居住地やビルを示すラベル.人口密集地から地方の村まで含む.小さく一人暮らしのスポット的な居住区ほどラベリングが難しくなる.

habitation1.jpg

4.2.4. 農耕地

アマゾンでは商用目的の作物を栽培する農耕地をラベリングすることは重要な技術となる.

agg1.jpg

4.2.5. 道路

road.jpg

4.2.6. 耕作地

移動耕作は農耕地の一要素であり,地方で個人や家族の運営する場合が多く見受けられる.
cultivation.jpg

4.2.7. 裸地

裸地(Bare Ground)は,人の影響ではなく自然に発生した木のない地帯全てに当てはまります.このような地域はアマゾンで自然に発生し,パンタナール湿地やセナード乾燥地帯といった比較的小さい地域もこれに含まれる.

bare.jpg

4.3. 珍しい地理ラベル

4.3.1. 焼畑

移動耕作地の一部に含まれる.濃い茶色や黒い色を示す.

slashburn1.jpg

4.3.2. 伐採

チークやマホガニーといった,高価な木を伐採した地帯を示すラベル.むき出しの茶色いパッチに隣接する曲がりくねった泥道として現れる.

logging1.jpg

4.3.3. 花咲

アマゾンで度々みられる自然現象.

bloom.jpg

4.3.4. 伝統的な採掘

アマゾンには大規模な伝統的な鉱山が存在する.

mine1.jpg

4.3.5. "技工的な"採掘

小さい規模の採掘がこれに当てはまる.アンデスの丘などの金が堆積している地帯に多い.この技工的な採掘には,しばしば違法行為が含まれ,周囲の土壌侵食を引き起こす.

artmine1.jpg

4.3.6. 吹き倒し(Blow Down)

ウィンドスロウとも呼ばれる,アンデス地方に見られる自然現象.アンデス山脈から局所的かつ突発的に乾いた冷たい風がふきおろし(100MPHを超える強力かつ高速の風),大きな熱帯雨林の木々をなぎ倒す.

blowdown.jpg

5. とりあえずデータを分析していく(EDA)

2つのnotebookについて紹介していきます。
5.1.はanokaのEDA
5.2.はPhilipp SchmidtのEDAです

5.1. ラベルのヒトグラム作成

元のnotebookはこちら
ここしばらく衛星写真を用いた画像解析がコンペで多く見受けられます.どうやら最近の流行りのようです.今回のコンペでは,少し珍しい全ての画像にマルチラベルを作成する分類問題です.
そこでここでは,全体的なラベルの様子をざっと分析していきます.

https://gist.github.com/TomHortons/d766738d4ce4bd564a96bbdd5529bfaa

衛星写真の番号とラベルは別ファイルで提供されています。
それらを読み込み、plotlyでヒストグラムを作成します。

スクリーンショット 2017-07-18 16.10.04 1.png

初めのデータ紹介であったように、熱帯雨林のラベル(primary)と気象状況ラベル(cloudy, clear)が最も多いです。
今回の評価指標はF2スコアなので、Recall重視です。
一般的な地理ラベルを抑えつつ、どれだけ珍しい地理ラベルを推定できるか、が勝負になります。

https://gist.github.com/TomHortons/68023b27e469bdb25e0249477ba9f123

次にラベルの組み合わせ頻度をチェックします。例えば、メインとなる熱帯雨林の中に、採掘場、焼畑、吹き倒しなどのラベルがくっついているはずです。このように、「共起(co-occurrence)」はそれ自体が情報となりえます。

スクリーンショット 2017-07-18 16.21.48.png

「農耕地と耕作地、道路」「道路と居住地、伝統的な採掘」など、組み合わせに傾向があることがわかります。

5.2. 画像のクラスタリングの植生分布の可視化

こちらが紹介元のnotebookです。
要素ごとにぶつ切りにしてgistへアップしているので、うまく動作しない場合はこのnotebookからコードをコピペしてください。

5.2.1. ラベルと画像のチェック

from glob import glob
import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
from subprocess import check_output
import matplotlib.pyplot as plt
from scipy.stats import bernoulli
import seaborn as sns
print(check_output(["ls", "../input"]).decode("utf8"))
%matplotlib inline


df = pd.read_csv('../input/train_v2.csv')
image_paths = sorted(glob('../input/train-jpg/*.jpg'))[0:1000]
image_names = list(map(lambda row: row.split("/")[-1][:-4], image_paths))
image_names[0:10]

plt.figure(figsize=(12,8))
for i in range(6):
    plt.subplot(2,3,i+1)
    plt.imshow(plt.imread(image_paths[i]))
    plt.title(str(df[df.image_name == image_names[i]].tags.values))

スクリーンショット 2017-07-18 18.26.11.png

見た所、画像のサイズは全て同じサイズのようです。
各画像にそれぞれ複数のタグがついていることがわかります。

5.2.2. 単純な分類機を作成してみる

https://gist.github.com/TomHortons/d7e9dc9382f6fcd4763163e1a7db9f99

上のコードはscikit learnのLogisticRegressionで分類した人のサンプルコードです。
データの前処理で、以下のあたりが少し煩雑な書き方になっています。
dfは先のサンプルでreadしたイメージデータです。

df['split_tags'] = df['tags'].map(lambda row: row.split(" "))
lb = MultiLabelBinarizer()
y = lb.fit_transform(df['split_tags'])
y = y[:n_samples]
X = np.squeeze(np.array([cv2.resize(plt.imread('../input/train-jpg/{}.jpg'.format(name)), (rescaled_dim, rescaled_dim), cv2.INTER_LINEAR).reshape(1, -1) for name in df.head(n_samples)['image_name'].values]))
X = MinMaxScaler().fit_transform(X)

トレーニングデータのタグをdf['tags'].map(lambda row: row.split(" "))で分解後、MultiLabelBinarizerでベクトル化します。
cv2.resizeで画像をリサイズし、np.squeezeでベクトルサイズを調整、最後にMinMaxScalerでデータを標準化します。

実行した結果は以下の通りです。

Average F2 test score 0.6798866894323994
F2 test scores per tag:
[('primary', 0.96369809349699664),
 ('clear', 0.87778940027894004),
 ('cloudy', 0.60324825986078878),
 ('agriculture', 0.38493549729504789),
 ('road', 0.26332094175960347),
 ('partly_cloudy', 0.22288261515601782),
 ('water', 0.1915041782729805),
 ('habitation', 0.17467248908296948),
 ('cultivation', 0.054811205846528627),
 ('bare_ground', 0.032894736842105261),
 ('haze', 0.0090090090090090089),
 ('slash_burn', 0.0),
 ('conventional_mine', 0.0),
 ('selective_logging', 0.0),
 ('blow_down', 0.0),
 ('blooming', 0.0),
 ('artisinal_mine', 0.0)]

珍しい地理情報がほぼ全滅しています。
F2スコアは0.67となっていますが、トレーニング用データを学習用とテスト用に分割しているので、実際のテストデータを使用するとこれよりはるかにスコアが低下するはずです。

5.2.3. 画像データをクラスタリングする

画像データをそのままクラスタリングしています。
https://gist.github.com/TomHortons/7786a01020de08ee4d14cfaee0ebe142
階層クラスタ、t-SNE、そして外れ値の発見を行います。

手順としては、
1. 画像をnp.arrayへreshape
2. ペアワイズ距離で画像間の距離を計算
3. seabornのclustermapで階層クラスタリング

クラスタした結果がこちらです。

ダウンロード (1).png

もう一つ、t-SNEでクラスタリングもしています。
先ほど作成した画像の距離データをt-SNEで3次元へ落とし込み、plotします。

https://gist.github.com/TomHortons/da0094c11bf77108c41529cf515f78df

実行するとこのようになります。

[t-SNE] Computing pairwise distances...
[t-SNE] Computing 91 nearest neighbors...
[t-SNE] Computed conditional probabilities for sample 600 / 600
[t-SNE] Mean sigma: 2.026685
[t-SNE] Iteration 25: error = 0.9391028, gradient norm = 0.0111124
[t-SNE] Iteration 50: error = 0.8462010, gradient norm = 0.0094869
[t-SNE] Iteration 75: error = 0.6071504, gradient norm = 0.0019813
[t-SNE] Iteration 100: error = 0.5155882, gradient norm = 0.0052735
[t-SNE] KL divergence after 100 iterations with early exaggeration: 0.515588
[t-SNE] Iteration 125: error = 0.4028054, gradient norm = 0.0006756
[t-SNE] Iteration 125: gradient norm 0.000676. Finished.
[t-SNE] Error after 125 iterations: 0.515588

スクリーンショット 2017-07-19 13.12.26.png

作成した距離情報から、外れ値を探していきます。
https://gist.github.com/TomHortons/3aec684789d64cc0631ad5834fa67222

maximally_dissimilar_image_idx = np.nanargmax(np.nanmean(sq_dists, axis=1))

先ほど作成したsq_distsから平均距離を算出し、その最大値(最小値)をピックアップするだけです。
下図の左が距離最大、右が距離最小です。
ダウンロード (1).png

ここでいう距離最小という意味は、全画像中でpairwise距離が最小、つまり最も平均的な画像という意味です。
逆に距離が最大となることは、珍しい画像と言えます。
ただし、単純に画像のpairwiseで算出しただけの距離なので、珍しい建造物がある、変わった形の道路がある、といった意味での「珍しさ」とは異なります。

最後に、t-SNEでマッピングした結果を画像でプロットしています。
https://gist.github.com/TomHortons/8afcf86ee53e073ad21b7d3d0eb4a9fa

ダウンロード (1).png

5.2.4. NDVI (Normalized Difference Vegetation Index)

衛星写真やドローンで撮影された航空写真を使い、植物の活性度を数値化することがしばしばあります。
一般的な指標としてNDVIが使用され、RGBのR、GBの比率から算出されます。
スクリーンショット 2017-07-19 14.46.06.png

今回のデータを使ってNDVIを計算しています。

https://gist.github.com/TomHortons/8d3de37dee1dc7962abe1380f0ff0395

NDVI自体の計算は非常にシンプルで、これ一行で計算が完了です。

ndvis = [(img[:,:,3] - img[:,:,0])/((img[:,:,3] + img[:,:,0])) for img in imgs]

あとはこのような感じで、NDVI画像、標準化したtiff画像、そのままのjpg画像と並べれば雰囲気がつかめると思います。

ダウンロード (2).png

NDVIでは緑色の部分だけが活性化されています。

さらに、全ての画像でNDVIを算出し、その平均値をヒストグラムで確認すれば、植物の活性分布がなんとなくつかめます。

import seaborn as sns
mndvis = np.nan_to_num([ndvi.mean() for ndvi in ndvis])
plt.figure(figsize=(12,8))
sns.distplot(mndvis)
plt.title('distribution of mean NDVIs')

ダウンロード (3).png

最後に、NDVIでランキング付けしてみます。

sorted_idcs = np.argsort(mndvis)
print(len(sorted_idcs))
plt.figure(figsize=(12,8))
plt.subplot(221)
plt.imshow(imgs[sorted_idcs[0]])
plt.subplot(222)
plt.imshow(imgs[sorted_idcs[50]])
plt.subplot(223)
plt.imshow(imgs[sorted_idcs[-30]])
plt.subplot(224)
plt.imshow(imgs[sorted_idcs[-11]])

ダウンロード (4).png

左上の白っぽい画像は、平均NDVIがマイナスの曇り画像、
右下の標準化されて植物が赤っぽく見える画像は、平均NDVIが最大の森林の画像です。