Lesson7

機械学習の実践

1. 学習の目標

このレッスンでは、実際に機械学習モデルを作っていきましょう。

2. 分類の機械学習

はじめに、分類の機械学習モデルを作っていきます。

2.1 実践1(航空写真の分類)

このチャプターで作るものと使用するモデル

航空写真にゴルフ場が映っているか否かを判定する分類の機械学習モデルを作ります。分類モデルには、前のレッスンでも用いた SVC を使います。

まずはJupyterLabの新しいノートブックを作成し、名前を Lesson7_1_nagaoka に変更してください。

必要なライブラリのインポート

ここでは、NumpyやPandas、Pillow、scikit-learn、Matplotlibを使います。次のようにしてインポートします。また画像やグラフをノートブック上で表示するため、%matplotlib inlineとしておきます。

# 必要なライブラリのインポート
import numpy as np
import pandas as pd
from PIL import Image
from sklearn.svm import SVC
from sklearn import metrics
import matplotlib.pyplot as plt

%matplotlib inline

学習に使うデータセットをインポートする

学習に使うデータセットをインポートします。あらかじめデータセットを用意しているので、下記をダウンロードしてください。

航空写真データをダウンロードする

今回の航空写真のデータは、新潟県長岡市の航空写真のオープンデータを使います。以下のサイトで配布されている、航空写真のJPEGファイルからピックアップしています。
データ配布元:地理情報(GIS)航空写真 | 新潟県長岡市

ダウンロードしたファイルを解凍すると nagaoka_gc というフォルダができます。この nagaoka_gc フォルダをそのままCloud9の Lesson7_1_nagaoka.ipynb ファイルと同じ場所にアップロードしてください。

配布データについて

nagaoka_gc フォルダの中には traintest フォルダがあります。

l7_02

以下のような写真が格納されています。

l7_01

また、学習データとして、あらかじめ「ファイル名」と「ゴルフ場が映っているか否かの値(0:映っていない、1:映っている)を示したCSV形式ファイルを用意しました。それぞれ traintest フォルダに格納してあります。

CSV形式ファイルの例

File name GC
08DE663C 0
08EE752D 1
…略… …略…

作った機械学習モデルに画像データを渡すと、そこにゴルフ場が写っているかを判定し、写っていれば「1」、写っていなければ「0」が返されるようになります。

訓練データの読み込み

配布データをアップロードしたら、まずは訓練データを読み込みます。nagaoka_gc/train/train_data.csv です。Pandasを使って読み込んで、先頭5行を表示してみましょう。File nameの項目が画像ファイル名で、実際に、同フォルダに、その名前のJPEG画像ファイルがあります。GCは、ゴルフ場が写っているときは「1」、写っていないときは「0」のデータです。

# 訓練データ用CSVの読み込み
train_data = pd.read_csv("nagaoka_gc/train/train_data.csv")
train_data.head()

出力結果:

l7_03

テストデータの読み込み

同様にしてテストデータも読み込みます。nagaoka_gc/test/test_data.csv です。

# テストデータ用CSVの読み込み
test_data = pd.read_csv("nagaoka_gc/test/test_data.csv")
test_data.head()

l7_04

画像データの読み込みと前処理

この実践では、画像をデータとして扱うため、scikit-learnで使えるように、読み込んだり前処理(画像の加工)したりすることを検討します。

Pillowのレッスンで説明しましたとおり、一般的なカラー画像を読み込むと 画像の高さ(縦), 画像の幅(横), 3(RGB) という3次元配列のデータとなっています。今回配布している航空写真データは、高さ600px・幅800pxのデータです。つまり、1つの画像につき 600 x 800 x 3 = 1440000 個も配列要素が存在することになり、Cloud9に大きな負荷がかかります。

そこで「色の情報は無視して 白黒の明度だけ でも分類の判断に役立つのではないか」と考え、画像のグレースケール化を検討します。

データの確認

まずは、ゴルフ場が写っていない画像を、ひとつ、グレースケールとして表示してみましょう。

# ゴルフ場が写っていない写真をグレースケール化して表示
sample_img1 = Image.open("nagaoka_gc/train/08DE663C.jpg").convert("L")
plt.imshow(sample_img1, cmap="gray")

出力結果:

l7_05

同様に、ゴルフ場を写っている画像を表示すると、次のようになります。

# ゴルフ場が写っている写真をグレースケール化して表示
sample_img2 = Image.open("nagaoka_gc/train/08EE752D.jpg").convert("L")
plt.imshow(sample_img2, cmap="gray")

出力結果:

l7_06

ゴルフ場のある場所は、白黒の明度だけでも特徴的であることがわかります。

画像データサイズの確認

念のため、グレースケール化した画像のサイズを確認します。

# グレースケール化した写真をndarrayに変換してサイズを確認
sample_img1_array = np.array(sample_img1)
sample_img1_array.shape

出力結果:

(600, 800)

出力結果から、600×800の要素数で構成されていることがわかります。要素数は 600 x 800 = 480000 となり、カラー画像を扱うよりもサイズを削減できるので、Cloud9にかける負荷も軽減できます。

そこで今回は、グレースケール化した画像を学習に使います。

データを訓練データとテストデータに分ける

データの分割について

今回配布した写真は、あらかじめ訓練データ(train)とテストデータ(test)に分けています。そのため、データの分割は必要ありません。説明変数(X)に画像データ、目的変数(y)にCSVファイルのGC列のデータを読み込めばOKです。

訓練用画像データの読み込み

以上を踏まえて、訓練用の画像データを読み込みます。

下記の例は X_train に訓練用の画像データ(JPEG形式ファイルから読み込んだもの)、y_trainにゴルフ場が写っているかを分類したデータ(0か1)を格納するプログラムです。

なお、配布した写真は、数が非常に少ないものとなっています。これは、配布元で配布されている写真に、ゴルフ場の写っているものが少なかったためです。画像データが少なすぎると学習精度が悪くなるため、今回は 水増し の処理を行います。

下記では、学習データを水増しする目的で、元データのほか「左右反転したもの」「上下反転したもの」「180度回転させたもの」と、ひとつの画像データから4パターンのデータを作成しています。それぞれを X_train というndarrayに保存しています。

# 訓練用の航空写真の読み込み

# ndarrayのデータを保管する領域の確保
train_len = len(train_data)
# 左右、上下、180度回転させたものを用意するため、4倍の容量を確保する
X_train = np.empty((train_len * 4, 480000), dtype=np.uint8)
y_train = np.empty(train_len * 4, dtype=np.uint8)

# 画像ひとつひとつについて繰り返し処理
for i in range(len(train_data)):

    # 基の画像をndarrayとして読み込んで訓練データに追加
    name = train_data.loc[i, "File name"]
    train_img = Image.open(f"nagaoka_gc/train/{name}.jpg").convert("L")
    train_img = np.array(train_img)
    train_img_f = train_img.flatten()
    X_train[i] = train_img_f
    y_train[i] = train_data.loc[i, "GC"]

    # 左右反転させたものを訓練データに追加
    train_img = np.fliplr(train_img)
    train_img_f = train_img.flatten()
    X_train[i + train_len] = train_img_f
    y_train[i + train_len] = train_data.loc[i, "GC"]

    # 上下反転させたものを訓練データに追加
    train_img = np.flipud(train_img)
    train_img_f = train_img.flatten()
    X_train[i + train_len * 2] = train_img_f
    y_train[i + train_len * 2] = train_data.loc[i, "GC"]

    # 180度回転させたものを訓練データに追加
    train_img = np.rot90(train_img, 2)
    train_img_f = train_img.flatten()
    X_train[i + train_len * 3] = train_img_f
    y_train[i + train_len * 3] = train_data.loc[i, "GC"]
