2022/10/27

はじめに

randomforestとかsvmとかを使って土地被覆分類をしたい！でも教師データは持ってないしGUIベースのソフトに詳しくない！！pythonだけで出来ねぇかな！？？！そんな人（主に自分）のためにメモ書きを残しました。QGISのSCPは3.22のサポート外だしArcGISは高いし…

使用ツールと全体の流れ

QGIS 3.22 (教師データの作成)
python 3.9.2 (教師データと.tifを使い分類)
モジュールはscikit-learn(サイキットラーン), gdal,pandasを使用
QGIS 3.22 (分類後の画像を見るため)

参考サイト

今回参考にするコード

import numpy as np
import gdal
import pandas as pd
from sklearn.ensemble import RandomForestClassifier

#######define inputs#######
inpRaster = 'C:/Path/to/Raster/raster.tif' 
outRaster = 'C:/Path/to/output/output.tif'
df = pd.read_csv('C:/Path/to/training/training_data.csv', sep=',')

########end of inputs#####

#read training data
#enter training data bands according to your csv columns name
data = df[['S2_1','S2_2','S2_3','S2_4','S2_5','S2_6','S2_7','S2_8','S2_9','S2_10']]
#enter training label according to your csv column name
label = df['COD']
del df

####no need to modify the code below###
#######################################


#open raster
ds = gdal.Open(inpRaster, gdal.GA_ReadOnly)

#get raster info
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount
geo_transform = ds.GetGeoTransform()
projection = ds.GetProjectionRef()

#read as array
array = ds.ReadAsArray()
ds = None

#modify structure
array = np.stack(array,axis=2)
array = np.reshape(array, [rows*cols,bands])
test = pd.DataFrame(array, dtype='int16')
del array


#set classifier parameters and train classifier
clf = RandomForestClassifier(n_estimators=50,n_jobs=-1)
clf.fit(data,label)
del data
del label


#predict classes
y_pred = clf.predict(test)
del test
classification = y_pred.reshape((rows,cols))
del y_pred

def createGeotiff(outRaster, data, geo_transform, projection):
    # Create a GeoTIFF file with the given data
    driver = gdal.GetDriverByName('GTiff')
    rows, cols = data.shape
    rasterDS = driver.Create(outRaster, cols, rows, 1, gdal.GDT_Int32)
    rasterDS.SetGeoTransform(geo_transform)
    rasterDS.SetProjection(projection)
    band = rasterDS.GetRasterBand(1)
    band.WriteArray(data)
    rasterDS = None

