sklearn: ドローン画像でPCA主成分分析 (PCA, GaussianMixture)

# 主成分分析[しゅせいぶんぶんせき]（PCA)は何だろう？
簡単に言うと、PCAはデータの次元減少の手法です。データは主な成分(Principle componets )までに減らして、データ内で重要なことを出てきます。例えば、りんご農園さんはたくさんなデータを解析してみます。
このリンゴデータ内で重要な成分は何でしょう？

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/999291/331c766d-e6d6-b92e-8a3f-18d9176b7709.png)


マーケティングや投資の解析では、よく使っている手法ですが、PCAは画像のセグメンテーションにも応用できます。
# データ
この記事は、[OpenAerialMap](https://map.openaerialmap.org/#/136.16455078125,29.477861195816843,5?_k=gdh1w5)からtifのファイルを利用しています。もし、コードに沿うようにしたい場合は、ファイルをダウンロードするか、自分のものを使ってください。

https://map.openaerialmap.org/#/136.16455078125,29.477861195816843,5?_k=gdh1w5

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/999291/bb540718-f799-110a-cdcc-5fa7775a9c61.png)



# コード
今回の記事のコードはcollabで実行できるノートブックを使います
利用しているデータはtifの形式なので、最初にRasterioをインポートします。
```python
!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)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/999291/6b5a7a1d-08ed-d831-284e-a95875ab334d.png)

## Sklearn PCA
PythonでPCAを実行することはとても簡単です。scikitlearnというモジュールはあなたに代わってすべてを行います。

https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html
```python
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を使っています。

```python

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)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/999291/d622c9fa-ce41-2a03-4c3b-298a55474dd6.png)

画像はPCAで処理するために、2Dのarrayに変換して、正規化も行わなければならないです。（正規化を抜ければ、実行できますが、結果が悪くなる）



## Sklearn GaussianMixture

次　PCAの結果を分類します。Sklearnのモジュールには、複数の分類の手法があり、今回 GaussianMixtureという手法を利用します。

https://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html

```python
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」の希望を使います

```python
from sklearn.mixture import GaussianMixture 
c=5 #クラス数
gmm = GaussianMixture(n_components=c, random_state=42).fit_predict(pca_img)
print(gmm.shape)
```
この結果は画像形式ではないので、表示するために、reshapeをします。[
```python
dims=img.shape
print(dims)
print(np.unique(gmm))
gmm_img=gmm.reshape(dims[0],dims[1])
```
プロットで表示しましょう。そして、綺麗に表示するためにカスタムのColorMapをも作ってみます。

```python
#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)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/999291/c236b008-7cfe-d912-45d9-6e3220104c54.png)

GaussianMixtureの手法は画像のPCAをけっこう分類できます。画像内で川範囲はほとんどが検出できたが、植生や畑や道路と同じにクラスに分類してしまいました。しかし、これは良い最初の結果であり、いくつかの調整と画像処理で、より良い結果を得ることができると思います。
しかし、今回の記事はここまです。次回で、この結果はShpのファイルに変換してみます。よろしくお願いします。