補足:Cloud9の環境は通常、メモリの容量が少ないため、dtype=np.uint8というキーワード引数を設定し、ムダにメモリを専有しないようにしています。
おさらい:文字列の中に変数の値を埋め込むときは、{ }を使って文字列の中に変数を埋め込みます。その際、クォーテーションのすぐ前にfの文字をつける必要があります。
テスト用画像データの読み込み

同様にして、テスト用画像データを読み込みます。ここでは、X_testy_testという変数に読み込みます。学習させるためのデータではないので、水増しは不要です。

# テスト用の航空写真の読み込み

# ndarrayのデータを保管する領域の確保
test_len = len(test_data)
X_test = np.empty((test_len, 480000), dtype=np.uint8)
y_test = np.empty(test_len, dtype=np.uint8)

# 画像ひとつひとつについて繰り返し処理
for i in range(test_len):

    # ndarrayとして読み込んで訓練データに追加
    name = test_data.loc[i, "File name"]
    test_img = Image.open(f"nagaoka_gc/test/{name}.jpg").convert("L")
    test_img = np.array(test_img)
    test_img_f = test_img.flatten()
    X_test[i] = test_img_f
    y_test[i] = test_data.loc[i, "GC"]

モデルを作って学習する

モデルを作って学習します。ここでは、線形分類モデルを作って学習します。画像データは少し大きいので、実際に実行すると、完了までに少し時間がかかります。なお、random_state=1 は、毎回同じ学習結果を出力するために追記しているキーワード引数です。

# 分類器の作成
classifier = SVC(kernel="linear", gamma="scale", random_state=1)
classifier.fit(X_train, y_train)

出力結果:

SVC(C=1.0, cache_size=200, class_weight=None, coef0=0.0,
    decision_function_shape='ovr', degree=3, gamma='scale', kernel='linear',
    max_iter=-1, probability=False, random_state=1, shrinking=True, tol=0.001,
    verbose=False)

期待する性能が出たかを評価する

できたモデルを使って、性能を評価しましょう。テストデータを、このモデルにテスト用の画像データを入れて、分類した結果を表示します。テスト用の画像データは、すでに X_test へ読み込んでいるので、それを指定します。

# 分類の実施と結果表示
y_pred = classifier.predict(X_test)
y_pred

出力結果:

array([0, 0, 0, 1, 0, 1], dtype=uint8)

対応する正解と比較してみましょう。正解はすでに y_test へ読み込んでいますから、その内容を確認します。

# 正解の表示
y_test

出力結果:

array([0, 1, 0, 1, 0, 1], dtype=uint8)

混同行列を使って正解数を確認しましょう。

# 混同行列で正答数の確認
print(metrics.confusion_matrix(y_test, y_pred))

出力結果:

[[3 0]
 [1 2]]

今回の混合行列は、次の値を示しています。ここで、陽性とは「1」のこと、陰性とは「0」のことです。この例であればゴルフ場があれば「1」、なければ「0」です。

  • 真陽性(TP:True Positive):テストデータが陽性で予測も陽性(正解)
  • 真陰性(TN:True Negative):テストデータが陰性で予測も陰性(正解)
  • 偽陽性(FP:False Positive):テストデータが陰性で予測は陽性(間違い)
  • 偽陰性(FN:False Negative):テストデータが陽性で予測が陰性(間違い)

陽性と陰性が以下のように並んでいます。

[[真陰性, 偽陽性],
 [偽陰性, 真陽性]]

今回のゴルフ場の例では、

[[ゴルフ場なしと予測し実際にゴルフ場がなかった数, ゴルフ場ありと予測するも実際にゴルフ場がなかった数(間違い)],
 [ゴルフ場なしと予測するも実際にゴルフ場があった数(間違い), ゴルフ場ありと予測し実際にゴルフ場があった数]
]

です。簡単にいうと、左上と右下が正解した数になります。

さて、「間違い」の部分に着目すると、この学習モデルでは「ゴルフ場ありと予測するも実際にゴルフ場がなかった」という間違いはありませんでした。しかし「ゴルフ場なしと予測するも実際にゴルフ場があった」という間違いが1つあった、ということがわかります。

分類精度を向上させるには、ゴルフ場の写っている写真データを増やす(基になるデータそのものを増やす、もしくは水増しの方法を増やす)ことが有効です。ただしCloud9の無料枠の環境では、画像を増やすとメモリ不足に陥って機械学習の処理がストップする危険性が高まります。
また SVC() のキーワード引数と値を調整するのも手です。具体的な解説まで踏み込むと、少々難しい話になりますので、ここでは省略します。
こういった精度向上の方針については本レッスンの「5.2 実践で行った内容の応用」でも触れていますので、人口予測の学習まで完了した後に、あらためて確認してください。

3. 予測の機械学習

次の実践として、数値予測をする機械学習モデルを作っていきます。

3.1 このチャプターで作るものと使用するモデル

過去の人口から、未来の人口を予測する機械学習モデルを作ります。予測には、回帰を使います。

このセクションでは、RESAS(地域経済分析システム)で提供されている、年度ごと都道府県別の人口データを用います。1960年から2009年までの人口データを訓練データとして学習し、2010年から2015年までを予測してみます。そしてこの期間の予測と実際の人口とが、どのぐらい合致しているのかをグラフで表示して確認します。

データ配布元:RESAS 地域経済分析システム

ここでは2つのサンプルを作ります。

  1. 東京都の人口を予測する
    まずは東京都の人口データのみを使って、東京都の人口を予測します。単回帰分析 を使い、「ある年の人口」から「その翌年の人口」を予測する線形回帰モデルを作ります。

  2. 関東地方の人口を予測する
    次に応用として 重回帰分析 を使い、関東地方の各県について、同じく「ある年の人口」から「その翌年の人口」を予測します。使用するモデルは線形回帰モデルだけでなく、「ランダムフォレスト回帰」と「サポートベクトル回帰」のモデルも試してみます。

3.2 実践2-1(単回帰分析)

では、最初に東京都のみの人口予測を行いましょう。

JupyterLabの新しいノートブックを作成し、名前を Lesson7_2_1_tokyo に変更してください。

使用するモデルについて

上述のとおり、線形回帰モデルを利用して単回帰分析を行います。

必要なライブラリのインポート

以下では、NumpyやPandas、scikit-learn、Matplotlibを使います。次のようにしてインポートします。また画像やグラフをJupyter Notebookのインラインで表示するため、%matplotlib inlineとしておきます。

# 必要なライブラリのインポート
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
%matplotlib inline

学習に使うデータセットをインポートする

学習に使うデータセットをインポートします。あらかじめデータセットを用意しているので、下記をダウンロードしてください。

人口統計データをダウンロードする

今回のデータは、RESASのサイトで「人口マップ」→「人口構成」を表示した際に取得できるデータを加工せずに配布しています。

ダウンロードしたファイルを解凍すると population_data というフォルダができます。この population_data フォルダをそのままCloud9の Lesson7_2_1_tokyo.ipynb ファイルと同じ場所にアップロードしてください。

配布データについて

population_data フォルダには japan_population.csv というCSVファイルのみが格納されています。こちらが今回のデータセットです。

l7_02_01

CSVデータの読み込み

アップロードしたCSVファイルを読み込み、head() で先頭5行を確認してみましょう。

このCSVファイルは、Pythonで処理しやすいよう、UTF-8の文字コードで保存してあります。そのため、Excelなどで開くと文字化けします。
# 人口構成のCSVファイルの読み込み
df = pd.read_csv("population_data/japan_population.csv")
df.head()

出力結果:

l7_02_02