#export classified image
createGeotiff(outRaster,classification,geo_transform,projection)```

https://www.youtube.com/watch?v=H9MnS2oFN7I&ab_channel=JesseBuyungo (1個目の解説動画)

とてつもなくこのコードに助けられました。danielm09さんありがとうございます

`trainingdata.csv`を作成

danielm09氏のGitHub issueでは、このcsvファイルはQGIS上の「sample raster values」により作成されたとあった。日本語版QGIS3.22では、この機能は「ベクタレイヤにラスタ値を追加」と訳されている。

マルチポイント形式の.shpを作成し、土地被覆とポイントデータを結びつける

QGISメイン画面左上の「新規shapefileレイヤ」から画像の画面に行き、ファイル名(必ず...から絶対パスを教えてあげること)・ジオメトリタイプ・新フィールドから属性テーブル　の三つを作成しOKを押す
フィールドリストに必要なのはidとその点の土地被覆の種類（土壌、森林、畑など）（以下genre）の二点

自分はw(water or not mangrove), sm (sparce mangrove), dm (dense mangrove), b(black) の4つで分類しました！

.shpを作成後、これをレイヤ上で選択した状態でQGIS上側の編集モード切り替え（鉛筆マーク）→点地物を作成→作成したポイントデータにidとgenreごとに統一した名称をつけていく
true colorが表示されている.tifなどを参考にしながら教師データを作成

.shpにRGBやNIRなどのバンド情報を与える

次に、danielm09氏も使用していた「ベクタレイヤにラスタ値を追加」ツールを使用する。

入力レイヤに先ほど作成したshpファイル、ラスタレイヤには.shpと結びつけたいバンド情報を持つ.tifを与える。サンプリングした出力の欄で …から「ファイルに保存」を選択して絶対パスを指定し、 結びつけた後の.shpの名称を入れ、実行ボタンを押そう。

完成した.shpを.csvに変換

ポイントデータが完成したら、レイヤ上で右クリック→エクスポート→新規ファイルに地物を保存を押す。

保存形式を「CSV」にし、CRSが元の.tifと一致していることを確認してexport！

ここでもファイル名を付ける際…のボタンから絶対パスを与えること！！

.fit させる

作成した教師データと元の.tifを使って学習させます

実際に使用したコードを分割して働きについてメモっていきます

import numpy as np
import gdal
import pandas as pd
from sklearn.ensemble import RandomForestClassifier

インポート。

inpRaster = 'C:/Path/to/Raster/raster.tif' 
outRaster = 'C:/Path/to/output/output.tif'
df = pd.read_csv('C:/Path/to/training/training_data.csv', sep=',')

inpRaster,outRaster,dfのパスを正しく与えましょう。作業ファイルからの間接パスでも行けると思います。C:/Path/to/output/output.tifっていう名前のファイルがなくても自動で作成してくれるけど.../output/っていうフォルダがなかった場合エラーを吐くのでフォルダは作ってあげましょう。（ここで1時間使った顔；〜；）

#read training data
#enter training data bands according to your csv columns name
data = df[["Red","Green","Blue","NIR","SWIR","SWIR2","NDVI","EVI","MNDWI"]]
#enter training label according to your csv column name
label = df['genre']

教師データのcsvファイルをよく見て、トレーニングデータとラベルがどの列（行かも）に参照するかを注意深くメモ書きして代入します。自分は事前にバンド名をわかりやすく加工しているためこのようなプログラムになりましたが、ここは各人によって異なります。

ds = gdal.Open(inpRaster, gdal.GA_ReadOnly)

rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount
geo_transform = ds.GetGeoTransform()
projection = ds.GetProjectionRef()

array = ds.ReadAsArray()

inpRasterを読み込み、行や列のサイズ、バンド数、CRSなどの情報を読み込んでarrayに起こしている。ここらへんは最後の書き出しの時に必要な情報になります。

array = np.stack(array,axis=2)
array = np.reshape(array, [rows*cols,bands])
test = pd.DataFrame(array, dtype='float')
test = test.set_axis(["Red","Green","Blue","NIR","SWIR","SWIR2","NDVI","EVI","MNDWI"],axis="columns")

arrayのデータ形を整えています。前半2行で行っている具体的な内容はyoutubeの解説動画に詳しいです。後半2行で、arrayをDataFrameにし、カラム名をdataのものと一緒にしています。これを変えないと、

UserWarning: X does not have valid feature names, but IsolationForest was fitted with feature name

というエラーを吐かれます。

clf = RandomForestClassifier(n_estimators=50,n_jobs=-1)
clf.fit(data,label)

dataとlabelを用いて学習させています！多分RandomForestClassifierの部分を他の分類器に変えれば色々応用が効きます！sklearn大好き♡

y_pred = clf.predict(test) 

classification = y_pred.reshape((rows,cols))

作成したclfという分類器を用いて、testデータを分類させています。また、1行目で作成されたy_predという変数は一次元のものなので、これを2次元の形に戻し、classificationに代入します。

for i in range(len(classification)):
    for j in range(len(classification[0])):
        if classification[i][j] == "w": # water or not mangrove 
            classification[i][j] = 0
        elif classification[i][j] == "sm": #spurce mangrove
            classification[i][j] = 1
        elif classification[i][j] == "dm": # dense mangrove
            classification[i][j] = 2
        elif classification[i][j] == "b": # black
            classification[i][j] = 3
        else:
            raise ValueError("not well classificated!")

classification内には分類結果が文字データとして格納されていますが、このままではQGIS上で表示させることができません。（当たり前だけど僕はここで3時間使いました）そのため、適当な数字に変換してあげます。

def createGeotiff(outRaster, data, geo_transform, projection):
    # Create a GeoTIFF file with the given data
    driver = gdal.GetDriverByName('GTiff')
    rows, cols = data.shape
    rasterDS = driver.Create(outRaster, cols, rows, 1, gdal.GDT_Int32)
    rasterDS.SetGeoTransform(geo_transform)
    rasterDS.SetProjection(projection)
    band = rasterDS.GetRasterBand(1)
    band.WriteArray(data)
    rasterDS = None

createGeotiff(outRaster,classification,geo_transform,projection)

.tifを作成する関数を定義し、これを実行してoutput完了です！

感想

エラー文をよく読んで検索して出てくるstackoverflowの英文にビビらず立ち向かおう