主成分分析[しゅせいぶんぶんせき]（PCA)は何だろう？
簡単に言うと、PCAはデータの次元減少の手法です。データは主な成分(Principle componets )までに減らして、データ内で重要なことを出てきます。例えば、りんご農園さんはたくさんなデータを解析してみます。
このリンゴデータ内で重要な成分は何でしょう？
マーケティングや投資の解析では、よく使っている手法ですが、PCAは画像のセグメンテーションにも応用できます。
データ
この記事は、OpenAerialMapからtifのファイルを利用しています。もし、コードに沿うようにしたい場合は、ファイルをダウンロードするか、自分のものを使ってください。
コード
今回の記事のコードはcollabで実行できるノートブックを使います
利用しているデータはtifの形式なので、最初にRasterioをインポートします。
!pip install rasterio
#インポート
import cv2
import rasterio as rio
import matplotlib.pyplot as plt
#ファイルを読み込み
DATA_LOC="D:/DATA/Okunawa/57e9afd92b67227a79b4fd5c.tif"
src=rio.open(DATA_LOC)
#画像準備
img=src.read().transpose(1,2,0)
img=cv2.resize(img,None,fx=0.25,fy=0.25)
#表示
plt.imshow(img)
Sklearn PCA
PythonでPCAを実行することはとても簡単です。scikitlearnというモジュールはあなたに代わってすべてを行います。
class sklearn.decomposition.PCA(n_components=None, *,
copy=True, whiten=False, svd_solver='auto', tol=0.0,
iterated_power='auto', n_oversamples=10,
power_iteration_normalizer='auto', random_state=None)
重要な変動要因は1つだけで、それはコンポーネントの数です(n_components)。それは何個の主成分まで減少することを示しています。このコンポーネントの数を大きくすると 実行時間が長くなるので、今回 n=3を使っています。
from sklearn.decomposition import PCA
n=3
pca=PCA(n)
pca_img=pca.fit_transform(img.reshape(-1,3)/255.)
print(pca_img.shape) #(34375730, 3)
plt.scatter(pca_img[:,0],pca_img[:,1],c=pca_img[:,2],s=0.3)
画像はPCAで処理するために、2Dのarrayに変換して、正規化も行わなければならないです。（正規化を抜ければ、実行できますが、結果が悪くなる）
Sklearn GaussianMixture
次 PCAの結果を分類します。Sklearnのモジュールには、複数の分類の手法があり、今回 GaussianMixtureという手法を利用します。
class sklearn.mixture.GaussianMixture(n_components=1, *,
covariance_type='full', tol=0.001, reg_covar=1e-06, max_iter=100,
n_init=1, init_params='kmeans', weights_init=None, means_init=None,
precisions_init=None, random_state=None, warm_start=False, verbose=0,
verbose_interval=10
GaussianMixtureでは設定できる変数はたくさんありますが、必要なのはコンポーネントの数だけです。今回、コンポーネントの数は、分類したいクラス数ということを表しています。今、利用している画像では、川、農地、森林、道路、黒背景に分類したいので、クラス数は5に設定してみます。そして、pcaの結果は分類するために、「fit_predict」の希望を使います
from sklearn.mixture import GaussianMixture
c=5 #クラス数
gmm = GaussianMixture(n_components=c, random_state=42).fit_predict(pca_img)
print(gmm.shape)
この結果は画像形式ではないので、表示するために、reshapeをします。[
dims=img.shape
print(dims)
print(np.unique(gmm))
gmm_img=gmm.reshape(dims[0],dims[1])
プロットで表示しましょう。そして、綺麗に表示するためにカスタムのColorMapをも作ってみます。
#Cmap作成
import matplotlib
from mpl_toolkits.axes_grid1 import make_axes_locatable
color_rgb= {'green':[0,255,0], 'blue':[0,0,255], 'yellow':[255,255,0], 'grey':[220,200,0], 'aqua':[61,142,149]}
def gen_cmap(labels,name):
my_cmap = matplotlib.colors.ListedColormap(np.array(list(labels.values()))/255., name=name)
return my_cmap
cmap=gen_cmap(color_rgb,"cmap")
#プロット
fig,axs=plt.subplots(ncols=3,figsize=(15,5))
#PCAグラフ
for i,col in enumerate(list(color_rgb.keys())):
print(i)
axs[0].scatter(pca_img[gmm==i][:,0],pca_img[gmm==i][:,1],s=0.3,color=col,label=f"k{i}")
axs[0].set_xlabel('PCA1')
axs[0].set_ylabel('PCA2')
axs[0].legend()
#画像
axs[1].imshow(img)
#gmm画像の結果
axs[2].imshow(gmm_img,cmap=cmap)
GaussianMixtureの手法は画像のPCAをけっこう分類できます。画像内で川範囲はほとんどが検出できたが、植生や畑や道路と同じにクラスに分類してしまいました。しかし、これは良い最初の結果であり、いくつかの調整と画像処理で、より良い結果を得ることができると思います。
しかし、今回の記事はここまです。次回で、この結果はShpのファイルに変換してみます。よろしくお願いします。