ArcGISで教師付分類

最近Deep learningとかRとかでいじってばかりで、久しぶりにArcGISやらQ-GISやらで画像そのものをいじっています。
(そもそも、画像頑張っていたころはER-Mapperとすぐ落ちる頃のGrassだったんですよね。ちょっとGRASSいじったら使い勝手が良くなっていてびっくりしました、相変わらず素人向けではないですが)

画像の分類についてはversion 10くらいから画像分類ウインドウがついて便利になりましたが、マニュアルにあまり記載されていない部分も多くありますし、多数のデータをいじるにはGUIで不便なのでToolboxなどの機能で実行したい場合も多くあります。ここではそうした方法を含めて手順を、備忘録がてらにメモいたします。

さらに、いつの間にかセグメンテーション(10.3から)だの、精度評価だの(10.4)、新機能がついていたとは驚きです。もうデスクトップ版の開発しないんじゃないかと思ってました。64bit版タッチのソフトも、従来のデスクトップ製品も、Online&Serverも全方向で開発がんばってるんですね。すみません。

(で、書き途中でBIOSが飛んでしまい、えらいめにあっています。メモしていたテキストもHDDごと取り出さないとダメそうなので、とりあえずここまで。)

2017/11/16追記:BIOSの件は、結局マザーボードが死んでいました・・・

レーニングデータの準備

レーニングデータは最終的には分類したいカテゴリーごとの値が属性に付与された、マルチパートポリゴンとして表現されている必要があります。
マニュアルでは画像分類のツールバーからトレーニングデータを用意するようにあります。
それでも良いのですが、ここでは、他のソフトを使ったり、人とデータのやりとりなどに便利ように、普通のShapeファイルを作成します。

  • Shapeファイルを作成し、属性に数字でカテゴリー値を入れておく。この際、隣接するポリゴンとの間を埋めるために(埋めなくても良いが、メモように)微細な隙間やズレの修正には[インテグレート] 、個別には[トレース形状に一致]、Infoの[トポロジ] を使用
  • 完成したら、ディゾルブでマルチパートポリゴンとのオプションを有効にして同じカテゴリー値をもつものを1つのポリゴンにする。
  • このファイルをそのまま画像分類のトレーニングサンプルマネージャで読み込んでクラス名を与えても良い

もしくは

  • あらかじめトレーニングサンプルマネージャで出力されるポリゴンと同じ列名をもつうようにポリゴンを修正してから読み込み。

教師データのバンド間散布図の作成

  • [ランダムポイントの作成]で教師データポリゴンの中に点を発生させる。ポリゴン全部のピクセル値を抽出するのは計算容量上困難なためである。(リモセンソフトならできるが、ArcGISはラスタ周りに弱いので)
  • [空間結合]で、ポイントに元のポリゴンのカテゴリを付与する。
  • 複数値からポイントに値を抽出

水深補正用消散係数の計算

logをとった画像からポイントに値を付与して、砂のカテゴリのみを抽出して、各バンド間の係数を計算(Excelno直線回帰の傾き)。
係数をかけたバンド比BIを計算して、この値を教師付分類用のラスターとする。なお赤にBand2-3を、緑にBand1-2を、青にBand1-3をいれると、海草を緑で砂を赤で抽出したっぽい色合いになります。

分類結果の出力

そのままGeoTIFFで出力するので良いです。しかし、Geotiffだと真っ黒や真っ白にしか画像のビューアーで表示されないため、画像として表示したいことが多くあると思います。

色を指定しないのであれば、飛ばしても良いですが、
カラーマップ インデックス(RGB 情報)を事前に指定しておくと色が変わらないので便利です。
この実施には、整数のセル値を持つ1 バンドのラスターデータである必要があります。小数値のラスタの場合はPythonで整数値の0でも足したら良いでしょう(Pythonは記入した桁数で計算するため)。

  • ラスターデータの属性テーブル(VATテーブル)を作成。(ArcToolbox[ラスター属性テーブルの構築(Build Raster Attribute Table)])
  • 属性テーブルに[フィールドの追加]でRed、Green、Blueの3つの[Double]のフィールドを順に作成。
  • [エディタ]ツールバーで[編集の開始]をして、RGBに任意の値を入れる。この際、255を1とする割合の値を入力(50の場合は50/255=0.196078)
  • 一度レイヤをTOC(コンテンツのテーブル)から削除して再度追加すると更新される。自動で更新されない場合は、[プロパティ] の [シンボル]の[表示]部分に新たに表示された[カラーマップ]をクリック

単純に


# Replace a layer/table view name with a path to a dataset (which can be a layer file) or create the layer/table view within the script
# The following inputs are layers or table views: "finished_trace_image\Trace_hand_TRN_LYL_2009dec_v1_diss_multi"
arcpy.CreateRandomPoints_management("C:/Users/gisdesktop/Dropbox/2017_heavy_drop/20171001BioMarina用データ/Shapes1","Trace_hand_TRN_LYL_2009dec_v1_diss_multi_point10k.shp","finished_trace_image/Trace_hand_TRN_LYL_2009dec_v1_diss_multi","0 0 250 250","10000","1 Meters","POINT","0")