ArcGISでPythonスクリプト書いてみた。


そろそろ来月に向け学会モードに突入。
まずは、下ごしらえということでバックデータを整備してます。
今回は、気象庁が作っているメッシュ気候図2000を使って温量指数のラスタマップ作り。なんですが販売されているデータには3次メッシュコードと気温のデータが入っているただのテキストデータなのでそのままだと位置情報データとしてまったく使い物になりません。というわけで、変換させないといけないわけですが、今回はその作業にPython使ってみます。

ArcGISとPython

ArcGISはいろいろな言語で拡張できるようになっているわけですが、その言語の一つにPythonがあります。拡張言語として独自の言語を使うGISもある中で、これはかなりえらい。
ただ、日本語の情報は皆無です。ので勉強するのにかなり苦労しました。そのかわりESRIのサイトにドキュメントが整備されています。でも、英語だけです。さすがESRIクオリティ。
それと、使えるメソッドとかオブジェクトをまとめた チートシートみたいなものもESRIが作ってくれてます。こいつはなかなかに便利。使えるメソッドとかがまとめられてるのでこいつを見ると、どういうことができそうとかイメージしやすいです。ただ、propertyなる見慣れない単語が。あ、attributeのことね。typoどころの問題じゃないけど、ま、そこもESRIクオリティ。
それとArcGISにはmodelbuilderという機能があります。これはさせたい作業のフローをつくると、その処理をソースコード吐き出してくれるという機能なんですが、これでPythonスクリプトを出してお勉強するといいかもです。
後、ESRIの海外サイトにも便利なPythonスクリプトがあげられているので、そいつらをいくつか眺めるのも勉強になりました。

変換スクリプト

かなりニッチ、というか僕くらいしか使う人はいなさそうですが、ArcGISをpythonでごにょごにょしたいという人のためにソースコード晒しときます。
というわけで、いきなりソースコードどん。

[python]
import arcgisscripting, csv, os, traceback, sys

##########################################################
##def func
##########################################################
def calc_WI(wd, second_mesh_code, third_mesh_code):
"""
Arguments:
– `wd`:ディレクトリ指定
– `second_mesh_code`:二次メッシュコードの文字列
– `third_mesh_code`:三次メッシュコードの文字列

"""
os.chdir(wd)
filename ="mean_temperature_" + second_mesh_code + ".csv"
tempdata = open (filename)
frag = 0

for line in csv.reader(tempdata):
if line[0] == third_mesh_code:
frag = 1
line = [float(each_temp)-50 for each_temp in line]
line[0] = line[0] + 50
after_thresholds = []
after_thresholds.append(line[0])
for each_after_value in range(1,13):
if line[each_after_value] < 0:
after_thresholds.append(0)
else:
after_thresholds.append(line[each_after_value])
wi =(sum(after_thresholds) – line[0]) / 10
line.append(wi)
return wi
break
## 気温情報がない場合は、99999を返す
if frag == 0:
wi = 99999
return wi

########################################################
##to ArcGIS
############################################################

gp = arcgisscripting.create()
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
gp.OverwriteOutput = 1

inshape = "D:/gis/test/mesh03_jgd_43.shp"
#gp.addfield(inshape, "WI", "long", "9")

rows = gp.UpdateCursor(inshape)
row = rows.Next()

try:
while row:
get_second_mesh_code = str(int(row.CODE/10000))
get_third_mesh_code = str(int(row.CODE))
wi = calc_WI("D:/gis/temp_mesh/",get_second_mesh_code, get_third_mesh_code)
row.WI = wi
gp.AddMessage(get_third_mesh_code)
print get_third_mesh_code
rows.UpdateRow(row)
row = rows.Next()

except:
tb = sys.exc_info()[2]
tbinfo = traceback.format_tb(tb)[0]
pymsg = tbinfo + "n" + str(sys.exc_type)+ ": " + str(sys.exc_value)
gp.AddError(pymsg)
print pymsg
[/python]

前半は温量指数を求めて返す関数の定義。後半はその関数を使ってすでにある3次メッシュポリゴンshpの属性テーブルを、定義した関数使って更新していく処理になります。メッシュポリゴンデータはさんがmurakamiさんが公開されているデータを拝借しました。ありがとうございます!!

属性テーブル更新には、各行をrowsオブジェクトにさせてから、処理を回します。sqlとかをPythonでごにょごにょしたことのある人なら常套手段なあれです。

というわけで、かなりニッチな話でした。ま、よく考えてみると国土数値情報とか標準地域メッシュベースのGISデータって結構ありそうなので、そいつらを変換する時とか参考にできるかもしれません。

それとpythonということで、いろいろ便利なモジュールが使えるかも。rpyとかかませば、クリック一つで解析終了みたいなウマーな仕組み作れるかもしれません。ま、多分Rの方でやっちゃう気がするけど。