ArcPy で画像処理をしてみよう～その３：ラスタ演算（Spatial Analyst）

当シリーズ：「ArcPy で画像処理をしてみよう」ではこれまで、ラスタオブジェクトの特徴やジオプロセシング関数での使い方について紹介してきましたが、いわゆる「画像処理」として、もっと数値計算的なものを期待していた方もいるかも知れません。今回はそのような期待に応えるべく、ラスタ演算について紹介します。
画像またはラスタを使って様々な演算を行うことをラスタ演算と言います。重ね合わせた複数の画像の同じ位置のピクセル同士で演算を行ったり、画像内の全ピクセルを使って統計量などを計算したり、対象ピクセルの周囲のピクセルで演算をしたりと様々なバリエーションがあります。そのような処理を Python で０から書いていくことも可能ですが、ArcGIS では、Spatial Analyst エクステンションで様々なラスタ演算を簡単に行うことができる演算子や関数を提供しています。Spatial Analyst を使用する場合でも、ラスタオブジェクトが重要な役割を果たします。画像データをラスタオブジェクトに格納し、計算結果もラスタオブジェクトに出力します。

ラスタ演算をしてみよう

今回もまずは ArcMap の Python ウィンドウで書いてみましょう。Spatial Analyst には「マップ代数演算」というラスタ演算用の電卓のようなジオプロセシングツールもありますが、Python ウィンドウの方が何より Python の体系と一緒に使えるため、断然使いやすくなります。ここでは、Spatial Analyst が得意とする DEM （ラスタ形式の地形データ）を使ったラスタ演算を行ってみましょう。

※ ArcMap メインメニューの [カスタマイズ] → [エクステンション] の [Spatial Analyst] にチェックが入っていることが前提ですので、確認しておいてください。

準備

前回行ったのと同様に、まずは適当な DEM ファイルを用意します。適当な DEM なんて持ってない！という方は以下のリンクから SRTM データをダウンロード→解凍して利用すると良いでしょう。

SRTM サンプルデータ

このデータ（.hgt ファイル）はそのまま ArcMap で読み込んだ状態で地理座標（緯度経度）を持っていますが、この状態のままでは後で使用する Slope 関数などで不都合がありますので、UTM 等の投影座標に投影変換して「dem.tif」という名前で出力したものを使用してください。

これを ArcMap に追加した状態から始めます。

Spatial Analyst の関数を使用する際には、以下のように arcpy.sa モジュールを読み込んでおけば、Spatial Analyst の関数（例：arcpy.sa.Slope() ）をいきなり関数名から（例： Slope() ）書くことができますので記述が短くなります。

>>> from arcpy.sa import *

まずは、用意した DEM ファイルをラスタオブジェクトにします。前回までは Raster() の引数としてファイルのパスを指定しましたが、ArcMap に DEM を追加している場合は、このようにラスタレイヤ名を入力に使用することもできます。

>>> inRaster = arcpy.Raster(”dem.tif”)

算術演算

それでは、簡単な例としてピクセル値を 2 倍にする計算をしてみましょう。これだけでできます。

>>> outRaster = inRaster * 2

ArcMap 上の Python ウィンドウを使用している場合は、この時点でテンポラリのラスタレイヤとして「outRaster」が自動的に追加されたはずです。（されない場合は、ArcMap メインメニューの [ジオプロセシング] → [ジオプロセシングオプション] のダイアログ下方にある [ジオプロセシング処理結果をマップに追加] を有効にしてください。）個別属性表示ツールでピクセル値を確認してみてください。ちゃんと 2 倍になっているでしょうか？

出力ラスタオブジェクトをファイルに保存するには、.save（ファイル名）を使用します。

>>> outRaster.save(r”c:tempoutput.tif”)

二値化（比較演算）

非常によく行われる画像処理である「二値化」つまり、条件によって 0 か 1 に分類する処理もとても簡単です。たとえば、ピクセル値が 1000 以上という条件で二値化するなら、このように書きます。

>>> outRaster = inRaster >= 1000

ちなみに左の「＝」は代入記号、右の「>=」は比較演算子です。どちらも数学の等号・不等号とは少し意味が異なる点に注意してください。つまり、画像 inRaster の値が 1000 以上なら 1、未満なら 0 を画像 outRaster に代入（出力）する、という式です。

このような算術演算や論理演算、論理演算の演算子として何が使用できるかは、以下を参照してください。