出力結果を見るように、「集計年」「都道府県のコード」「都道府県名」「総人口」で構成されています。全件数も確認しておきましょう。

# 全件数の確認
df.shape

出力結果:

(3055, 4)
東京都のデータだけを取り出す

では、東京都の人口を予想するモデルを作っていきましょう。第1段階として、読み込んだCSVデータから、東京都のものだけを取り出します。

tokyo = df[df["都道府県名"] == "東京都"]
tokyo.head(10)

出力結果:

l7_02_03

人口予測は、ある年の人口から翌年の人口を予想する方式とします。そこで、1960年の人口から1961年の人口を予想する、1961年の人口から1962年の人口を予想する…。というように、説明変数(X)を n年の人口 とし、目的変数(y)を n+1年の人口 とします。

そこで、次のようにして、説明変数 X と目的変数 y を作成します。

# 1960年から2015年までの東京都の総人口数データを機械学習にかける。
# ある年の総人口数を説明変数X、その翌年の総人口数を目的変数yに設定する
# (Xとyの両方とも、縦56行×横1列のndarrayで作成する)
X = np.empty((56, 1), dtype=np.uint32)
y = np.empty((56, 1), dtype=np.uint32)

# 人口はデータの左から3番目
for i in range(56):
    X[i][0] = tokyo.iloc[i, 3]
    y[i][0] = tokyo.iloc[i + 1, 3]

X の値を、簡単に先頭10件ほどで確認しておきます。

X[0:10]

出力結果:

array([[ 9683802],
       [ 9967000],
       [10224000],
       [10470000],
       [10668000],
       [10869244],
       [11018000],
       [11162000],
       [11286000],
       [11367000]], dtype=uint32)

同じようにして y の値も確認します。

y[0:10]

出力結果:

array([[ 9967000],
       [10224000],
       [10470000],
       [10668000],
       [10869244],
       [11018000],
       [11162000],
       [11286000],
       [11367000],
       [11408071]], dtype=uint32)

データの前処理について

今回配布したデータには欠損値も無く、加工も必要ないので、前処理についてはスキップします。

データを訓練データとテストデータを作る

ここまでで生成した Xy を、訓練データとテストデータに分けます。今回は、訓練データを1960年から2009年まで(配列の49番目まで)、2010年から2015年まで(配列の50番目以降)をテストデータとします。

# 1960年から2009年までを訓練データ、
# 2010年以降をテストデータとして分割する
X_train = X[:50]
X_test = X[50:]
y_train = y[:50]
y_test = y[50:]

モデルを作って学習する

では、モデルを作って学習します。ここでは、線形回帰モデルを作成し、訓練データを読み込ませます。

# 線形回帰モデルの作成と学習の実行
model = LinearRegression()
model.fit(X_train, y_train)

出力結果:

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)

期待する性能が出たかを評価する

できたモデルを使って、性能を評価しましょう。テストデータ(X_test)を、このモデルに入れて、結果を計算して表示します。

# テストデータで「翌年の総人口」予測の実施
y_pred = model.predict(X_test)

予測結果(配列の要素)は実数値になっています。次のようにして整数値に変換して、内容を確認しましょう。

# 予測結果は実数値のため、整数値に変換
y_pred = y_pred.astype(np.uint32)
y_pred

出力結果の一例:

array([[13149404],
       [13185920],
       [13219966],
       [13289004],
       [13376011],
       [13485972]], dtype=uint32)

CSVに記載されている正解値は、y_test に読み込んでありますから、それを表示して比較します。

# 正解の表示
y_test

出力結果:

array([[13198000],
       [13234000],
       [13307000],
       [13399000],
       [13515271],
       [13624000]], dtype=uint32)
グラフで確認する

数値での比較だとわかりにくいので、グラフを描いて比較しましょう。まずは、1961年から2010年までの実測値と、2011年からの予測値を連結したデータ y_pred_gr を作ります。

# 正解値とグラフで比較するため
# 実測値と予測値とを連結させた配列y_pred_grを作成
y_pred_gr = np.concatenate([y_train, y_pred])

そしてこのグラフを描画します。赤線が予想値、青線が実測です。若干、予測値のほうが低いですが、人口の伸び率は、実測値と似ていることがわかります。

# 正解値と予測値のグラフ表示
plt.plot(range(56), y_pred_gr, label='Predicted', color='red')
plt.plot(range(56), y, label='Actual', color='blue')
plt.xlabel('Years')
plt.ylabel('Population')
plt.title("Tokyo's population")
plt.grid(True)
plt.legend(loc = "upper left")

出力結果:

l7_02_04

3.3 実践2-2(重回帰分析)

続いて、東京都だけでなく、関東地方の各県の人口を予測するモデルを作ります。

Lesson7_2_1_tokyo.ipynb ファイルと同じ場所にJupyterLabの新しいノートブックを作成し、名前を Lesson7_2_2_kantoに変更してください。

概要と使用するモデルについて

このモデルでは、入力が「東京都の人口」のように1つの説明変数ではなく、「東京都の人口、神奈川県の人口、埼玉県の人口、千葉県の人口…」と、説明変数が複数あります。

必要なライブラリは、東京都の人口を予測する場合と同じです。ここでは、さらに理解を深めるため、線形回帰モデル以外に、以下の回帰モデルを試してみます。

  • ランダムフォレスト回帰(RandomForestRegressor)
  • サポートベクトル回帰(SVR)
必要なライブラリのインポート

今回の重回帰分析に必要なライブラリをインポートします。

# 必要なライブラリのインポート
import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.svm import SVR

学習に使うデータセットをインポートする

データセットは、東京都の予測で使ったものを引き続き利用します。

# 人口構成のCSVファイルの読み込み
df = pd.read_csv("population_data/japan_population.csv")
df.head()

出力結果:

l7_02_02

関東地方のデータだけを取り出す

関東地方のデータだけを取り出します。東京都の単回帰分析では 都道府県名が"東京都" という抽出方法を行いましたが、関東地方の各県のデータ抽出では「都道府県コード」を使います。関東地方各県の都道府県コードは以下のとおりです。

都道府県コード 都道府県名
8 茨城県
9 栃木県
10 群馬県
11 埼玉県
12 千葉県
13 東京都
14 神奈川県

都道府県コードであれば、関東地方各県を取り出す際に 都道府県コードが8以上14以下 という条件指定で抽出できるので便利です。

なお、今回は、関東地方各県のデータの取り出し処理と「データの前処理」を同時に行います。

データの前処理をする

都道府県コードの数値には特別な意味はありませんが、都道府県コードをそのまま説明変数にしてしまうと、数値の大小に意味があるものとしてモデルが認識してしまうかもしれません。そこで、関東地方各県のコード番号を カテゴリ値 とした上で ダミー変数 化して構成することにします。

説明変数(X)として、下記の表のような7列の配列を作成して、処理します。

今回、神奈川県の項目は持たせず「すべてが0のときは神奈川県」というルールにしました。これは、すべてのカテゴリのダミー変数を用意すると、互いの変数に束縛関係が働き、回帰分析の精度低下の可能性があるためです(多重共線性)。なお多重共線性への対処は、必須ではありませんが、ここでは一例として、この方法を採用しました。

茨城県 栃木県 群馬県 埼玉県 千葉県 東京都 人口 (説明)
1 0 0 0 0 0 2047024 1960年の茨城県の人口
0 1 0 0 0 0 1513624 1960年の栃木県の人口
0 0 1 0 0 0 1578476 1960年の群馬県の人口
0 0 0 1 0 0 2430871 1960年の埼玉県の人口
0 0 0 0 1 0 2306010 1960年の千葉県の人口
0 0 0 0 0 1 9683802 1960年の東京都の人口
0 0 0 0 0 0 3443176 1960年の神奈川件の人口
1 0 0 0 0 0 2053000 1961年の茨城県の人口
0 1 0 0 0 0 1513000 1961年の栃木県の人口

