繁忙期に役立つ？かもしれないArcGIS Pro Pythonでの複数データの一括処理(2)
はじめに
こんにちは。
前回は、フォルダ内にある複数のデータファイルを一括で処理するという方法を説明しました。ラスタデータを対象にしていましたが、少しパラメーターを調整するだけで、shapeファイルなどの処理にも応用できるものになります。
さて、今回はファイルではなくて、1つのデータの中にある各要素、1つのポイントだったり、１つのポリゴンだったり、に順次当たりながら処理する方法について説明します。ArcGISだと各要素はfeature（フィーチャ）と言うので、ここでもfeatureということにします。ここで1つのデータ（ベクトルレイヤとして使われることが多い）に相当するのはfeatureclassということになります。featureclassやfeatureの詳細についてはこちら。
処理の概要
今回は次の処理をするコードを書きます。
ラスタ画像の範囲ポリゴン、１つずつに対して、重なる都道府県ポリゴンを調べ、その数を求めて、テーブルに入力する（そのポリゴン自体の所属都道府県１を含む）
結果の有効性は、とりあえず脇に置いておいてください。そもそも都道府県単位（＋4kmバッファー）で作ったラスタ画像なので。
基本的に上記で十分ですが、前後に次の内容を付加します。それぞれ一回限りで、手動でやれば良い処理なので、なくても良いのですが、一括処理の参考になるので付加したものです。
表示レイヤを全てオフにし、処理が終わったら全部再度オンにする
数を入力するためのフィールドを（もしなければ）追加する
テーブルの内容を全て取り出す
データとしては前回作成した、ラスタ画像の範囲を示した長方形のポリゴン（以降DemExtent）と、コンサベーションGISコンソーシアムジャパンhttp://cgisj.jp/index.htmlで公開している「行政界2021年版」情報を使用します。
行政界2021年版は市区町村単位のポリゴンになっているのですが、都道府県単位でDissolve（集約）し、マルチパートフィーチャーとしたもの（以降、PrefecturePolygon）を使用します。Dissolve方法については関連サイトなどをご覧ください。いずれのポリゴンもdefaultのgeodatabase内にあるものとします。
Pythonコードの例
必要に応じコメントを付記しています、ご参考まで。
一部の解説はコードのあとに書きます。
%%time
import arcpy
from arcpy import env
env.workspace = r'd:\dem\Default.gdb''
# layer表示を一旦全部オフにする
# map表示時に起こり得る遅延の問題を解消する
# 処理には関係ないのでなくても良い
arcgisProj = arcpy.mp.ArcGISProject('CURRENT') #現在開いているproject
mapObj = arcgisProj.listMaps('マップ')[0] #「マップ」と言う名のついているmap object
layers = mapObj.listLayers() #map内の全layerのリスト
# layer表示を1つずつオフにする
for layer in layers:
layer.visible = False
# field一覧を見る
# 重なった都道府県の数を格納するためのfieldがあるかどうか調べて
# もしなかったら足す
for field in arcpy.ListFields('DemExtent'):
if field!='Prefcount':
arcpy.AddField_management('DemExtent','PrefCount','SHORT')
# featureを参照用として１つずつあたる
# 参照して、属性を表示する（内容的にはテーブルの内容と同じ）
# shape@XYは図形の重心を得ることができる
with arcpy.da.SearchCursor('DemExtent',['DemName','Projection','Min','Max','PrefCount','shape@XY']) as cursor:
for row in cursor:
print(f'ファイル名: {row[0]}, 空間参照系: {row[1]},標高最小: {row[2]},標高最大; {row[3]},都道府県数: {row[4]},重心: {row[5]}')
# featureを更新用に１つずつあたる
with arcpy.da.UpdateCursor('DemExtent',['SHAPE@','DemName','PrefCount']) as cursor:
for row in cursor:
featuregeometry=row[0]
selectedpolygons=arcpy.SelectLayerByLocation_management('PrefecturePolygon',"INTERSECT",featuregeometry,selection_type="NEW_SELECTION")
count = selectedpolygons[2]
# 選択結果のlayerのリストに含まれる内容
# 0:選択が適用されたlayer,選択が適用された状態の更新された入力。
# 1 出力layer名,更新された入力の名前
# 2 選択数
row[2]=count
cursor.updateRow(row) #そのfeatureデータの更新
# 都道府県のfeatureが選択状態になっているので解除しておく
arcpy.SelectLayerByAttribute_management('PrefecturePolygon',"CLEAR_SELECTION")
# layer表示を再度全部オンにする
for layer in layers:
layer.visible = True構文としてキーになるのは次の6つかと思います。
今回のコードでは相互に関わっている部分があるので、構文からでなく、処理の内容から説明します。
arcpy.mp.ArcGISProject、arcgisProj.listMaps、mapObj.listLayers
arcpy.ListFields
arcpy.da.SearchCursor
arcpy.da.UpdateCursor
arcpy.SelectLayerByLocation_management
arcpy.SelectLayerByAttribute_management
layer表示を一旦全部オフにする
はじめに、
arcpy.mp.ArcGISProject projectを取得
arcgisProj.listMaps map一覧を取得
mapObj.listLayers layer一覧を取得
をしています。
その後、取得したlayer一覧から、for layer in layers:で各layerの情報を取り、visibleをFALSEにして、非表示にしています。
後続処理で順次特定のfeatureを選択していくのですが、map上で表示していると、この進行状況が見えるようになる（それはそれで面白いのですが）ので、ここで表示を消してしまって、そこで生じうるタイムロスを少しでも防ぐ意味があります。もちろん事前に手動でオフにしておけばいい話なので、なくても良いです。
改めて気づくことですが、この3つのobjectを使うと、mapを操作できることがわかります。使用例は少ないと思いますが、layer内容や凡例などをコードから逐一変更する事もできます。
field一覧を見る
arcpy.ListFields('DemExtent')でDemExtentという名のteatureclassのfieldをすべて調べます。field object（field名）が取得できるので、追加しようとしているfield（ここではPrefCount）が存在してなければ足すようにします。これも一度手動で行えばすむ部分です。
featureを参照用として１つずつあたる
arcpy.da.SearchCursorではfeatureclassのfeatureを1つずつあたっています。第1引数でfeatureを、第2引数でfield群をリストで指定します。
fieldの中でshape@XYだけが異質ですが、これはfeatureの図形オブジェクトから重心のX,Yをタプル(X.Y)で返すものです。
フィールド値はfor row in cursor:に続くrow[0]などで得ることができます。数字はfield群のリストで指定した順番となり、0から始まります。例えば、Minの値を見たいのであればrow[2]となります。
なお、with構文とともに使うと終了処理を意識しないですみますので、セットで使いましょう。
SearchCursorは読み取り専用で、これだけでは、テーブルを見ているだけで実効性はないですが、この操作は入れ子にすることができ、2つ以上のfeature間で値を参照しあうことなどもできるようになるので、もう少し複雑なことができます。
featureを更新用に１つずつあたる
arcpy.da.UpdateCursorでもfeatureclassのfeatureを1つずつあたっていきますが、SearchCursorと違い、UpdateCursorではデータ更新ができます。ただ引数はSearchCursorと同様、第1引数でfeatureを、第2引数でfield群をリストで指定します。
データ更新の手順自体は少し注意が必要で、まずrow[2]=countと値をセットし、最後にcursor.updateRow(row)とすることで完了します。
シンプルにfield同士を計算しても良いのですが、それではfield演算と何ら変わらないので、この例では少し細かいことをしています。
まず、リストにSHAPE@を入れています。これはfeatureの図形です。
次に featuregeometry=row[0]としてその図形の情報を取ります。属性とはいっても、実際は図形なので形（geomerty）を保持しています。
そのgeomertyと、都道府県の情報（PrefecturePolygon）の2つを、空間検索機能であるarcpy.SelectLayerByLocation_managementに投入し、geomertyから見た、都道府県情報との重なりを把握します。
戻り値はいくつかありますが、ここではシンプルにその合致数を把握しています。それが、count = selectedpolygons[2]の部分です。
選択結果は複数featureになっていることが多いので、これも同じくSearchCursorなどを入れ子にして使うことで、「DemExtentに重なる個々の都道府県境界の…」などのような、より複雑な分岐を持つ処理もできます。
なお、最後の1つ前の処理として、arcpy.SelectLayerByAttribute_managementがありますが、これは最後の処理後に選択状態で残っているfeatureを解除するためのものです。深い意味はありません。
結果
色々書きましたが、結果はこんな感じになります。PrefCountに書き込んだ数字をラベルとして表示しています。長野県は内陸ということもあり非常に多くの県と接していますが、ここでは9（自分の県1を除くと8）となっています。表題サムネイルにもあげてますが、兵庫県も同じく9になります。
画像は標高緯度経度情報の長野県のデータ、背景は地理院地図
標高緯度経度データの作成に当たっては、国土地理院長の承認を得て、同院発行の基盤地図情報数値標高モデルデータを使用しました。（測量法に基づく国土地理院長承認（使用）R3JHs592。）
右端のPrefCountをラベルにして表示
まとめ
今回はfeatureclassのfeatureを１つずつ調べる方法を紹介しました。
前回はファイル単位（featureclass単位）の処理、今回はファイルの中のデータ単位（feature単位）の処理ということで、この2つを知っておけば、大半のことはどうにでもなります。
残るはラスタを1ピクセルごとに操作するという、根気のいるものがありますが。
ただ、これら逐一処理は確実ではあるのもの、どうしても処理時間がかかってしまうので、既存で関数があればそちらを使うほうが良いと思います。
せめてもの抵抗ですが、データをI/Oアクセスの速いメディア、M.2 SSDなどに保存しておいて、処理すると多少早くなります。
「どうにも、うまく自分の要望を叶えられる機能がないな」と思ったときなどに、参考にしていただけると幸いです。
今後の予定としては、１つ戻ってラスタ処理などをpythonで行う事例を紹介しようと思います。
最後に
ここまでご覧いただき、ありがとうございました。
普段は北海道に本拠地を置くNPOに所属し、環境保全を主な題材としてGISやリモセンに関する仕事をしています。
コンサベーションGISコンソーシアムジャパン の活動もその1つです。
この分野の仕事や活動でなにかお困りのある方、ご相談ごとのある方など、是非コメントをお寄せいいただくか、上記WEBサイト掲載のメールアドレスまでメールをお寄せください。
コンサベーションGISコンソーシアムジャパンの活動及びこのnoteでの活動はボランティアで行っているのですが、何分資金が不足しています。
もしサポートしてもいいなという方がいらっしゃいましたら、よろしくお願いします。
この記事が気に入ったらサポートをしてみませんか？