マップ代数演算の演算子の操作

関数の利用

続いて、Spatial Analyst 特有の地形解析の関数を使ってみましょう。

>>> outSlice = Slice(inRaster,10) #高さを 10 階層に分類した段彩図

>>> outHill = Hillshade(inRaster) #陰影起伏図

>>> outSlope = Slope(inRaster) #斜度（°）

Slope 関数と先ほどの比較演算子を応用して、DEM から斜度 30 度以上の急傾斜地を抽出する処理を１行で記述することもできます。

>>> outRaster = Slope(inRaster) >= 30

非常に簡単な例でしたが、マップ代数演算は Python ウィンドウで行った方がずっと簡単、ということは理解していただけたのではないでしょうか。

関数は ArcToolbox の中の「Spatial Analyst ツールセット」として GUI でも利用できるようになっています。どのような関数があり、どのように使用するかは以下のページから確認することができます。

Spatial Analyst ツールの完全なリスト

バンド演算をスクリプトで書いてみよう
もちろん、Landsat のような衛星画像で解析を行う場合にも ArcPy のラスタ演算を使用することができます。特に、衛星画像のバンド同士で演算を行うバンド演算は、元の画像上では見えにくい地表の特徴を強調するためによく使われます。

バンド演算

ここでは、衛星画像の「赤」バンドと「近赤外」バンドを使用した正規化植生指数（NDVI）の計算を行うスクリプトを書いてみましょう。衛星画像がない方はブログ記事「無料でダウンロードできる Landsat 画像の活用」を参考に「コンポジットバンド」まで行った状態の Landsat 画像を準備してください。

Python スクリプトで書く場合には、まず arcpy.CheckOutExtension(“Spatial”) により Spatial Analyst エクステンションのライセンスを取得する必要がある点が Python ウィンドウの場合と異なります。

スクリプトの説明

（ア）必要なモジュールをインポートします。

（イ） Spatial Analyst ライセンスを取得します。

（ウ）入出力ファイルを定義します。

（エ）マルチバンドの入力ファイルをラスタオブジェクトに格納し、さらに arcpy.Describe().children を使用して、取得したバンドの情報を bands に格納します。バンド名の取得の仕方はこのようにややトリッキーになります。

（オ） inRaster と bands からバンド３とバンド４それぞれの全パス名を生成し、ラスタオブジェクトにします。

（カ）ラスタオブジェクト化されたバンド３とバンド４を使って NDVI を計算します。Float 関数を使用するのは、ゼロで割るのを回避するためです。

（キ）出力ファイルに保存して終了です。

NDVI は -1 ～ 1 の値をとるデータで、値が大きいほど植物の活性度が高いという指標です。さらにこれを適当な閾値で二値化すれば、緑地の抽出に応用することができます。この記事の前半を思い出して（カ）の末尾に >0.5 などと追加するだけですね。

いかがでしたか？今回は Spatial Analyst のラスタ演算機能を ArcPy で使用してみましたが、意外に簡単だったのではないでしょうか。今回は単純なものだけでしたが、Spatial Analyst に含まれる多変量解析や水文解析、日射量解析などの高度なツールも Python の関数として使用できます。そういった複雑なロジックを自分で組まなくても利用できるのは大きな魅力と言えるでしょう。

しかしながら、自分で考えたアルゴリズムなどを実装するには不便な面もあります。ラスタオブジェクトの演算は画像という単位で処理をしてくれて非常に楽なのですが、通常の画像処理のプログラミングで行うような個々のピクセルの値へのアクセスはできません。そこで次回は、本シリーズの最後として、「NumPyArray」を利用した画像処理を紹介したいと思います。

[注意事項]

本記事に関するご質問は、保守 Q&A サポートサービス以外の有償サービスにて対応いたします。本記事に関して生じたいかなる損害についても、弊社では責任を負いかねますことを予めご了承願います。

本スクリプトツールは、サンプルツールとして提供しておりますので、すべてのお客さまの環境での動作を保障するものではありません。また、スクリプトの内容に関するご質問はお受けできませんので、ご了承ください。
なお、本サンプルを使用して生じたいかなる障害についても、弊社では責任を負いかねますことを予めご了承願います。

弊社サポートサイトには多くのサンプルコードが掲載されております。 Python を使用したスクリプトツールを作成される際は、是非弊社サポートサイトをご活用ください（サポートサイトのご利用はサポートページへのログインが必要です）。