元々のCSVデータには左から2列目(配列の1列目)に「都道府県コード」の項目があり、値8(茨城県)~値14(神奈川県)までが、関東1都6県です。また、東京都の単回帰分析と同様、目的変数(y)は n+1年後の人口 とします。

そこで、次のようなコードを書いて、Xy を作成します。

# 1960年から2015年までの関東地方各県の総人口数データを機械学習にかける。
# ある年の都道府県コードと総人口数を説明変数X、その翌年の総人口数を目的変数yに設定する
# (Xは縦56×7行の横2列、yは縦56×7行の横1列でndarrayを作成する)
X = np.zeros((56 * 7,7), dtype=np.uint32)
y = np.zeros(56 * 7, dtype=np.uint32)

cnt = 0
for i in range(56 * 47):
    pref_id = df.iloc[i, 1]
    population = df.iloc[i, 3]
    next_population = df.iloc[i + 47, 3]

    if pref_id >= 8 and pref_id <= 14:
        if pref_id < 14:
            X[cnt][pref_id - 8] = 1

        X[cnt][6] = population
        y[cnt] = next_population
        cnt += 1

# Xを確認
X[0:10]

出力結果:

array([[      1,       0,       0,       0,       0,       0, 2047024],
       [      0,       1,       0,       0,       0,       0, 1513624],
       [      0,       0,       1,       0,       0,       0, 1578476],
       [      0,       0,       0,       1,       0,       0, 2430871],
       [      0,       0,       0,       0,       1,       0, 2306010],
       [      0,       0,       0,       0,       0,       1, 9683802],
       [      0,       0,       0,       0,       0,       0, 3443176],
       [      1,       0,       0,       0,       0,       0, 2053000],
       [      0,       1,       0,       0,       0,       0, 1513000],
       [      0,       0,       1,       0,       0,       0, 1581000]],
      dtype=uint32)

y についても確認しておきましょう。

# yを確認
y[0:10]

出力結果:

array([2053000, 1513000, 1581000, 2497000, 2356000, 9967000, 3606000,
       2056000, 1513000, 1584000], dtype=uint32)

データを訓練データとテストデータを作る

次に、この Xy から訓練データたテストデータを作ります。1960年から2009年までを訓練データ、2010年から2015年までをテストデータとして構成します。1年が7件分あるので、1960年から2009年までのデータの位置は 50年 × 7県 = 350 番目までとなります。

# 1960年から2009年までを訓練データ、
# 2010年以降のデータをテストデータとして分割する
X_train = X[:350]
X_test = X[350:]
y_train = y[:350]
y_test = y[350:]

モデルを作って学習する

データができたので、モデルを作って学習し、予測します。まずは、LinearRegression を使った線形回帰モデルを試してみましょう。

# LinearRegressionで回帰モデルを作成
model1 = LinearRegression()
model1.fit(X_train, y_train)

出力結果:

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)

期待した性能が出たかを評価する

テストデータ(X_test)を、作成したモデルに入れて、予測を実行してみます。

# LinearRegressionで予測実行
y_pred1 = model1.predict(X_test)

結果を確認すると、次のようになります。

# LinearRegressionの予測結果の表示
y_pred1 = y_pred1.astype(np.uint32)
y_pred1

出力結果の一例:

array([ 2976504,  2011342,  2011854,  7234625,  6247455, 13182582,
        9095368,  2967043,  2003902,  2005009,  7248612,  6248144,
       13219973,  9106668,  2954454,  1996155,  1998231,  7255390,
        6231681, 13254834,  9116352,  2944770,  1990345,  1990484,
        7267011,  6232650, 13325526,  9129909,  2935086,  1984534,
        1983705,  7285410,  6240397, 13414617,  9148308,  2925379,
        1978971,  1978006,  7304327,  6253631, 13527211,  9170788],
      dtype=uint32)

正解と比較してみましょう。

# 正解の表示
y_test

出力結果:

array([ 2960000,  2000000,  2001000,  7209000,  6217000, 13198000,
        9060000,  2947000,  1992000,  1994000,  7216000,  6200000,
       13234000,  9070000,  2937000,  1986000,  1986000,  7228000,
        6201000, 13307000,  9084000,  2927000,  1980000,  1979000,
        7247000,  6209000, 13399000,  9103000,  2916976,  1974255,
        1973115,  7266534,  6222666, 13515271,  9126214,  2905000,
        1966000,  1967000,  7289000,  6236000, 13624000,  9145000],
      dtype=uint32)

割と近い予測結果となりました。

ランダムフォレストとSVRの回帰モデルを試す

同様にして、ランダムフォレスト回帰とサポートベクトル回帰のモデルを、それぞれ試してみましょう。

ランダムフォレスト回帰

まずはランダムフォレスト回帰(RandomForestRegressor)を試します。n_estimators=100 は、警告の表示を回避するために追記しているという認識で構いません。

# RandomForestRegressorで回帰モデルを作成
model2 = RandomForestRegressor(n_estimators=100)
model2.fit(X_train, y_train)

出力結果:

RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,
                      max_features='auto', max_leaf_nodes=None,
                      min_impurity_decrease=0.0, min_impurity_split=None,
                      min_samples_leaf=1, min_samples_split=2,
                      min_weight_fraction_leaf=0.0, n_estimators=100,
                      n_jobs=None, oob_score=False, random_state=None,
                      verbose=0, warm_start=False)

テストデータで予測を実行すると、次のようになります。

# RandomForestRegressorで予測実行
y_pred2 = model2.predict(X_test)

# RandomForestRegressorの予測結果の表示
y_pred2 = y_pred2.astype(np.uint32)
y_pred2

出力結果の一例:

array([ 2970892,  2011011,  2012057,  7260296,  6286004, 13082563,
        9026135,  2968499,  2004395,  2004801,  7260296,  6286004,
       13082563,  9026135,  2961032,  1996628,  2002517,  7260296,
        6261015, 13082563,  9026135,  2952997,  1991491,  1991445,
        7264396,  6261015, 13082563,  9026135,  2947377,  1986717,
        1985324,  7310847,  6278435, 13082563,  9026135,  2944633,
        1983093,  1982940,  7310847,  6320482, 13082563,  9026135],
      dtype=uint32)
サポートベクトル回帰

同様に、サポートベクトル回帰(SVR)のモデルを試します。

# SVRで回帰モデルを作成
model3 = SVR(gamma='scale')
model3.fit(X_train, y_train)

出力結果:

SVR(C=1.0, cache_size=200, coef0=0.0, degree=3, epsilon=0.1, gamma='scale',
    kernel='rbf', max_iter=-1, shrinking=True, tol=0.001, verbose=False)

テストデータを用いた予測結果は、以下の通りです。

# SVRで予測実行
y_pred3 = model3.predict(X_test)

# SVRの予測結果の表示
y_pred3 = y_pred3.astype(np.uint32)
y_pred3

出力結果の一例:

array([3928213, 3928193, 3928193, 3928350, 3928321, 3928364, 3928382,
       3928213, 3928192, 3928192, 3928350, 3928321, 3928364, 3928383,
       3928213, 3928192, 3928192, 3928350, 3928320, 3928363, 3928383,
       3928212, 3928192, 3928192, 3928350, 3928320, 3928362, 3928383,
       3928212, 3928192, 3928192, 3928351, 3928321, 3928361, 3928383,
       3928212, 3928192, 3928192, 3928351, 3928321, 3928359, 3928383],
      dtype=uint32)

