さて、また間が空いてしまいましたが、今回もArcPyレシピを書いてみたいと思います。テーマは「フィーチャクラスの座標系の取得です」。例えば、GDBに複数のフィーチャクラスが入っている際、座標系が異なるフィーチャクラスが入ってしまっていた、というご経験があるのではないでしょうか。今回のレシピはそういったトラブルを回避するために役立ちます。
環境
Windows7 64bit
ArcGIS10.2.0
Python2.7.2
データ
- GDBの作成
- フィーチャクラスの作成
→座標系は右記にしました(平面直角座標系9系(JGD2000)、平面直角座標系8系(JGD2000)、平面直角座標系7系(JGD2000))
サンプルコード
対象のGDBに存在するフィーチャクラスの座標系のコードとフィーチャクラス名を多次元リストで取得するサンプルです。座標系のコードはこのサイトをご参照願います。
GISのための測地成果、測地系、楕円体、投影座標系、EPSGコードのまとめ - 自然環境保全のための周辺技術
# -*- coding: utf-8 -*- import arcpy from os.path import join def get_projection_code(path): """フィーチャクラスの座標系のリストを取得します arguments: path -- 任意フォルダを指定してください。 指定したフォルダにprojection.csvを出力します。 """ arcpy.env.workspace = path projection_list = [] try: for gdb in arcpy.ListWorkspaces(): arcpy.env.workspace = join(path,gdb) if gdb.find('Raster') > -1: featureclasses = arcpy.ListDatasets() else: featureclasses = arcpy.ListFeatureClasses() for fc in featureclasses: projection = arcpy.Describe(fc).spatialReference.factoryCode if projection not in projection_list: projection_list.append([fc,projection]) #ラスタの場合でも、ラスタ座標系ではなくXY座標系にしているので要注意 return projection_list except Exception as e: raise if __name__ == '__main__': lists = get_projection_code(ur"D:\test") for l in lists: print l
結果を見てみると、見た目は文字コードの問題でこのようになっていますが(見た目だけのものなので、特に問題なし)、ばっちり座標系のコードは取得できました。自分はよく座標系の設定をミスってしまうので、こういうチェックとかがあると非常に便利だと思います。簡単ではありますが、本日は以上です!