空撮画像から作成される農場の 3 次元モデルを分析して収量予測や雑草の識別を行う場合，モデルに含まれる森林，家屋，倉庫，作業道路などを除外する必要がある．本研究では，農場の 3 次元モデルから草地圃場のみを選択するための処理を紹介する．モデル作成用ソフトウェアに備わっている点群分類機能と数値表層モデルの利用を前提とし，プログラム言語 Python によって選択処理をコンピューターへ実装する．処理は 3 段階からなる：1） 点群分類の実施による森林と建築物の大部分を除いた数値表層モデルとオルソモザイク画像の構築，2） 数値表層モデルの傾斜属性の利用による森林部分に残存するノイズの完全な除去，3） 作業道路の色情報の利用による解析対象外の圃場等の除去．1） を実施するために Agisoft Metashape Professional を利用し，2） と 3） を実施するために Python スクリプトを作成した．その際，処理を視覚的に実行するためのユーザーインターフェースも組み入れた．異なる2 つの草地で 3 つの時期に 3 種の小型マルチコプターを用いて撮影された空撮画像においても本処理は有効であり，これら一連の画像処理を利用すれば不定形の草地圃場を選択することも容易である．
処理
1．点群分類機能を有する3 次元モデル構築ソフトウェアの利用
写真測量に基づいて 3 次元モデルを構築するためのソフトウェアは数多く存在する（表 1）．ほとんどの有償品の価格は数十万円と高価なので，個人的な購入は困難かもしれない．無償品はインターネットを介して入手可能である．なお，WebODM は安価に導入できるし，手動での導入およびオフラインでの利用ならば無償である．著者らは表中のソフトウェアのごく一部しか試用できていないが，点群分類機能を有するソフトウェアとして，WebODM，Metashape Professional および Pix4Dmapper を確認した．本研究では Metashape Professional 1.7.2 を利用する．本ソフトウェアは，特定の分類群を除去した後の 3 次元モデル構築機能も有する．
本研究において点群分類を利用する目的は，農場画像から森林部分をおおまかに除去することである．高密度点群データでは森林部分は他に比べて高さが突出する．分類が可能であれば，その部分を含めずに 3 次元モデルを構築できる．以下において北海道内の酪農家が所有する放牧草地（芦別市，ペレニアルライグラス（Lolium perenne L.），ケンタッキーブルーグラス（Poa pratensis L. subsp. pratensis），シロクローバ（Trifolium repens L.）混生）の画像に対する処理を示す．8 月下旬に放牧草地全体を含むように小型マルチコプター Phantom 4（DJI 社，Shenzhen，China）を使用して撮影した（渡辺 2018）．撮影高度は約 90 m で，撮影面積は約 75 ha である．まず，空撮写真から Metashape Professional を使用して高密度点群データを得る（図 1A）．このデータには畑，建物，路面，および画面に収まらない草地等が含まれるが，選択対象をそれらからほぼ東側にある放牧草地とした．そこは森林に区切られた三圃場からなり，面積は約 25 ha である．
得られた点群データに対して同ソフトウェアにより分類を実施する．メニューから［ツール］→［高密度クラウド］と進むと分類方法として［ポイントの分類］と［地上ポイントの分類］の二種類が用意されている．［ポイントの分類］の分類対象は，地面（分類後に茶色で表示される），高植生（緑色），建物（赤色），路面（黒色），車（水色），人工物（薄橙色）とされており，それらを個別に選択できる．これらすべての対象を選択した場合の分類結果を図 1B1 に示す．一方，［地上ポイントの分類］は点群から地面のみを分類することに相当する．こちらを選ぶと分類のための条件設定ウィンドウが開き，地表面の最大傾斜角，地表面の最大変化高度，点群を分割するセルサイズを入力できる．初期設定値はそれぞれ 15.0 度，1 m，50 m である．その条件における分類結果を図 1B2 に示す．
［ポイントの分類］の結果は完全ではないものの良好である．森林部のほとんどが高植生を示す緑色として区別され，草地内の立木も緑色で表示された．草地は大部分が地面を示す茶色で示されたが，建物を示す赤色や路面を示す黒色で表示された部分もあった．放牧草地内の植被の少ない箇所が建物や路面として分類されたものである．路面が黒色でほぼ正確に分類され，建物は緑色で表示される箇所が誤分類であった．他方，［地上ポイントの分類］の結果は非常に良好である．草地や路面は地面を示す茶色で，森林や建物は灰色で表示された．
次に森林部を除去した DSM（図 1C）とオルソモザイク画像を構築する（図 1D）．［ポイントの分類］を利用した構築の際には，除去したい対象をいろいろな組み合わせで選択できる．ここでは，地面（≒草地）として分類された部分のみを選択した場合（図 1C1，D1），および高植生（≒森林）として分類された部分を除いた場合（図1C2，D2）に構築されるモデルを示す．また，［地上ポイントの分類］を利用し，高植生を除いた場合に構築されるモデルを図1C3，D3 に示す．
どの場合でも，DSM とオルソモザイク画像から森林部分および放牧草地内の立木の大部分を除去できたが，森林であった部分に除去できないノイズが多数残存する．注目すべきは，建物や路面の除去についてである．
点群モデルにおいては主要な建物は屋根の色である水色で，路面は灰色で判別できる．［ポイントの分類］に基づくと画像からそれらを除去できる場合（図 1D1）と除去できない場合（図 1D2）がある．しかし，除去できる場合では放牧草地部分の植被が少ない箇所まで除去してしまう．これは，点群分類時にその箇所を建物や路面として誤分類するからである．一方，［地上ポイントの分類］に基づくと画像から建物は除去できるが路面は除去できない（図 1D3）．草地植生の画像解析の際には人工物等は除去されていることが望ましいと考えられる．植被のある場合は，草地以外の部分をほとんど除去できる［ポイントの分類］が有用で，その後の処理を容易とする．植被の少ない状態（裸地や刈取り直後）の草地の場合は，誤分類がありうる［ポイントの分類］ではなく，［地上ポイントの分類］を利用したほうが安全である．以下では，どちらの場合でも対応できるように［地上ポイントの分類］を利用して構築したモデルを扱う．
表 2 に［地上ポイントの分類］を含めた Metashape Professional バッチファイルを示す．コード群は Extensible Markup Language で書かれているので，任意のファイル名に拡張子「.xlm」を付して保存しておけば 3 次元モデル構築時に読み込み利用できる．分類に関するコードは 8 行目である．注意すべきは 11 行目と 18 行目で，DSM およびオルソモザイク画像構築時に内挿補間を無効とする．内挿補間が有効であると除去したい部分が近傍の色で埋められてしまうので，それを防ぐためである．構築されたモデルは任意のファイル名を付して保存する．保存の際の設定は以下のとおりである．まず，DSM とオルソモザイク画像の地上解像度を統一する．また，DSM の保存の際の設定で，パレットを「グレースケール」に，モードを「傾斜」にする．［ツール］→［設定］の［外観］タブを選び［オルソビュー］画面内で設定しておけば，保存時にそれらが反映される．さらに，保存時オプションのラスター変換にパレットを，背景に白色を選ぶ．これらはバッチ処理においても設定できる．
2．点群分類後オルソモザイク画像からの対象草地の選択
以上までで農場全体の画像から森林部分のほとんどを除去することができるが，放牧草地以外の部分を除去するためにさらに処理を要する．まず，森林部分に残存するノイズを除去するとともに画像内から最大の領域を選択する．そのために DSM を利用して不要部分を除くためのマスクを作成し，オルソモザイク画像へ重ねる．DSM は傾斜という属性で出力させたので，草地内牧草よりも草高（樹高）が高い林縁や路面脇の植生は明るいエッジとして表出し検出が可能である．次いで，オルソモザイク画像の路面を利用して余分な圃場のマスクを作成し，元のオルソモザイク画像へ重ねることで目的とする放牧草地三圃場を選択する．路面は植生と比較して青色域と赤色域の反射率が高い（例えば亀山 2002）ので，どちらかを検出し画像内の圃場間の切りしろとして利用する．
コンピューターへの実装として Python によるスクリプトを作成した（表 3）．スクリプトは二つの部分からなる．一つは画像処理のための部分（2，3，5-26 行目）で，コンピュータービジョンライブラリー OpenCV と数値計算ライブラリー NumPy を利用する（3 行目）．マスク作成時の変数 thr_1，cls，thr_2（22，25 行目）に直接数値を入力するのであれば，この部分で処理を完結できる．しかし，適切な値を決定するためには上記変数，特に thr_1 とthr_2 に多数の値を試行しなければならず煩雑である．煩雑緩和のためのコードが残りの部分（1，4，27-37 行目）で，グラフ描画ライブラリー matplotlib と対話型ユーザーインターフェース作成ライブラリー ipywidgets を利用する（1 行目）．閾値決定のためのスライダー等インターフェースはデータ分析ツール Jupyter Notebook ないしJupyterLab 上での表示を前提としているので，スクリプト全体を分析ツールのどちらかで実行する．
本スクリプトは計算環境としてUNIX 互換オペレーティングシステムを想定しており，Debian GNU/Linux 10.9（Debian Project 2021）上で作成した．Microsoft Windows（Microsoft 社，Redmond，Washington）上で作成・改変・実行・保存するときは，画像ファイルへのアクセスのためのコード（17，20，26 行目）先頭の「/」を「C:¥」のように記憶領域に付与された英文字を追加した形式にする．なお，ファイルへのアクセスのためにはそれらコードを相対パスで記述してもよい．また，画像ファイル名の拡張子「.xxx」は処理する画像形式に相応するものとする．スクリプト作成時の各ソフトウェアやライブラリーのバージョンは Python 3.7.0，NumPy 1.20.1，OpenCV3.4.1， matplotlib 3.3.4，ipywidgets 7.6.3，Jupyter Notebook 6.3.0，JupyterLab 3.0.14 である．これらの各オペレーティングシステムへの導入法についてはインターネット上に情報が溢れているが，著者らは Anaconda（Anaconda 社，Austin，Texas）の利用を推奨する．この方法によれば，実行環境のインストール先は利用者領域（ユーザーランド）となるので，システム環境を乱すことはない．実行環境が不要となった場合やインストールに失敗した場合は，該当のフォルダを削除するだけで済む．
1）DSM 画像とオルソモザイク画像のサイズの統一
以下においては DSM 画像とオルソモザイク画像という二種類を適宜処理するが，それぞれの画像サイズが異なると重ね合わせ処理等が不能となる．そのような処理の実態は行列の演算なので，サイズが 1 画素でも異なると演算が不可能となるからである．双方の画像サイズを同一にするため OpenCV と NumPy を利用した．表 3 の 3 行目がライブラリーの読み込み，17-21 行目が処理のためのコードである．なお，OpenCV で処理できる画像サイズの最大値は，初期設定では230 画素である．コンピューターに十分な容量のメモリーを搭載していて，初期設定値を超える画素数の画像を処理可能とするためには，2 行目のコード後半部 OPENCV_IO_MAX_PIXELS に240 画素を意味する pow（2，40）を設定する．
Metashape Professional において DSM はその周囲に特定幅の白色ないし黒色背景が付加されて出力される．それを除くためのコードが 17-19 行目である．白色背景付きの DSM 画像の場合，グレースケールとして画像ファイルを読み込む（17 行目）．画像の光度を反転し，非黒色部分を識別し，非黒色部分画像を含む最大方形の座標を求める（18 行目）．その方形サイズで DSM 画像を切り取る（19 行目）．これら処理で DSM 画像のサイズがオルソモザイク画像のサイズと近くなるが，一般的に数画素ほど DSM 画像の方が小さくなる．それを補正するためのコードが 20，21 行目である．なお，図 1 の例では，余白を除去された DSM の画像サイズがオルソモザイク画像のサイズより幅方向で1 画素分小さかった．まず，オルソモザイク画像を読み込む（20 行目）．オルソモザイク画像のサイズ情報を取得し，切り取られた DSM 画像（19 行目）のサイズを揃える（21 行目）ことで，オルソモザイク画像と同サイズの DSM 画像を得る（図 2A1）．
2）不要部分をマスクするための DSM 画像の利用
DSM 画像を利用すれば，森林部分に残存するノイズをマスクするだけでなく林縁部や路面脇の比較的高い植生を利用して目的とする草地以外の部分をできるだけマスクすることもできる．残存ノイズをマスクするだけであればオルソモザイク画像のみを利用しても可能であるが，草地に近接するノイズや目的外の圃場が残存する場合が多い．今回利用する DSM 画像は傾斜をグレースケールで表現するので，草地植生とそれより高い植生ないし低い裸地等との境界は比較的明るくなる．そのような明るい部分とその他を区分することで有効なマスクを作成する．この処理のため OpenCV と NumPy を利用した．表 3 の 3 行目がライブラリーの読み込み，5-16 行目と 22-24 行目が処理に相当する．
まず，21 行目で得られた DSM 画像の光度を反転し，特定の閾値で二値化することで引数 bin としての白黒画像を得る（22 行目）．これら処理で草地部分は白色となる．光度を反転する理由は，この後の領域選択という処理が黒色部からの白色部選択という実装になっているからである．なお，閾値の初期値として 210 を設定した（31 行目）が，以下で得られる結果が不十分であればスライダーによって変更する．さらに，圃場内に残存するノイズを減らすためにクロージング処理を行う（23 行目）．
白黒画像を関数 mask(bin) へ渡し（5 行目）処理を続ける．同一色画素を結ぶ輪郭を求める（6 行目）．輪郭によって画像の白色部分は閉じた領域となるので，小さなノイズから大きな草地まで相当数の領域が浮かび上がる．それらの面積を把握し（7 行目），最大面積の領域を選択する（8 行目）．元の白黒画像を BGR（青，緑，赤）色空間へ変換し（9 行目，OpenCV では色の配列データが RGB ではなく BGR の順番で格納される），そこへ最大領域を持つ輪郭を描くとともにその内部を赤色に塗りつぶす（10 行目，図 2B1）．さらに，オルソモザイク画像から赤色部に相当する領域を選択するための処理を続ける．BGR 画像を HSV（色相，彩度，明度）色空間へ変換してから（11 行目）赤色部のみを選択する（12-14 行目）．赤色部以外をマスクに割当てるため白黒画像とし，BGR 色空間へ変換する（15 行目，図 2C1）．これは，オルソモザイク画像（21 行目）と同じ色空間をもつ行列とするためである．オルソモザイク画像とマスク画像とを合成し，それを戻り値とする（16 行目）．マスク部分は白色を示す 255 の値を，それ以外は黒色を示すゼロの値を持つので，オルソモザイク画像と合成するとマスク部分は結局 255 の白色を，それ以外はオルソモザイク画像の色を示すこととなる（図 2D1）．以上によって，森林部分に残存したノイズや西側の畑および南側の草地等を除去したオルソモザイク画像を得る．この画像を次の処理のために格納する（24 行目）．この段階で解析のために十分な画像が得られるのであれば，以下の処理は不要である．
3）不要部分をマスクするための路面の色情報の利用
図 2D1 には西側に解析から除きたい草地圃場部分が残存する．これは道路によって放牧草地とつながっているためにマスクできなかった部分である．以下は，路面の色情報を利用してこの部分をマスクするための処理である．この処理は既述と類似で，利用する画像と色情報が異なるだけである．表 3 の3 行目がライブラリーの読み込み，5-16 行目および 25，26 行目が処理に相当する．
まず，図 2D1 の画像について BGR チャンネルのうち単一色域のみを選択する（25 行目の grs_1[:, :, cls]，図 2A2）．デフォルトとして青色域を設定したが，マスク作成に有効な色域が赤である場合も想定し選択ボックスをインターフェースに追加した（33 行目）．この画像の光度を反転し，特定の閾値で二値化することで引数 bin としての白黒画像を得る（25 行目）．これら処理で路面部分は白色となる．閾値の初期値として 130 を設定した（32 行目）が，得られる結果が不十分であればスライダーによって変更する．以後の処理は既述と同様で，5-16 行目の関数が引き受ける．最大領域の選択と赤色化（図 2B2）およびマスク作成（図 2C2）を経て目的のオルソモザイク画像を得る（図 2D2）．以上で道路を挟んだ西側の草地を除くことができた．戻り値としてのオルソモザイク画像はさらなる解析のために格納し，必要に応じて適切なフォルダへ保存する（26 行目）．
4）閾値等決定のためのユーザーインターフェース
処理過程において，DSM を利用したマスク作成時の二値化のための閾値，路面の色情報を利用したマスク作成時の色域と二値化のための閾値を決定しなければならない．これらのための matplotlib と ipywidgets を利用してユーザーインターフェースを実装した（図 3）．表 3 においてこれらライブラリーを読み込むコードは 1 行目である．インターフェースは操作部分と結果部分とからなる．コードとして前者は表 3 の 31-37 行目が相当し，スライダーや選択ボックスの定義（31-34 行目），それらと各閾値等変数の対応（35 行目），各ウィジットの配置と表示（36，37 行目）である．後者は 27-30 行目が相当し，図とフォントそれぞれの大きさ（27 行目），表示する図のリスト（28 行目），図の配置と表示（29，30 行目）である．分析ツール上のセル内に表 3 のスクリプトすべてを記述・実行すると結果がセル下部に表示される．その後，スライダー等を操作すれば再計算結果が下部の画像へ反映される．結果表示までの時間はコンピューターの能力や処理する画像のサイズに依存するが，著者らの計算環境においては数秒から数十秒までの範囲であった．なお，表示画像として処理前のオルソモザイク画像も含めた（図 3 左）が，不要であれば 27-30 行目のコードを書き直す．
DSM 利用によるマスク作成のための閾値は，スライダー Threshold_1 をマウスあるいはキーボードの左右矢印で動かし画像を見ながら決定する．このスライダーを動かすと中央および右の画像が変化する．この閾値を決定後，Channel ボックスで色域を選択しスライダー Threshold_2 を動かすと右の画像が変化する．こちらは路面の色情報利用によるマスク作成のための操作で，画像を見ながら対話的に閾値を決定する．これらの処理で双方の閾値が反映された画像が保存される．
3．他の草地圃場画像への適用結果
本スクリプトの他の異なる 2 圃場での 3 つの時期におけるオルソモザイク画像に対する有効性を検討した（図 4）．市販の小型マルチコプターに搭載されているカメラは通常 RGB 画像を撮影する．そのため，構築されるオルソモザイク画像の色調は圃場，利用するカメラ，あるいは空撮時の天候および牧草の高さや密度などに依存して変動する．ここでは，異なる色調のモデル 6 種についてスクリプトの適用結果を述べる．なお，検討はすべて単一圃場の画像について行った．実際の空撮調査では草地圃場は個々に撮影されることが多いと考えられるからである．また，供試圃場は北海道内の法人および酪農家が所有する採草地で，圃場 1 は音更町のオーチャードグラス（Dactylis glomerata L.）主体草地，圃場 2 は芦別市のチモシー（Phleum pratense L.）主体草地であった．撮影に用いた小型マルチコプターは， Phantom 4，Phantom 4 Pro，またはfMavic Pro ( 共にDJI 社，Shenzhen，China)で，異なる 3 つの時期に撮影した．撮影高度は約 90 m で，撮影面積は圃場 1 で約 34 ha，圃場 2 で約 9 ha であった．
6 種の高密度点群（図 4A1-A6）から点群分類で森林部や人工物を大まかに除いたオルソモザイク画像（図 4B1-B6）に対して本スクリプトで画像処理したところ，森林部の残存ノイズのマスクについてはすべての圃場画像に対して良好な結果を得ることができ，路面で接続する草地や地面以外はほぼ完全にマスクされた（図 4C1-C6）．圃場 2 の上から三番目の画像（図 4C6）に関してはこの段階で対象圃場のみを選択できたので，次の処理は必須ではない．マスク作成時に設定した閾値の範囲は 213-235 であり，初期値 210 付近に分布した．路面の青色域反射を利用した対象草地以外の部分をマスクする処理に関しても，検討した圃場画像で良好な結果を得た．このように異なる圃場，撮影時期，撮影機種においても，本スクリプトを用いて目的とする草地圃場のみをほぼ完全な形で選択することができた（図 4D1-D6）．マスク作成時に設定した閾値の範囲は 117-237 であり，立体情報の利用に比べ広範囲の数値に分布した．