上述の y_test と比較すると、残念ながら LinearRegression のほうが良い予測結果を得られた結果となりました。本来、ランダムフォレストやサポートベクトルは非常に強力なモデルですが、データ量が少ない場合は予測精度が下がる傾向にあります。

このように「予測や分類をしたい内容やデータの種類・数によって、もっとも良い予測精度を得られるモデルというのは変わってくる」という点は、ぜひ覚えておいてください。

なお、ここでは説明を省略しますが、性能の評価にはグラフ化したり評価関数を使ったりすると、より明確に違いがわかります。

4. まとめ

このレッスンでは機械学習の「分類」「回帰(予測)」の実践例を実習していただきました。プログラム構築の流れはレッスン6で説明したとおりとなり、どのような題材でも基本的な部分は変わりません。さまざまな機械学習に取り組んでみてください。深層学習(ディープラーニング)を学ぶのも良いでしょう。

5. さらに学習したい方へ

学習内容は以上となりますが、さらにPythonや機械学習を学びたい方向けに、各種サイトや今後の学習内容を紹介します。

5.1 公式ドキュメント

まずはPython3.7の公式ドキュメントです。

本コースで学習したPython文法は最低限です。「オブジェクト指向プログラミング」や「リストの内包表記」などを学んで習得できると、今後のPythonプログラミングがより便利に、なおかつ色々なことができるようになります。ぜひ公式サイトなどで、さらなるPythonの文法知識を身に着けてください。

また、各種ライブラリのドキュメントサイトも紹介します。日本語化されていないのもほとんどですので、英語の苦手な方は苦労するかもしれませんが、ぜひ参考にしてみてください。

5.2 実践で行った内容の応用

このレッスンでは「航空写真データにゴルフ場が写っているかの分類」、「統計データを用いて東京都や関東地方各県の人口変動の予測」といった機械学習を行いました。これらの事例を応用してみましょう。

精度を向上させる

今回の事例では、モデルを作成する際にキーワード引数をほとんど利用せず、デフォルト値で行いました。前のセクションで掲載した公式ドキュメントを参照すると、実にたくさんのキーワード引数の存在が見てとれます。これらのキーワード引数と値を設定して、モデルのチューニングを行うと、精度向上の可能性があります。

たとえば sklearn.linear_modelLinearRegression() は以下のようなキーワード引数を持っています。

キーワード引数 設定できる値 概要
fit_intercept True か False Falseを指定した場合、y = aX + bb が0になる前提で学習を行う。
(引数を指定しない場合は自動でTrue)
normalize True か False Trueを指定した場合、モデルが自動で説明変数をスケーリングする。
(引数を指定しない場合は自動でFalse)
copy_X True か False Trueを指定した場合、説明変数のデータをメモリ内で複製し、
その複製データを使って学習を行う。
(引数を指定しない場合は自動でTrue)
n_jobs 整数値 コンピュータが持つCPUのうち、何個を学習のために使用するかを指定。
-1を指定した場合は全CPU。
(引数を指定しない場合は自動で1)

仮に fit_intercept をFalse、n_jobs は2にしたい場合、以下のように記述します。

from sklearn.linear_model import LinearRegression

# 中略

lr = LinearRegression(fit_intercept=False, n_jobs=2)

n_jobs の値を変更するだけでも計算速度や精度の向上につながります。しかし残念ながらCloud9の無料枠ではCPUが1個のため、Cloud9上のJupyterLabで1以外を指定しても効果が現れませんので注意してください。

他のモデルでもキーワード引数はいくつか設定されています。機械学習の理論的なところを理解してからでないと扱うのが難しいものもありますが、いろいろなキーワード引数と値を設定して試行錯誤してみるだけでも予測や分類の結果が変わりますので、試してみるとよいでしょう。

また、航空写真データのほうでは、さらに画像を水増しすると精度が上がるかもしれません。しかし、Cloud9の無料枠はメモリやCPUのスペックが乏しいため、水増しすればするほど、メモリ不足などCloud9の内部的なエラーによってJupyterLab上の処理が強制的にストップしやすくなりますので注意してください。

他の予測や分類を行う

たとえば、航空写真の分類でしたら「住宅が写っているかの分類」が行えます。統計データでしたら他の県や地方についての人口変動の予測が行えます。もちろん、今回配布しているデータだけでなく、配布元からさまざまなデータを取得して、何らかの機械学習を行うことも可能です。面白そうな分類や予測を検討して、機械学習を試してみてください。

発展的な画像解析を行う

こちらもCloud9上ではスペックの関係で実行するのは難しいですが、「対象物が写真に映っているor映っていない」を発展させて「対象物が写真内のどこに映っているか」を解析することも可能です。ただし、それを実現するプログラムを自分で作るのは非常に難しいため、外部ライブラリを使うと良いでしょう。

参考:オブジェクト検出 YOLO(ReNom)

5.3 データ解析のコンペティションに参加

KaggleSIGNATE など、データ解析に関するオンラインコミュニティサイトがあります。こういったサイトでは、データ解析のコンペティションが開催されているので多くのユーザと精度の高さを競うことができ、さらに、他の参加者のソースコードの閲覧もできます。初心者向けに練習用のコンペティションもありますので、臆せず参加してみてください。

5.4 ほかの機械学習アルゴリズムやディープラーニングを学ぶ

このレッスンでは「ランダムフォレスト」や「サポートベクトルマシン」について軽く触れました。こういったモデルの種類(アルゴリズム)を深く学んだり、深層学習(ディープラーニング)について学んだりしてみてください。高度な数学の知識が必要にはなってきますが、データによって適切なアルゴリズムは異なりますので、さまざまなアルゴリズムを学ぶと、より高い精度のモデルを作成できるようになります。コンペティションでも上位に入れるかもしれません。

なお、ディープラーニングを行う場合は、scikit-learnよりも以下のようなライブラリを使うと良いでしょう。ぜひ学んでみてください。

5.5 応用的な画像加工の方法を学ぶ

Pillowのレッスンでは、PillowやNumPyの配列(ndarray)を利用した画像加工の基礎を学びました。ここでは、より応用的な画像加工の方法を紹介します。

なお、ここでは、以下のロケットの写真画像を利用します。この写真を右クリック→「名前をつけて保存」(ファイル名は H2A.jpg のままで保存)してください。ここでは縮小して表示していますが、元々の写真画像のサイズは 800 x 1067 ピクセルとなっています。

ダウンロードした画像ファイル(H2A.jpg)をCloud9上にアップロードしてください。また、画像ファイルと同じフォルダに新規ノートブックを作成し、Images というノートブック名に変更してください。

とりあえず、写真を表示してみます。併せて、ndarrayに変換したデータを用意しておきましょう。

# ライブラリのインポートなど
from PIL import Image
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

# 画像の読み込み
rocket_org = Image.open("H2A.jpg")

# 画像の表示
plt.imshow(rocket_org)

# 画像をndarrayに変換
rocket_org_array = np.array(rocket_org)

出力結果:

画像の色空間

Pillowで読み込んだデータは、デフォルトで RGB の形式になっています。たとえば、縦が750px、横が500pxの画像を読み込んでndarrayに変換すると、そのデータの形式は (750, 500, 3) という3次元配列になります。3番目の 3R(赤)・G(緑)・B(青) の3色を意味しており、縦750 x 横500 に並んだ各ピクセルのデータ(2次元配列)について、RGBそれぞれの色の強さが 0〜255 の256段階のデータとして保存されています。

また、グレースケール(モノクロ)の画像の場合は、白黒1色の情報のみを保持しており、 0(黒) から 255(白) までの256段階で色の強さをデータ化しています。Pillowで読み込んだ画像をグレースケール化する際は、convert("L") を使います。

