繁忙期に役立つ？かもしれないArcGIS Pro Pythonでの複数データの一括処理(1)
はじめに
（追記）正確には一括処理というより、逐次処理を何回も行った結果としての一括した処理結果を得ていると書くべきなのですが、ぱっと見て、意味がわかりにくいので、一括処理と書いています。
こんにちは。お久しぶりです。
12月も半ばを過ぎ、いよいよ繁忙期です。
建設コンサル系、地図情報系業務に関わっている関係で、年度末にかけて、解析系の作業はこれからが本番だったりしますが、この時期は少しでも時間を大切にしたいもの。
時間に余裕があるときと違い、条件が何度も変わったり、データの更新があったときなど、手作業でやっていたら、混乱必至になる場合があります。
コードを書いてできることはいっぱいあるのですが、ここでは、複数データの一括処理（ループなどで回す処理）に限定して、ご紹介します。
なお、ここで書いたコードは、あくまで個人が参考用に作成したものであり、公式見解ではありません。また、ソフトのバージョンで動作が変わるかもしれません。使用上の責任は負えないのでご注意ください。
コードの作成環境
作業環境は以下の通りです。
ソフト ArcGIS Pro 3.0.3
コードの作成環境 Jupyter Notebook
ここではArcGIS Proに備わっている、Jupyter Notebookを使います。
同じ操作は、ArcGIS Proに備わってるPython Windowでも可能ですし、外部のスクリプト開発環境（例：Visual Studio Code）でも可能です。
ただ、作業環境構築の手間がないのと、セルごとに実行を確認できる利便性などから、敷居は一番低いと思います。
更に言うと、ArcGIS ProではGUIで視覚的にModelを作成できるModelBuilderという優れた機能があります。こちらもかなり便利なツールで、ループ処理も可能です。
ただ、少し複雑な処理の場合、Modelを分けないといけなかったり、結局個々にプログラムを書かないといけなかったりして、可読性が減ることがあります。
そもそもジオプロセッシングツールの実行だけですむようなものもあります。
いずれを使うかは、解析の内容によって使い分けるといいと思います。
処理の概要
今回は次の処理をするコードを書きます。
処理自体はかなり特殊なことをしています。
また、他の方法でもできる部分がありますがそれは後述します。
あるフォルダ以下の複数のラスタ画像情報を逐一取得
その画像の範囲を含む長方形のポリゴン（envelope、extentなどといいます）を1つずつ作成する
ファイル名、投影座標系（地理座標系）、ラスタ画像の最小値、最大値を属性としてポリゴンに付加する
データは何でも良いのですが、コンサベーションGISコンソーシアムジャパンhttp://cgisj.jp/index.htmlで公開している「標高緯度経度」情報を対象にします。このデータの特徴と今回の保存の仕方は以下のとおりです。
地理座標系、測地系JGD2011（EPSG：6668）
ERDAS IMAGINE形式
データの範囲は原則、都府県、北海道のみは振興局単位。離島が存在している場合は離島の範囲になっている場合あり
ローカルフォルダは都道府県単位としていて、そこに1つもしくは複数のデータを保存しています
Pythonコードの例
必要に応じコメントを付記しています、ご参考まで。
一部の解説はコードのあとに書いています。
%%time
#マジックコマンド 処理時間を記載
import arcpy
import os
from arcpy import env
DemDir = r"d:\dem\imagine" #画像のあるフォルダ rをつけることで￥でも使える
env.workspace = r'd:\dem\Default.gdb' #指定しないとデフォルトのジオデータベース
# ファイルが複数フォルダ、サブフォルダにまたがっている場合はwalkを使う
# （globでも良いがwalk推奨）
walk = arcpy.da.Walk(DemDir,datatype='RasterDataset',type='IMG')
# 図形形状含めデータ属性を格納するリスト
attributes=[]
# 既にデータがあったら消す（ここはなくても良い）
if arcpy.Exists('DemExtent'):
arcpy.Delete_management('DemExtent')
# walkに含まれるファイルの情報でループを回す
for dirpath, dirnames, filenames in walk:
for filename in filenames:
rasterlink=os.path.join(dirpath, filename) #ファイルをフルパスで指定する
desc=arcpy.Describe(rasterlink) #rasterの概要情報
name=desc.name #データ名称（ここではファイル名と同じ）
srs=desc.spatialReference #データの投影（地理）座標系(この例では地理座標系、測地系JGD2011、SRID、EPSG 6668)
srsname=srs.name #データの空間参照系名（今回は、GCS_JGD_2011）
rasterObj = arcpy.Raster(rasterlink) #rasterオブジェクト生成
DemHeightMin=rasterObj.minimum #データの最小値
DemHeightMax=rasterObj.maximum #データの最大値
# 4隅の頂点座標
XMin,YMin,XMax,YMax=desc.extent.XMin,desc.extent.YMin,desc.extent.XMax,desc.extent.YMax
# 取り込むデータをリストにセット
# 頂点座標リスト、データ名称、投影（地理）座標系、最小値、最大値
rowvalue=[[(XMin,YMin),(XMax,YMin),(XMax,YMax),(XMin,YMax)],name,srsname,DemHeightMin,DemHeightMax]
# 複数あるのでリストに追加
attributes.append(rowvalue)
# データ格納用に空のフィーチャークラスを作成（ジオデータべース、名称、形状、投影（地理）座標系もセット）
arcpy.management.CreateFeatureclass(env.workspace,'DemExtent', 'POLYGON', spatial_reference=6668)
# フィーチャークラスを指定
feature_class = env.workspace+'/DemExtent'
# 以上２つはこうやっても良い
# result=arcpy.management.CreateFeatureclass(env.workspace,'DemExtent', 'POLYGON', spatial_reference=6668)
# feature_class = result[0]
# ィーチャークラスの属性フィールド名称とタイプを追加する
arcpy.AddField_management('DemExtent','DemName','TEXT')
arcpy.AddField_management('DemExtent','Projection','TEXT')
arcpy.AddField_management('DemExtent','Min','FLOAT')
arcpy.AddField_management('DemExtent','Max','FLOAT')
# 属性リストを元に図形と属性を追記していく
# SHAPE@が図形そのものを表す情報で、DemName以降がテーブルの属性
for row in attributes:
with arcpy.da.InsertCursor(feature_class,['SHAPE@','DemName','Projection','Min','Max']) as cursor:
cursor.insertRow(row) ここでキーになるのは次の5つかと思います。
walk
arcpy.Describe
arcpy.Raster
arcpy.management.CreateFeatureclass
arcpy.da.InsertCursor、cursor.insertRow
walk
あるフォルダ以下の指定タイプのデータのフォルダ、ファイル名情報を取得します。この例ではrasterdataset（ラスタデータ）、IMG（ERDAS IMAGINE形式）のデータを探しています。
arcpy.Describe
データの概要情報を入手します。ファイル名、投影（地理）座標系、画像の範囲などです。画像の範囲はその画像を含む最小の長方形（座標に平行）の4隅となります。
arcpy.Raster
walkで取得した場所にあるファイルからデータのオブジェクトを取得します。今回はその結果からデータの最小値、最大値を取得するだけですが、データ処理などに使うこともできます。
arcpy.management.CreateFeatureclass
データ格納用に空のフィーチャークラスを作成します。引数として、格納するジオデータべース、フィーチャークラスの名称、データ形状（Point、Polyline、Polygonなど）、投影（地理）座標系などを与えることができます。ジオデータベースでなくshapeファイル形式を指定することもでき、その時はフォルダとSHAPEファイル名をそれぞれ指定する必要があります。
その後、arcpy.AddField_managementで属性フィールドを追加します。データ型も含め指定します。
arcpy.da.InsertCursor、cursor.insertRow
arcpy.management.CreateFeatureclassで作った空のフィーチャークラスに１つずつフィーチャー（ここではポリゴン）を加えていきます。
まず、arcpy.da.InsertCursorで入れるデータを引数としてセットし、cursor.insertRowで実際にデータを書き込みます。
ここでは図形も引数の1つとなっていて、SHAPE@はフィーチャのジオメトリオブジェクトになります。
ジオメトリオブジェクトとしては、ポイントオブジェクトを作ってから入れる方法のほか、ここに示したように直接、座標指定だけで設定することもできます。 [(XMin,YMin),(XMax,YMin),(XMax,YMax),(XMin,YMax)]の部分です。
オブジェクトを作らなくていいので、処理は速いようです。
ただ、この方法だとドーナツポリゴンのような複雑な形状、マルチパートフィーチャーは指定できないようです。詳細は以下を参考にしてください。
その他の属性も、同じように続けて書きます。
['SHAPE@','DemName','Projection','Min','Max']
ここでは、複数のデータをリストに格納し、それを再度ループで取得しているので、カッコの入れ子が細かくなっています。注意しましょう。
1行目に記載しているのはマジックコマンドで処理時間を計測するものです。今回は95のラスタデータが対象になっているのですが、私の環境では約7秒かかっています。処理として速いかどうかは比較していないのでわかりませんが、1つ1つ手でやること、何度でも繰り返し使えることを考えれば有効ではないかと思います。
結果
結果、ポリゴンがこんな感じで生成されます。
属性も同時に生成されています。
背景は地理院地図
下の図は、これを長野県の範囲にズームしてみたものです。
長野県部分のみ青枠で強調しています。また白黒のデータは標高緯度経度ラスタデータです。県境から4kmの外側の範囲までをデータにしています。
画像の外周に沿った形で、青い長方形ができていることがわかります。
画像は標高緯度経度情報の長野県のデータ、背景は地理院地図
標高緯度経度データの作成に当たっては、国土地理院長の承認を得て、同院発行の基盤地図情報数値標高モデルデータを使用しました。（測量法に基づく国土地理院長承認（使用）R3JHs592）。
まとめ
今回はArcGIS Pro Pythonでの複数データ処理の一例を紹介しました。
初回したのはラスタデータの処理ですが、対象がshapeファイルでも、同じようなコードで応用できます。
逐一処理をする方法を知っておけば、原則的にはその後何でもできますので、知っておくといいと思います。
勿論、実際には処理速度、メモリ使用量の問題などがあるので100%万能というわけではありませんが。
実はこの結果のポリゴンはArcMapならラスターカタログ、ArcGISPROだとモザイクータセットで同じような結果を得ることができます。
詳細は省略しますが、もしある処理を複数のラスタデータに一括適用したい場合とか、もっと大量のデータがあるような場合には是非使用を検討してみてください。
そうではなくて単にデータの範囲と概要だけを見たいだけということであれば、こちらは軽い処理ですので、有効かと思います。
次回の予定としては、個々のデータファイルではなくて個々のレコード（フィーチャー）を逐一処理するような方法を紹介しようと思います。
最後に
ここまでご覧いただき、ありがとうございました。
いつもは北海道に本拠地を置くNPOに所属し、環境保全を主な題材としてGISやリモセンに関する仕事をしています。
コンサベーションGISコンソーシアムジャパン の活動もその1つです。
この分野の仕事や活動でなにかお困りのある方、ご相談ごとのある方など、是非コメントをお寄せいいただくか、上記WEBサイト掲載のメールアドレスまでメールをお寄せください。
コンサベーションGISコンソーシアムジャパンの活動及びこのnoteでの活動はボランティアで行っているのですが、何分資金が不足しています。
もしサポートしてもいいなという方がいらっしゃいましたら、よろしくお願いします。
この記事が気に入ったらサポートをしてみませんか？