シリーズでお送りしてきた「ArcPy で画像処理をしてみよう」もいよいよ最終回となりました。当シリーズではこれまでその１「ラスタオブジェクトとプロパティ」でラスタオブジェクト、その２「ジオプロセシング」でジオプロセシング、その３「ラスタ演算（Spatial Analyst）」でラスタ演算とステップを踏んできましたので、最後は最も Python らしい世界へご案内したいと思います。
すでに ArcPy をお使いの皆様の中には、オンラインヘルプの [ジオプロセシング] → [ArcPy] → [ArcPy 関数] → [Raster] の下に「NumPyArrayToRaster」「RasterToNumPyArray」という名前の関数があるのに気付いている方もいらっしゃるかも知れません。これらは その名の通り、NumPyArray とラスタオブジェクトの間の変換を行う関数です。では「NumPyArray」 とは？
NumPyArray とは？
ArcGIS のバージョン 9.3 以降では、NumPy というモジュールが、一緒にインストールされるようになっています。NumPy は数値計算用の Python モジュールで、科学技術分野では比較的よく使われているものです。その中核とも言えるデータ構造が Array（アレイ）、つまり「配列」です。簡単に言うと、複数の数値をひとまとめにして効率よく取り扱うための仕組みです。大量のデータを含む画像や時系列データなどを扱う際に活躍します。NumPy モジュールの中で作成・使用する配列なので、NumPyArray というわけです。
NumPyArray 自体は単純な配列ですので、ArcGIS ／ Python の中でもいろいろな場面で使用できますが、今回もやはりラスタの話ですので、2次元の配列を主に取り扱います。そもそも画像／ラスタとは2次元の配列です。これまで紹介してきたラスタオブジェクトやそれを利用するツールでは、ユーザはそういうことを意識する必要がないようになっているのですが、画像をより原始的に（数値データの集合として）取り扱いたい場合や、他のモジュールと連携を図る場合に配列として取り扱う必要がある場合があります。
NumPy についての詳細は以下のサイトなどを参照してください。（これ以外にも日本語サイトも多数存在します。）
Numpy and Scipy Documentation
RasterToNumPyArray 関数を使ってみよう
前置きが長くなりましたが、今回も ArcMap メニューの [ジオプロセシング] → [Python] から、Python ウィンドウを起動してください。また前回同様、SRTM から作成した dem.tif を ArcMap に追加しておきます。
まず Python ウィンドウにて、以下の記述でラスタオブジェクトに格納します。
>>> inRaster = arcpy.Raster(“dem.tif”)
続いて RasterToNumPyArray という、ラスタオブジェクトを NumPyArray に変換する関数を使用して、本日の主役である NumPyArray を生成します。
>>> inArray = arcpy.RasterToNumPyArray(inRaster)
試しに NumPyArray の shape プロパティ（配列の縦横のサイズ）を取得してみましょう。
>>> inArray.shape
(1260, 1041)
配列ですから、前から何番目という番号（インデックス）を指定すれば、そこに入っている値にアクセスすることができます。たとえば、画像の左上を原点（0,0）として、下方向へ 800 、右方向へ 750 の位置にあるピクセルの値を表示してみましょう。
>>> inArray[800,750]
3715
ちなみにこれは、富士山の頂上付近のピクセルです。
いくつかのピクセルをまとめて取得することもできます。例えば、先ほどの「下方向へ 800」の代わりに、「700 ～ 900」と幅を持たせてピクセルを取得します。
>>> inArray[700:900,750]
array([1179, 1189, 1196,（・・略・・） 1139, 1133], dtype=int16)
このように大量のピクセル値が取得されたはずです。こうして取得された値の集合も NumPyArray になります。
・PyLab モジュールによる可視化
せっかくですので、これをグラフにプロットしてみましょう。NumPyArray を使用する利点として、当の「NumPy」はもちろん、データを可視化するモジュール「PyLab」と連携できる点があります（ArcGIS 10.1 からは PyLab も同時にインストールされるようになりました）。まずはこれらをインポートします。np、pl のような別名を指定する文献が多いようですので、ここでもそれに従います。
>>> import numpy as np
>>> import pylab as pl
先ほどと同様に、inArray の一部を取り出して y という名前で格納します。これがグラフの Y 軸になります。
>>> y = inArray[700:900,750]
グラフの X 軸には整数を順に並べた配列を使用します。NumPy の arange 関数の引数に配列 y の長さ（値の個数）を入れます。
>>> x = np.arange(len(y))
X 軸と Y 軸のデータが揃ったので、pylab の plot 関数でプロットのオブジェクトを生成します。
>>> pl.plot(x,y)
[<matplotlib.lines.Line2d object at 0x25CA3D90>]
最後にオブジェクトを表示する show 関数を実行します。
>>> pl.show()
以下のような富士山の断面図が描けたでしょうか。
PyLab では断面図のほかにコンターや散布図など様々な形態の図を作成することができます。PyLab（matplotlib） については、以下のページが参考になります。
http://matplotlib.org/index.html
NumPyArrayToRaster 関数を使ってみよう
続いて、NumPyArray をラスタオブジェクトに変換する関数、NumPyArrayToRaster() を使って、2 次元配列からラスタデータを作成してみましょう。
まず、先ほど作成した inArray に対して何か演算を行います。例えばセル値を 2 倍にする計算をする場合、以下のような配列計算で簡単に書くこともできますし、
>>> outArray = inArray * 2
以下のような二重の for 文で個々のセル値について演算を行うこともできます。ここでは入力画像のピクセル値を 2 倍する計算を全ピクセルについて繰り返すだけなのであまり意味がありませんが、各ピクセルの近傍ピクセルを使った演算など、個々のピクセルにアクセスするような処理ではこのような書き方も役立ちます。
>>> outArray = np.zeros(inArray.shape) #出力用配列を作成（初期値０）
>>> for i in np.arange(inArray.shape[0]):
… for j in np.arange(inArray.shape[1]):
… outArray[i,j] = inArray[i,j] * 2
…
outArray をラスタオブジェクトにします。
>>> outRaster = arcpy.NumPyArrayToRaster(outArray)
ラスタオブジェクトになったので、ArcMap のビューに表示されたはずですが、元のデータとは重ならず、はるか遠い位置に表示されてしまいます。作成されたラスタオブジェクトには地図座標が定義されていないからです。NumPyArray は単なる数値の集合ですので、地図上の位置の情報は別途与えてやる必要があるのです。NumPyArrayToRaster 関数は入力配列以外に、画像の左下の位置と横縦のセルサイズ、NoData 値を引数として与えることができます。ここでは入力配列の元になったラスタオブジェクト inRaster の情報を取得してそのまま使用します。
まずは inRaster の左下隅の座標を与えるポイントオブジェクトを作成します。
>>> llc = arcpy.Point(inRaster.extent.XMin, inRaster.extent.YMin)
セルサイズと空間参照も inRaster から取得します。
>>> cx = inRaster.meanCellWidth
>>> cy = inRaster.meanCellHeight
>>> sr = inRaster.spatialReference
左下隅座標ポイントとセルサイズを NumPyArrayToRaster の引数に追加します。
>>> outRaster = arcpy.NumPyArrayToRaster(outArray, llc, cx, cy)
空間参照は別途 DefineProjection_management 関数で与えます。
>>> arcpy.DefineProjection_management(outRaster, sr)
最後にファイル出力します。
>>> outRaster.save(r”c:tempdem_out.tif”)
いかがでしたか？ここでは 「NumPyArrayToRaster」「RasterToNumPyArray」 関数を紹介しましたが、NumPyArray は arcpy.da モジュールの関数でフィーチャクラスやテーブルとの間の相互変換も可能です（関数名は Raster を FeatureClass や Table に置き換えたもの）。GIS データを NumPyArray にして NumPy の機能による演算や PyLab の様々な可視化ライブラリを利用することも可能になりますし、逆に様々な数値計算の結果の NumPyArray を ArcMap で使用できる GIS データとして出力することができるのです。
＊＊＊＊
さて、シリーズでお送りした「ArcPy で画像処理をしてみよう」はこれで一旦終了です。近年 GIS やリモートセンシング系のソフトウェアでは、それらの機能を Python モジュールとして利用できるようになってきました。ArcGIS 内部のツールや Python の標準モジュールのみならず、そういった外部のソフトウェアの機能とも連携して、ジオプロセシング ツールの開発や業務の効率化にお役立て頂ければ幸いです。
[注意事項]
本記事に関するご質問は、保守 Q&A サポート サービス以外の有償サービスにて対応いたします。 本記事に関して生じたいかなる損害についても、弊社では責任を負いかねますことを予めご了承願います。
本スクリプトツールは、サンプルツールとして提供しておりますので、すべてのお客さまの環境での動作を保障するものではありません。また、スクリプトの内容に関するご質問はお受けできませんので、ご了承ください。
なお、本サンプルを使用して生じたいかなる障害についても、弊社では責任を負いかねますことを予めご了承願います。
弊社サポート サイトには多くのサンプル コードが掲載されております。 Python を使用したスクリプト ツールを作成される際は、是非弊社サポート サイトをご活用ください（サポート サイトのご利用はサポート ページへのログインが必要です）。
■関連リンク
・トレーニング【Pythonによるジオプロセシング スクリプト入門】
・第9回 GISコミュニティフォーラム Python 関連セッション：「初心者でも簡単にできる！Python スクリプトの作成と共有方法」（プレフォーラム・セミナー 5月29日（水）16：00-16：40）