この convert() は、"L" 以外にも引数として指定できる文字列があり、RGBやグレースケールとは異なる「色の形式(色空間といいます)」に変換できます。

主な色空間への変換

convert()の引数 概要
L グレースケール(1色のみの情報をもつ)
RGB RGB形式(3種類の情報をもつ)
RGBA RGBに透明度(アルファチャンネル)を加えた4種類の情報をもつ
CMYK C(シアン) M(マゼンタ) Y(黄色) K(黒) の4種類の情報をもつ
HSV H(色の色相) S(色の彩度) V(色の明度) の3種類の情報をもつ

HSVという色空間については、以下のページを参照してください。

HSV色空間 - Wikipedia

また、Pillow以外の画像加工ライブラリ(たとえばOpenCV)では BGR という色空間を扱えるものがあります。こちらは、RGBの色空間とほぼ同じですが、情報の並び順が B(青)・G(緑)・R(赤) となっているだけです。RGBをBGRに変換する場合は、ndarray化した上でRとBを入れ替えればOKです。

このコードは入力・実行する必要ありません

# ndarray化した画像データが img_ary に格納されているとして、
# 以下のようにしてRGBをBGRに変換する
temp = img_ary[:, :, 0]
img_ary[:, :, 0] = img_ary[:, :, 2]
img_ary[:, :, 2] = temp

画像の幅と高さを取得

Pillowで読み込んだ画像データの状態で、画像の幅と高さを取得したい場合は、画像オブジェクトの widthheight のプロパティを調べるとわかります。

画像の幅

rocket_org.width

出力結果:

800

画像の高さ

rocket_org.height

出力結果:

1067

また、size のプロパティを参照することで、幅と高さをまとめて取得できます。

rocket_org.size

出力結果:

(800, 1067)

画像の連結

2つ以上の画像を連結させたい場合は、土台となるPillowの Image オブジェクトを new() で新規作成し、その土台の上に paste() で画像を貼り付けます。以下のサンプルは、同じ画像を2つ横に並べた画像を作成しています。

# 土台の作成
#(new()の最初の引数は色空間の指定、2つ目の引数は(幅, 高さ)のタプル)
rocket2 = Image.new("RGB", (rocket_org.width * 2, rocket_org.height))

# 画像の貼り付け
#(paste()の最初の引数は貼り付ける画像、2つ目の引数は左上の座標(横, 縦)のタプル)
rocket2.paste(rocket_org, (0, 0))
rocket2.paste(rocket_org, (rocket_org.width, 0))

# 新しい画像の表示
plt.imshow(rocket2)

出力結果:

画像の切り抜きと分割

画像を切り抜くにはPillowの crop() を使います。crop() の引数には、左上の座標(横, 縦) と 右下の座標(横, 縦) を並べたタプルを指定します。

# ロケットの部分の切り抜き
rocket_crop = rocket_org.crop((400, 200, 500, 650))

# 切り抜いた画像の表示
plt.imshow(rocket_crop)

出力結果:

crop() とfor文を活用することで、いくつかの画像に分割できます。以下のサンプルは、画像を横に4枚で分割し、一番左の分割画像を表示しています。

# 画像を縦4枚に分割
rocket_crop_list = []
for i in range(4):
    rocket_crop_list.append(rocket_org.crop((0, 200 * i, rocket_org.width, 200 * (i + 1))))

# 下端の画像のみ表示
plt.imshow(rocket_crop_list[3])

出力結果:

複数の画像をGIFアニメーションで保存

save() を使う際、指定のキーワード引数を追加することで、複数の画像を「GIF形式のアニメーション画像のファイル」として保存できます。以下のようにして save() のキーワード引数を指定します。

最初のフレームの画像.save("保存するファイル名",
                     save_all=True,
                     append_images=[2枚目以降の画像],
                     loop=ループ回数,
                     duration=コマ送りのミリ秒数)

save_all=Trueは複数の画像をGIFアニメーションで保存するための基本設定で、append_imagesに2つ目以降の画像をリスト形式で指定、loopはループする回数(指定しない場合は1巡したらアニメーション終了)、durationは次のコマを表示するまでのミリ秒数(指定しない場合は一瞬で次の画像に書き換わる)です。

以下のサンプルは crop() で切り抜いた4枚のデータを rocket_anim.gif で保存する処理を行っています。

# 分割した画像をGIFアニメーションとして保存
rocket_crop_list[0].save("rocket_anim.gif",
                       save_all=True,
                       append_images=rocket_crop_list[1:4],
                       loop=0,
                       duration=500)

作成された rocket_anim.gif をダウンロードし、PCで開いて確認してみてください。

出力結果:

RGBの値を操作して画像加工

Pillowのレッスンでは、ndarray化したRGBの値を操作して、画像を暗くする加工方法を紹介しました。他にも、RGBの値を操作して、さまざまな加工が行えます。

画像の単色化(減色)

「赤の情報だけ残したい」といった単色化を行いたい場合は、以下の手順で処理を行います。

  1. すべて初期値0のndarray(元画像と同じ大きさ)を作成
  2. 元画像から特定のカラーチャンネルのデータをスライスして、1の配列に格納
  3. imshow()で表示

以下のサンプルでは、赤・緑・青それぞれの色のみを残す単色化の処理を行っています。

# 赤のみの単色化
rocket_r = np.zeros((rocket_org.height, rocket_org.width, 3), dtype=np.uint8)
rocket_r[:, :, 0] = rocket_org_array[:, :, 0]

# 緑のみの単色化
rocket_g = np.zeros((rocket_org.height, rocket_org.width, 3), dtype=np.uint8)
rocket_g[:, :, 1] = rocket_org_array[:, :, 1]

# 青のみの単色化
rocket_b = np.zeros((rocket_org.height, rocket_org.width, 3), dtype=np.uint8)
rocket_b[:, :, 2] = rocket_org_array[:, :, 2]

# 単色化した画像を横に並べて表示
rocket_monos = Image.new("RGB", (rocket_org.width * 3, rocket_org.height))
rocket_monos.paste(Image.fromarray(rocket_r), (0, 0))
rocket_monos.paste(Image.fromarray(rocket_g), (rocket_org.width, 0))
rocket_monos.paste(Image.fromarray(rocket_b), (rocket_org.width * 2, 0))
plt.imshow(rocket_monos)

出力結果:

こちらを応用すれば、ひとつの色を減色する処理も可能です。以下のサンプルは、青だけ減色した画像に加工しています。

# 青だけ抜いた(赤と緑のみの)画像に加工
rocket_rg = np.zeros((rocket_org.height, rocket_org.width, 3), dtype=np.uint8)
rocket_rg[:, :, 0] = rocket_org_array[:, :, 0]
rocket_rg[:, :, 1] = rocket_org_array[:, :, 1]
plt.imshow(rocket_rg)

出力結果:

二値化処理

基本的にRGBやグレースケールは、256段階で色の強さを表します。「各ピクセルが一定以上の色の強さかどうか」というというのを調べたい(その結果を画像として表示したい)場合は、以下のような配列の計算を行います。以下のサンプルは、各ピクセルの色の強さが128以上かどうかを調べ、128以上なら「白」・128未満なら「黒」で表示しています(ちなみに、このような処理は二値化処理と呼ばれています)。

# 色の強さ 128 を基準とした二値化
rocket_gr_array = np.array(rocket_org.convert("L"))
rocket_wb1 = rocket_gr_array // 128
plt.imshow(rocket_wb1, cmap="gray")

出力結果:

各ピクセルの値を128でわり算(小数点以下は切り捨て)を行い、その商を配列の各要素に代入しなおしています。グレースケール画像として表示する場合、「0と1しか含まれていないデータ」を白と黒のみ(灰色なし)の画像として表示できます。

もちろん、128以外の数値で二値化処理が可能です。

# 色の強さ 160 を基準とした二値化
rocket_wb2 = rocket_gr_array // 160

# 色の強さ 192 を基準とした二値化
rocket_wb3 = rocket_gr_array // 192

# 二値化した画像を横に並べて表示
rocket_wbs = Image.new("L", (rocket_org.width * 3, rocket_org.height))
rocket_wbs.paste(Image.fromarray(rocket_wb1), (0, 0))
rocket_wbs.paste(Image.fromarray(rocket_wb2), (rocket_org.width, 0))
rocket_wbs.paste(Image.fromarray(rocket_wb3), (rocket_org.width * 2, 0))
plt.imshow(rocket_wbs, cmap="gray")

出力結果:

特定の条件を満たすピクセルを抽出

さきほど紹介した二値化処理を応用することで「特定の色の強さ以上のピクセルデータを抽出」できます。以下の手順で処理を行います。

  1. 特定の色のみのデータを用意する
  2. 二値化処理を実施し、各ピクセルの色の強さがある値以上かどうか調べる
  3. 1の配列と2の配列をかけ算する
  4. 3の配列を imshow() で表示する

3の部分ですが、二値化した配列は0か1しか含まれていないのは上述のとおりです。この配列を画像データとかけ算することで、二値化した配列のピクセルが 0 の部分に対応する画像データのピクセルは、かけ算により 0 となります。また、二値化した配列のピクセルが 1 の場合は、元の画像データの数値がそのまま残ります。このような計算により、特定の色の強さ以上のピクセルデータを抽出することが可能です。

以下のサンプルを参考にしてください。

# 赤のみの画像データの元(rocket_r)を変更しないよう、複製物の作成
rocket_r_trimmed = np.copy(rocket_r)

# 赤の色について、色の強さ 240 を基準とした二値化
rocket_r_wb = rocket_r_trimmed[:, :, 0] // 240

# 赤のみの画像データと二値化したデータのかけ算
rocket_r_trimmed[:, :, 0] *= rocket_r_wb

# 赤の強さ 240 以上のみのピクセルを残した画像データの表示
plt.imshow(rocket_r_trimmed)

出力結果:

特定の条件を満たすピクセルをマスキング

上記の例は、条件を満たすピクセルをそのまま表示しましたが、「別の色をつけて表示したい」という場合はマスキングの処理が必要になります。この場合、Pillowの Image.composite() を使います。

マスキングのために必要な画像は、

  • 土台
  • マスキングで表示したいデータ
  • どこのピクセルを表示するかを示したデータ(マスクデータ)

以上の3つです。以下のサンプルでは、土台は単純に Image.new() で作成しています(すべてのピクセルが真っ黒になります)。また、表示したいデータは元の写真データ(rocket_org)、マスクデータは「赤の色の強さが192以上か、それ未満かで二値化処理したデータ(rocket_r_trimmed)」を指定しています。この3つの画像を使うことで、「元の写真を表示する際、赤の色の強さが192以上のピクセルは切り抜いて表示する(切り抜かれた部分は土台が表示されるので真っ黒になる)。そうでない部分は元の画像のまま表示する」という処理をしているサンプルです。

# 土台
rocket_base = Image.new("RGB", (rocket_org.width, rocket_org.height))

# マスキング用のデータ
rocket_r_trimmed_img = Image.fromarray(rocket_r_trimmed[:, :, 0]).convert("L")

# マスク処理を実行して結果を表示
rocket_masked = Image.composite(rocket_base, rocket_org, rocket_r_trimmed_img)
plt.imshow(rocket_masked)

出力結果:

画像の重ね合わせ

最後に、微妙に異なる2つの画像を重ね合わせる方法を紹介します。まずは、画像ファイルを格納している、以下の圧縮ファイルをダウンロードしてください。

画像ファイルをダウンロード

解凍して現れる home1.pnghome2.png を、Cloud9上の Images.ipynb と同じフォルダにアップロードしてください。

まずは2つの画像をそれぞれ表示してみます。

home1.png

im1 = Image.open("home1.png")
plt.imshow(im1)

出力結果:

home2.png

im2 = Image.open("home2.png")
plt.imshow(im2)

出力結果:

2つの画像の違いは「青い風船(片方の写真には映っていない)」のみです。この2つの画像を「重ね合わせて」みます。

配列のひき算

まずは、単純に2つの画像をひき算してできたデータを表示します。

# 画像をndarray化
im1a = np.array(im1)
im2a = np.array(im2)

# 2つの配列をひき算した結果を表示
im2a -= im1a
plt.imshow(im2a)

出力結果:

撮影時のタイミングにおける光の誤差の影響でノイズが表示されていたり、風船の色も一部が緑色になったりしていますが、風船のみがはっきりと見え、背景部分はほぼ真っ黒です。このように、単純な引き算により、画像内の差分の部分を抽出できます。

Image.blend()

もうひとつ、Pillowの Image.blend() を利用して、重ね合わせた画像を表示します。引数は、1つ目と2つ目が重ね合わせたい画像、3つ目は透明度を指定します。透明度は0から1.0までの小数値で指定します(数値が小さいほど透明になります)。たとえば3つ目に 0.3 と入力すると、2つ目の画像の透明度が 0.3、1つ目のほうは 1 - 0.3 = 0.7 の透明度で表示されます。

im3 = Image.blend(im1, im2, 0.3)
plt.imshow(im3)

出力結果:

風船が薄く表示されます。

他にも、このセクションで紹介した内容を応用すれば、たとえば差分の部分に他の色をつけて表示する、といったことが可能です。

その他

ここで紹介したは画像加工に関する処理は、ほんの一部です。まだまだ紹介しきれなかった画像加工の方法もあれば、Pillow以外のライブラリであれば可能な加工処理(例:OpenCVを使えば「特定の条件を満たす領域にマーキングする処理」が可能)もあります。もっと深く学習されたい方は、書籍やネットなどで自主的に学習を進めましょう。

5.6 Tellusを活用してデータ解析を行う

基礎や理論を学ぶだけでなく、学んだことを実践に活かすことが重要です。Tellusを活用して、実践的なデータ解析を行いましょう。

オウンドメディア「宙畑」について

オウンドメディア「宙畑」では、Tellusを使った解析事例も紹介しています。詳しくは下記サイトよりご覧ください。

https://sorabatake.jp/tag/tellus/

6. Tellusをはじめよう

受講お疲れ様でした。

このコースを受講した皆さんは、開発環境を用いて衛星データを解析する基礎が身につきました。

この記事では「衛星データを用いて画像処理を試してみたい」という方に向けて、Tellusを利用することでどのようなことを実現できるか、その利用方法を含めて以下の内容をご紹介します。

  • Tellusの解析環境の種類
  • ユーザー登録の方法
  • Tellusの搭載データ
  • Tellus OSでできること
  • 統合開発環境の申請手順
  • 統合開発環境でできること

なお、Tellusは2021年3月までは、どなたでも原則無料でご利用いただけます。

6.1 Tellusの解析環境の種類

Tellusは2つのアクセス方法があります。

image6.png

1つはブラウザ上で表示した地図にさまざまなデータを重ね、操作できる Tellus OS、もう1つはプログラミングでデータを解析する 統合開発環境 となります。

Tellus OSではツールによる操作以外にScript Editor上でJavaScriptを用いて解析を行うことができます。2019年12月現在、この機能はβ版ですが、OS上で簡易な解析を行う場合には便利です。

Tellusの開発環境では、各種データとJupyter NotebookがAPIで連携しています。そのため、Jupyter Lab上からプログラミングによりAPIを利用することで、データを取得し、OS上より柔軟かつ高度な解析ができます。

OS上でのJavaScriptを用いた解析との差としては、データ解析に適した豊富なモジュールやパッケージを備えたPythonやRで使用できるという点です。開発環境では無料の仮想環境と専有環境が提供されるだけでなく、さらに高スペックなGPU開発環境も利用することが可能となっています。衛星データの解析には高い処理能力を要求されることがありますが、Tellusで提供している開発環境では、そのような状況にも対応が可能です。

6.2 ユーザー登録の方法

実際にユーザー登録をしましょう。

(1)Tellusのホームページ へアクセスします。

image19.png

(2)「ユーザー登録」をクリックし、登録ページに移ります。入力フォームに必要事項を記載し、約款および個人情報の取扱いについて「同意します」にチェックを入れ、「登録」ボタンをクリックします。

image7.png

(3)登録を済ませると、ユーザー登録受け付けについてのメールがTellusから自動で配信されます。メールに記載されたURLをクリックし、登録を完了させます。ユーザー登録が完了した後にも、登録を完了したことをお知らせするメールが配信されるので確認してください。なお、メールが配信されていない場合には迷惑メールフォルダに格納されている可能性があります。

(4)これでTellus OSを使う準備ができました。再びトップページへ移動し、ログインを済ませて、「Tellus OS」をクリックします。

image23.png

(5)Tellus OSが起動します。無事に起動できれば、下のような画面が現れます(画像では各ツールの説明も加えてあります)。

image21.png

6.3 Tellus OSの搭載データ

現在、Tellus OS上では約40の異なったデータを扱うことができます。データはその種類によって、

  • 衛星データ
  • 地上データ

の2つに大別されています。各々のデータに関する詳細は データカタログ をご覧ください。

加えて、Tellus OSでは、衛星データに付帯する「メタデータ」という情報に基づいて画像を取得できます。

Tellus OSの画面左上に虫眼鏡のアイコンがあります。その隣にあるのがメタデータ検索のアイコンです。そこをクリックすると、新しいウィンドウが開きます。このウィンドウから任意のメタデータを指定できます。検索では、衛星の種類や期間の指定だけでなく、他にも複数の検索条件を設定でき、また検索結果に対しては観測日や被雲率などで並べ替えも可能です。

image4.png

6.4 Tellus OSでできること

Tellus OSでは、基本的にツールをクリックして解析を進めることになります。たとえば、OS上では以下のようなことができます。

(a) さまざまな衛星画像の観察

image22.png

(b) 特定の場所に植物が繁茂しているかの判別

image12.png

上図(右)ではPresetでFalse Colorを使用したことで、植生部分が緑色で表示されている

(c) 任意の地点からの植物割合の算出

image9.png

上図のように対象をポリゴンで囲む。囲まれた部分の面積から割合を算出できる

(d) 異なる衛星(地上)データの比較

image17.png

上図の左はNDVIを用いた図。右は土地利用図。どちらの図でも緑地が緑色で抽出されているのがわかる

(e) 異なる衛星(地上)データの重ね合わせ

image15.png

上図は人流マップとアメダス(降雨量)を重ねたもの

(f) 外部ファイルを読み込んで衛星(地上)データと比較

image3.png

上図は鎌倉市が提供しているオープンデータ(公共トイレの場所を示したGeoJSONデータ)をTellus上に表示したもの

詳しくは 公式オウンドメディアの宙畑の記事 をご覧ください。

ツールを利用する以外では、Script Editor(β版)を用いてJavaScriptでTellus OSを操作することもできます。

image1.png

Tellus OSの画面上部のバーにある「Script Editor」ボタンをクリックするとエディター画面が開きます。

image5.png

Script Editorでは、Exampleタブ内にプロジェクトのサンプルを用意しています。自分で書いたスクリプトを保存することも可能です。

6.5 統合開発環境の申請手順

ここまではTellus OSに関する説明をしました。ここからは統合開発環境について説明します。まずは統合開発環境を使えるようにするための申請の手順をご紹介します。

(1)ログインした状態で ユーザーマイページ へアクセスし、ページ内にある「開発環境の申し込み」をクリックします。申し込みフォームでは、申し込みプランを選択し、必要事項を入力し、「送信する」をクリックします。

image2.png

(2)開発環境の申請後、開発環境を利用できるようになるまでにはしばらく時間がかかります。窓口からの開発環境提供のお知らせを、しばらくお待ちください。開発環境が提供されますと、ユーザーマイページの「開発環境利用状況」から、環境ホスト名やトークン情報を確認できます。

image10.png

(3)提供された開発環境へアクセスするには2通りの方法があります。

  • 開発環境利用状況に記載された「環境ホスト名/IP」の情報をコピーし、直接ブラウザのURL窓に貼り付けて、アクセスします。
  • Tellus OSに入り、カメラのアイコン(スクリーンショット)の右にある「開発環境を開く」のアイコンをクリックします。

image13.png

これで統合開発環境へアクセスできます。

6.6 統合開発環境でできること

JupyterLabを利用した統合開発環境では、さまざまなことができます。下に挙げているのはあくまで一例であり、ユーザーの目的に応じて柔軟な解析が可能になっています。

(a) 光学画像の取得(拡大画像も取得)

image11.png

ASNARO-1の衛星画像。左(拡大前)と右(拡大後)

(b) 光学画像の波長帯の調節

image12.png

AVNIR-2衛星画像。左(True color)、右(False color)

(c) 複数の画像の結合

image18.png

SLATSの衛星画像を連続して結合したもの

(d) 衛星画像のクリッピング(一部切り出し)

image8.png

(e) 衛星画像の時系列比較

image14.png

左(2007年8月)、右(2008年8月)の同地点の水田の状態を比べたもの

(f) 衛星から標高データの取得

image16.png

任意地点に対して標高情報を持ったタイルを抽出。それらの高さをヒストグラムで表示

これらの方法についても 別のチュートリアル記事 で詳細な説明があります。ぜひ、そちらをご一読ください。

これでTellus OSと開発環境に関する説明は終わりです。本コースで身につけたスキルを活かし、解析に取り組んでみてください。

6.7 Tellus Trainer

さらに機械学習やTellusについて学んでみたい、という方に向けて、Tellusの利用方法を学ぶe-Learning「Tellus Trainer」のコンテンツも提供しています。

以下はTellus Trainerの概要です。

コンテンツ内容

  • データサイエンス/AI概論
  • データサイエンス講座中級編
  • 衛星データの基礎知識
  • Tellus上の衛星データ画像解析講座
  • 衛星データ×機械学習概論
  • 機械学習における物体検出概論
  • 物体検出演習

参加費

無料

運営

さくらインターネット株式会社
株式会社SIGNATE

身につく知識・スキル

  • 光学衛星、SAR衛星など、衛星データの基本原理、解析方法
  • Tellusの利用を想定した衛星データの解析方法
  • Pythonを使った機械学習スキル

Tellus Trainerの受講方法

受講は下記ページからお願いします。

https://tellusxdp.github.io/tellus-trainer/index.html

7. 最後に

以上で、テキストは終了です。お疲れさまでした!

学習に終わりはありません。本気で仕事に活用していくのであれば、本コースだけでは足りないこともあります。新しい技術や知識をどんどんと取り込み、機械学習やAIの世界を拡げていってください。皆さんはすでに最初の一歩、いや二歩三歩ぐらいは進めることができたはずです。技術は日々進化を遂げていますが、基本的な仕組みの部分はこれからも変わりません。

これからの時代を生き抜くうえで、本コースで学んだことが、皆さんの一助となれば幸いです。

Copyright 2020 SAKURA Internet Inc. All rights reserved